Actual source code: baijsolvnat4.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  1:  #include <../src/mat/impls/baij/seq/baij.h>
  2:  #include <petsc/private/kernels/blockinvert.h>

  4: /*
  5:       Special case where the matrix was ILU(0) factored in the natural
  6:    ordering. This eliminates the need for the column and row permutation.
  7: */
  8: PetscErrorCode MatSolve_SeqBAIJ_4_NaturalOrdering_inplace(Mat A,Vec bb,Vec xx)
  9: {
 10:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ*)A->data;
 11:   PetscInt          n  =a->mbs;
 12:   const PetscInt    *ai=a->i,*aj=a->j;
 13:   PetscErrorCode    ierr;
 14:   const PetscInt    *diag = a->diag;
 15:   const MatScalar   *aa   =a->a;
 16:   PetscScalar       *x;
 17:   const PetscScalar *b;

 20:   VecGetArrayRead(bb,&b);
 21:   VecGetArray(xx,&x);

 23: #if defined(PETSC_USE_FORTRAN_KERNEL_SOLVEBAIJ)
 24:   {
 25:     static PetscScalar w[2000]; /* very BAD need to fix */
 26:     fortransolvebaij4_(&n,x,ai,aj,diag,aa,b,w);
 27:   }
 28: #elif defined(PETSC_USE_FORTRAN_KERNEL_SOLVEBAIJUNROLL)
 29:   fortransolvebaij4unroll_(&n,x,ai,aj,diag,aa,b);
 30: #else
 31:   {
 32:     PetscScalar     s1,s2,s3,s4,x1,x2,x3,x4;
 33:     const MatScalar *v;
 34:     PetscInt        jdx,idt,idx,nz,i,ai16;
 35:     const PetscInt  *vi;

 37:     /* forward solve the lower triangular */
 38:     idx  = 0;
 39:     x[0] = b[0]; x[1] = b[1]; x[2] = b[2]; x[3] = b[3];
 40:     for (i=1; i<n; i++) {
 41:       v    =  aa      + 16*ai[i];
 42:       vi   =  aj      + ai[i];
 43:       nz   =  diag[i] - ai[i];
 44:       idx +=  4;
 45:       s1   =  b[idx];s2 = b[1+idx];s3 = b[2+idx];s4 = b[3+idx];
 46:       while (nz--) {
 47:         jdx = 4*(*vi++);
 48:         x1  = x[jdx];x2 = x[1+jdx];x3 = x[2+jdx];x4 = x[3+jdx];
 49:         s1 -= v[0]*x1 + v[4]*x2 + v[8]*x3  + v[12]*x4;
 50:         s2 -= v[1]*x1 + v[5]*x2 + v[9]*x3  + v[13]*x4;
 51:         s3 -= v[2]*x1 + v[6]*x2 + v[10]*x3 + v[14]*x4;
 52:         s4 -= v[3]*x1 + v[7]*x2 + v[11]*x3 + v[15]*x4;
 53:         v  += 16;
 54:       }
 55:       x[idx]   = s1;
 56:       x[1+idx] = s2;
 57:       x[2+idx] = s3;
 58:       x[3+idx] = s4;
 59:     }
 60:     /* backward solve the upper triangular */
 61:     idt = 4*(n-1);
 62:     for (i=n-1; i>=0; i--) {
 63:       ai16 = 16*diag[i];
 64:       v    = aa + ai16 + 16;
 65:       vi   = aj + diag[i] + 1;
 66:       nz   = ai[i+1] - diag[i] - 1;
 67:       s1   = x[idt];  s2 = x[1+idt];
 68:       s3   = x[2+idt];s4 = x[3+idt];
 69:       while (nz--) {
 70:         idx = 4*(*vi++);
 71:         x1  = x[idx];   x2 = x[1+idx];x3    = x[2+idx]; x4 = x[3+idx];
 72:         s1 -= v[0]*x1 + v[4]*x2 + v[8]*x3   + v[12]*x4;
 73:         s2 -= v[1]*x1 + v[5]*x2 + v[9]*x3   + v[13]*x4;
 74:         s3 -= v[2]*x1 + v[6]*x2 + v[10]*x3  + v[14]*x4;
 75:         s4 -= v[3]*x1 + v[7]*x2 + v[11]*x3  + v[15]*x4;
 76:         v  += 16;
 77:       }
 78:       v        = aa + ai16;
 79:       x[idt]   = v[0]*s1 + v[4]*s2 + v[8]*s3  + v[12]*s4;
 80:       x[1+idt] = v[1]*s1 + v[5]*s2 + v[9]*s3  + v[13]*s4;
 81:       x[2+idt] = v[2]*s1 + v[6]*s2 + v[10]*s3 + v[14]*s4;
 82:       x[3+idt] = v[3]*s1 + v[7]*s2 + v[11]*s3 + v[15]*s4;
 83:       idt     -= 4;
 84:     }
 85:   }
 86: #endif

 88:   VecRestoreArrayRead(bb,&b);
 89:   VecRestoreArray(xx,&x);
 90:   PetscLogFlops(2.0*16*(a->nz) - 4.0*A->cmap->n);
 91:   return(0);
 92: }

 94: PetscErrorCode MatSolve_SeqBAIJ_4_NaturalOrdering(Mat A,Vec bb,Vec xx)
 95: {
 96:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ*)A->data;
 97:   const PetscInt    n=a->mbs,*vi,*ai=a->i,*aj=a->j,*adiag=a->diag;
 98:   PetscInt          i,k,nz,idx,jdx,idt;
 99:   PetscErrorCode    ierr;
100:   const PetscInt    bs = A->rmap->bs,bs2 = a->bs2;
101:   const MatScalar   *aa=a->a,*v;
102:   PetscScalar       *x;
103:   const PetscScalar *b;
104:   PetscScalar       s1,s2,s3,s4,x1,x2,x3,x4;

107:   VecGetArrayRead(bb,&b);
108:   VecGetArray(xx,&x);
109:   /* forward solve the lower triangular */
110:   idx  = 0;
111:   x[0] = b[idx]; x[1] = b[1+idx];x[2] = b[2+idx];x[3] = b[3+idx];
112:   for (i=1; i<n; i++) {
113:     v   = aa + bs2*ai[i];
114:     vi  = aj + ai[i];
115:     nz  = ai[i+1] - ai[i];
116:     idx = bs*i;
117:     s1  = b[idx];s2 = b[1+idx];s3 = b[2+idx];s4 = b[3+idx];
118:     for (k=0; k<nz; k++) {
119:       jdx = bs*vi[k];
120:       x1  = x[jdx];x2 = x[1+jdx]; x3 =x[2+jdx];x4 =x[3+jdx];
121:       s1 -= v[0]*x1 + v[4]*x2 + v[8]*x3 + v[12]*x4;
122:       s2 -= v[1]*x1 + v[5]*x2 + v[9]*x3 + v[13]*x4;
123:       s3 -= v[2]*x1 + v[6]*x2 + v[10]*x3 + v[14]*x4;
124:       s4 -= v[3]*x1 + v[7]*x2 + v[11]*x3 + v[15]*x4;

126:       v +=  bs2;
127:     }

129:     x[idx]   = s1;
130:     x[1+idx] = s2;
131:     x[2+idx] = s3;
132:     x[3+idx] = s4;
133:   }

135:   /* backward solve the upper triangular */
136:   for (i=n-1; i>=0; i--) {
137:     v   = aa + bs2*(adiag[i+1]+1);
138:     vi  = aj + adiag[i+1]+1;
139:     nz  = adiag[i] - adiag[i+1]-1;
140:     idt = bs*i;
141:     s1  = x[idt];  s2 = x[1+idt];s3 = x[2+idt];s4 = x[3+idt];

143:     for (k=0; k<nz; k++) {
144:       idx = bs*vi[k];
145:       x1  = x[idx];   x2 = x[1+idx]; x3 = x[2+idx];x4 = x[3+idx];
146:       s1 -= v[0]*x1 + v[4]*x2 + v[8]*x3 + v[12]*x4;
147:       s2 -= v[1]*x1 + v[5]*x2 + v[9]*x3 + v[13]*x4;
148:       s3 -= v[2]*x1 + v[6]*x2 + v[10]*x3 + v[14]*x4;
149:       s4 -= v[3]*x1 + v[7]*x2 + v[11]*x3 + v[15]*x4;

151:       v +=  bs2;
152:     }
153:     /* x = inv_diagonal*x */
154:     x[idt]   = v[0]*s1 + v[4]*s2 + v[8]*s3 + v[12]*s4;
155:     x[1+idt] = v[1]*s1 + v[5]*s2 + v[9]*s3 + v[13]*s4;
156:     x[2+idt] = v[2]*s1 + v[6]*s2 + v[10]*s3 + v[14]*s4;
157:     x[3+idt] = v[3]*s1 + v[7]*s2 + v[11]*s3 + v[15]*s4;

159:   }

161:   VecRestoreArrayRead(bb,&b);
162:   VecRestoreArray(xx,&x);
163:   PetscLogFlops(2.0*bs2*(a->nz) - bs*A->cmap->n);
164:   return(0);
165: }

167: PetscErrorCode MatSolve_SeqBAIJ_4_NaturalOrdering_Demotion(Mat A,Vec bb,Vec xx)
168: {
169:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ*)A->data;
170:   const PetscInt    n  =a->mbs,*ai=a->i,*aj=a->j,*diag=a->diag;
171:   PetscErrorCode    ierr;
172:   const MatScalar   *aa=a->a;
173:   const PetscScalar *b;
174:   PetscScalar       *x;

177:   VecGetArrayRead(bb,&b);
178:   VecGetArray(xx,&x);

180:   {
181:     MatScalar       s1,s2,s3,s4,x1,x2,x3,x4;
182:     const MatScalar *v;
183:     MatScalar       *t=(MatScalar*)x;
184:     PetscInt        jdx,idt,idx,nz,i,ai16;
185:     const PetscInt  *vi;

187:     /* forward solve the lower triangular */
188:     idx  = 0;
189:     t[0] = (MatScalar)b[0];
190:     t[1] = (MatScalar)b[1];
191:     t[2] = (MatScalar)b[2];
192:     t[3] = (MatScalar)b[3];
193:     for (i=1; i<n; i++) {
194:       v    =  aa      + 16*ai[i];
195:       vi   =  aj      + ai[i];
196:       nz   =  diag[i] - ai[i];
197:       idx +=  4;
198:       s1   = (MatScalar)b[idx];
199:       s2   = (MatScalar)b[1+idx];
200:       s3   = (MatScalar)b[2+idx];
201:       s4   = (MatScalar)b[3+idx];
202:       while (nz--) {
203:         jdx = 4*(*vi++);
204:         x1  = t[jdx];
205:         x2  = t[1+jdx];
206:         x3  = t[2+jdx];
207:         x4  = t[3+jdx];
208:         s1 -= v[0]*x1 + v[4]*x2 + v[8]*x3  + v[12]*x4;
209:         s2 -= v[1]*x1 + v[5]*x2 + v[9]*x3  + v[13]*x4;
210:         s3 -= v[2]*x1 + v[6]*x2 + v[10]*x3 + v[14]*x4;
211:         s4 -= v[3]*x1 + v[7]*x2 + v[11]*x3 + v[15]*x4;
212:         v  += 16;
213:       }
214:       t[idx]   = s1;
215:       t[1+idx] = s2;
216:       t[2+idx] = s3;
217:       t[3+idx] = s4;
218:     }
219:     /* backward solve the upper triangular */
220:     idt = 4*(n-1);
221:     for (i=n-1; i>=0; i--) {
222:       ai16 = 16*diag[i];
223:       v    = aa + ai16 + 16;
224:       vi   = aj + diag[i] + 1;
225:       nz   = ai[i+1] - diag[i] - 1;
226:       s1   = t[idt];
227:       s2   = t[1+idt];
228:       s3   = t[2+idt];
229:       s4   = t[3+idt];
230:       while (nz--) {
231:         idx = 4*(*vi++);
232:         x1  = (MatScalar)x[idx];
233:         x2  = (MatScalar)x[1+idx];
234:         x3  = (MatScalar)x[2+idx];
235:         x4  = (MatScalar)x[3+idx];
236:         s1 -= v[0]*x1 + v[4]*x2 + v[8]*x3  + v[12]*x4;
237:         s2 -= v[1]*x1 + v[5]*x2 + v[9]*x3  + v[13]*x4;
238:         s3 -= v[2]*x1 + v[6]*x2 + v[10]*x3 + v[14]*x4;
239:         s4 -= v[3]*x1 + v[7]*x2 + v[11]*x3 + v[15]*x4;
240:         v  += 16;
241:       }
242:       v        = aa + ai16;
243:       x[idt]   = (PetscScalar)(v[0]*s1 + v[4]*s2 + v[8]*s3  + v[12]*s4);
244:       x[1+idt] = (PetscScalar)(v[1]*s1 + v[5]*s2 + v[9]*s3  + v[13]*s4);
245:       x[2+idt] = (PetscScalar)(v[2]*s1 + v[6]*s2 + v[10]*s3 + v[14]*s4);
246:       x[3+idt] = (PetscScalar)(v[3]*s1 + v[7]*s2 + v[11]*s3 + v[15]*s4);
247:       idt     -= 4;
248:     }
249:   }

251:   VecRestoreArrayRead(bb,&b);
252:   VecRestoreArray(xx,&x);
253:   PetscLogFlops(2.0*16*(a->nz) - 4.0*A->cmap->n);
254:   return(0);
255: }

257: #if defined(PETSC_HAVE_SSE)

259: #include PETSC_HAVE_SSE
260: PetscErrorCode MatSolve_SeqBAIJ_4_NaturalOrdering_SSE_Demotion_usj(Mat A,Vec bb,Vec xx)
261: {
262:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
263:   unsigned short *aj=(unsigned short*)a->j;
265:   int            *ai=a->i,n=a->mbs,*diag = a->diag;
266:   MatScalar      *aa=a->a;
267:   PetscScalar    *x,*b;

270:   SSE_SCOPE_BEGIN;
271:   /*
272:      Note: This code currently uses demotion of double
273:      to float when performing the mixed-mode computation.
274:      This may not be numerically reasonable for all applications.
275:   */
276:   PREFETCH_NTA(aa+16*ai[1]);

278:   VecGetArray(bb,&b);
279:   VecGetArray(xx,&x);
280:   {
281:     /* x will first be computed in single precision then promoted inplace to double */
282:     MatScalar      *v,*t=(MatScalar*)x;
283:     int            nz,i,idt,ai16;
284:     unsigned int   jdx,idx;
285:     unsigned short *vi;
286:     /* Forward solve the lower triangular factor. */

288:     /* First block is the identity. */
289:     idx = 0;
290:     CONVERT_DOUBLE4_FLOAT4(t,b);
291:     v =  aa + 16*((unsigned int)ai[1]);

293:     for (i=1; i<n; ) {
294:       PREFETCH_NTA(&v[8]);
295:       vi   =  aj      + ai[i];
296:       nz   =  diag[i] - ai[i];
297:       idx +=  4;

299:       /* Demote RHS from double to float. */
300:       CONVERT_DOUBLE4_FLOAT4(&t[idx],&b[idx]);
301:       LOAD_PS(&t[idx],XMM7);

303:       while (nz--) {
304:         PREFETCH_NTA(&v[16]);
305:         jdx = 4*((unsigned int)(*vi++));

307:         /* 4x4 Matrix-Vector product with negative accumulation: */
308:         SSE_INLINE_BEGIN_2(&t[jdx],v)
309:         SSE_LOAD_PS(SSE_ARG_1,FLOAT_0,XMM6)

311:         /* First Column */
312:         SSE_COPY_PS(XMM0,XMM6)
313:         SSE_SHUFFLE(XMM0,XMM0,0x00)
314:         SSE_MULT_PS_M(XMM0,SSE_ARG_2,FLOAT_0)
315:         SSE_SUB_PS(XMM7,XMM0)

317:         /* Second Column */
318:         SSE_COPY_PS(XMM1,XMM6)
319:         SSE_SHUFFLE(XMM1,XMM1,0x55)
320:         SSE_MULT_PS_M(XMM1,SSE_ARG_2,FLOAT_4)
321:         SSE_SUB_PS(XMM7,XMM1)

323:         SSE_PREFETCH_NTA(SSE_ARG_2,FLOAT_24)

325:         /* Third Column */
326:         SSE_COPY_PS(XMM2,XMM6)
327:         SSE_SHUFFLE(XMM2,XMM2,0xAA)
328:         SSE_MULT_PS_M(XMM2,SSE_ARG_2,FLOAT_8)
329:         SSE_SUB_PS(XMM7,XMM2)

331:         /* Fourth Column */
332:         SSE_COPY_PS(XMM3,XMM6)
333:         SSE_SHUFFLE(XMM3,XMM3,0xFF)
334:         SSE_MULT_PS_M(XMM3,SSE_ARG_2,FLOAT_12)
335:         SSE_SUB_PS(XMM7,XMM3)
336:         SSE_INLINE_END_2

338:         v += 16;
339:       }
340:       v =  aa + 16*ai[++i];
341:       PREFETCH_NTA(v);
342:       STORE_PS(&t[idx],XMM7);
343:     }

345:     /* Backward solve the upper triangular factor.*/

347:     idt  = 4*(n-1);
348:     ai16 = 16*diag[n-1];
349:     v    = aa + ai16 + 16;
350:     for (i=n-1; i>=0; ) {
351:       PREFETCH_NTA(&v[8]);
352:       vi = aj + diag[i] + 1;
353:       nz = ai[i+1] - diag[i] - 1;

355:       LOAD_PS(&t[idt],XMM7);

357:       while (nz--) {
358:         PREFETCH_NTA(&v[16]);
359:         idx = 4*((unsigned int)(*vi++));

361:         /* 4x4 Matrix-Vector Product with negative accumulation: */
362:         SSE_INLINE_BEGIN_2(&t[idx],v)
363:         SSE_LOAD_PS(SSE_ARG_1,FLOAT_0,XMM6)

365:         /* First Column */
366:         SSE_COPY_PS(XMM0,XMM6)
367:         SSE_SHUFFLE(XMM0,XMM0,0x00)
368:         SSE_MULT_PS_M(XMM0,SSE_ARG_2,FLOAT_0)
369:         SSE_SUB_PS(XMM7,XMM0)

371:         /* Second Column */
372:         SSE_COPY_PS(XMM1,XMM6)
373:         SSE_SHUFFLE(XMM1,XMM1,0x55)
374:         SSE_MULT_PS_M(XMM1,SSE_ARG_2,FLOAT_4)
375:         SSE_SUB_PS(XMM7,XMM1)

377:         SSE_PREFETCH_NTA(SSE_ARG_2,FLOAT_24)

379:         /* Third Column */
380:         SSE_COPY_PS(XMM2,XMM6)
381:         SSE_SHUFFLE(XMM2,XMM2,0xAA)
382:         SSE_MULT_PS_M(XMM2,SSE_ARG_2,FLOAT_8)
383:         SSE_SUB_PS(XMM7,XMM2)

385:         /* Fourth Column */
386:         SSE_COPY_PS(XMM3,XMM6)
387:         SSE_SHUFFLE(XMM3,XMM3,0xFF)
388:         SSE_MULT_PS_M(XMM3,SSE_ARG_2,FLOAT_12)
389:         SSE_SUB_PS(XMM7,XMM3)
390:         SSE_INLINE_END_2
391:         v += 16;
392:       }
393:       v    = aa + ai16;
394:       ai16 = 16*diag[--i];
395:       PREFETCH_NTA(aa+ai16+16);
396:       /*
397:          Scale the result by the diagonal 4x4 block,
398:          which was inverted as part of the factorization
399:       */
400:       SSE_INLINE_BEGIN_3(v,&t[idt],aa+ai16)
401:       /* First Column */
402:       SSE_COPY_PS(XMM0,XMM7)
403:       SSE_SHUFFLE(XMM0,XMM0,0x00)
404:       SSE_MULT_PS_M(XMM0,SSE_ARG_1,FLOAT_0)

406:       /* Second Column */
407:       SSE_COPY_PS(XMM1,XMM7)
408:       SSE_SHUFFLE(XMM1,XMM1,0x55)
409:       SSE_MULT_PS_M(XMM1,SSE_ARG_1,FLOAT_4)
410:       SSE_ADD_PS(XMM0,XMM1)

412:       SSE_PREFETCH_NTA(SSE_ARG_3,FLOAT_24)

414:       /* Third Column */
415:       SSE_COPY_PS(XMM2,XMM7)
416:       SSE_SHUFFLE(XMM2,XMM2,0xAA)
417:       SSE_MULT_PS_M(XMM2,SSE_ARG_1,FLOAT_8)
418:       SSE_ADD_PS(XMM0,XMM2)

420:       /* Fourth Column */
421:       SSE_COPY_PS(XMM3,XMM7)
422:       SSE_SHUFFLE(XMM3,XMM3,0xFF)
423:       SSE_MULT_PS_M(XMM3,SSE_ARG_1,FLOAT_12)
424:       SSE_ADD_PS(XMM0,XMM3)

426:       SSE_STORE_PS(SSE_ARG_2,FLOAT_0,XMM0)
427:       SSE_INLINE_END_3

429:       v    = aa + ai16 + 16;
430:       idt -= 4;
431:     }

433:     /* Convert t from single precision back to double precision (inplace)*/
434:     idt = 4*(n-1);
435:     for (i=n-1; i>=0; i--) {
436:       /*     CONVERT_FLOAT4_DOUBLE4(&x[idt],&t[idt]); */
437:       /* Unfortunately, CONVERT_ will count from 0 to 3 which doesn't work here. */
438:       PetscScalar *xtemp=&x[idt];
439:       MatScalar   *ttemp=&t[idt];
440:       xtemp[3] = (PetscScalar)ttemp[3];
441:       xtemp[2] = (PetscScalar)ttemp[2];
442:       xtemp[1] = (PetscScalar)ttemp[1];
443:       xtemp[0] = (PetscScalar)ttemp[0];
444:       idt     -= 4;
445:     }

447:   } /* End of artificial scope. */
448:   VecRestoreArray(bb,&b);
449:   VecRestoreArray(xx,&x);
450:   PetscLogFlops(2.0*16*(a->nz) - 4.0*A->cmap->n);
451:   SSE_SCOPE_END;
452:   return(0);
453: }

455: PetscErrorCode MatSolve_SeqBAIJ_4_NaturalOrdering_SSE_Demotion(Mat A,Vec bb,Vec xx)
456: {
457:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
458:   int            *aj=a->j;
460:   int            *ai=a->i,n=a->mbs,*diag = a->diag;
461:   MatScalar      *aa=a->a;
462:   PetscScalar    *x,*b;

465:   SSE_SCOPE_BEGIN;
466:   /*
467:      Note: This code currently uses demotion of double
468:      to float when performing the mixed-mode computation.
469:      This may not be numerically reasonable for all applications.
470:   */
471:   PREFETCH_NTA(aa+16*ai[1]);

473:   VecGetArray(bb,&b);
474:   VecGetArray(xx,&x);
475:   {
476:     /* x will first be computed in single precision then promoted inplace to double */
477:     MatScalar *v,*t=(MatScalar*)x;
478:     int       nz,i,idt,ai16;
479:     int       jdx,idx;
480:     int       *vi;
481:     /* Forward solve the lower triangular factor. */

483:     /* First block is the identity. */
484:     idx = 0;
485:     CONVERT_DOUBLE4_FLOAT4(t,b);
486:     v =  aa + 16*ai[1];

488:     for (i=1; i<n; ) {
489:       PREFETCH_NTA(&v[8]);
490:       vi   =  aj      + ai[i];
491:       nz   =  diag[i] - ai[i];
492:       idx +=  4;

494:       /* Demote RHS from double to float. */
495:       CONVERT_DOUBLE4_FLOAT4(&t[idx],&b[idx]);
496:       LOAD_PS(&t[idx],XMM7);

498:       while (nz--) {
499:         PREFETCH_NTA(&v[16]);
500:         jdx = 4*(*vi++);
501: /*          jdx = *vi++; */

503:         /* 4x4 Matrix-Vector product with negative accumulation: */
504:         SSE_INLINE_BEGIN_2(&t[jdx],v)
505:         SSE_LOAD_PS(SSE_ARG_1,FLOAT_0,XMM6)

507:         /* First Column */
508:         SSE_COPY_PS(XMM0,XMM6)
509:         SSE_SHUFFLE(XMM0,XMM0,0x00)
510:         SSE_MULT_PS_M(XMM0,SSE_ARG_2,FLOAT_0)
511:         SSE_SUB_PS(XMM7,XMM0)

513:         /* Second Column */
514:         SSE_COPY_PS(XMM1,XMM6)
515:         SSE_SHUFFLE(XMM1,XMM1,0x55)
516:         SSE_MULT_PS_M(XMM1,SSE_ARG_2,FLOAT_4)
517:         SSE_SUB_PS(XMM7,XMM1)

519:         SSE_PREFETCH_NTA(SSE_ARG_2,FLOAT_24)

521:         /* Third Column */
522:         SSE_COPY_PS(XMM2,XMM6)
523:         SSE_SHUFFLE(XMM2,XMM2,0xAA)
524:         SSE_MULT_PS_M(XMM2,SSE_ARG_2,FLOAT_8)
525:         SSE_SUB_PS(XMM7,XMM2)

527:         /* Fourth Column */
528:         SSE_COPY_PS(XMM3,XMM6)
529:         SSE_SHUFFLE(XMM3,XMM3,0xFF)
530:         SSE_MULT_PS_M(XMM3,SSE_ARG_2,FLOAT_12)
531:         SSE_SUB_PS(XMM7,XMM3)
532:         SSE_INLINE_END_2

534:         v += 16;
535:       }
536:       v =  aa + 16*ai[++i];
537:       PREFETCH_NTA(v);
538:       STORE_PS(&t[idx],XMM7);
539:     }

541:     /* Backward solve the upper triangular factor.*/

543:     idt  = 4*(n-1);
544:     ai16 = 16*diag[n-1];
545:     v    = aa + ai16 + 16;
546:     for (i=n-1; i>=0; ) {
547:       PREFETCH_NTA(&v[8]);
548:       vi = aj + diag[i] + 1;
549:       nz = ai[i+1] - diag[i] - 1;

551:       LOAD_PS(&t[idt],XMM7);

553:       while (nz--) {
554:         PREFETCH_NTA(&v[16]);
555:         idx = 4*(*vi++);
556: /*          idx = *vi++; */

558:         /* 4x4 Matrix-Vector Product with negative accumulation: */
559:         SSE_INLINE_BEGIN_2(&t[idx],v)
560:         SSE_LOAD_PS(SSE_ARG_1,FLOAT_0,XMM6)

562:         /* First Column */
563:         SSE_COPY_PS(XMM0,XMM6)
564:         SSE_SHUFFLE(XMM0,XMM0,0x00)
565:         SSE_MULT_PS_M(XMM0,SSE_ARG_2,FLOAT_0)
566:         SSE_SUB_PS(XMM7,XMM0)

568:         /* Second Column */
569:         SSE_COPY_PS(XMM1,XMM6)
570:         SSE_SHUFFLE(XMM1,XMM1,0x55)
571:         SSE_MULT_PS_M(XMM1,SSE_ARG_2,FLOAT_4)
572:         SSE_SUB_PS(XMM7,XMM1)

574:         SSE_PREFETCH_NTA(SSE_ARG_2,FLOAT_24)

576:         /* Third Column */
577:         SSE_COPY_PS(XMM2,XMM6)
578:         SSE_SHUFFLE(XMM2,XMM2,0xAA)
579:         SSE_MULT_PS_M(XMM2,SSE_ARG_2,FLOAT_8)
580:         SSE_SUB_PS(XMM7,XMM2)

582:         /* Fourth Column */
583:         SSE_COPY_PS(XMM3,XMM6)
584:         SSE_SHUFFLE(XMM3,XMM3,0xFF)
585:         SSE_MULT_PS_M(XMM3,SSE_ARG_2,FLOAT_12)
586:         SSE_SUB_PS(XMM7,XMM3)
587:         SSE_INLINE_END_2
588:         v += 16;
589:       }
590:       v    = aa + ai16;
591:       ai16 = 16*diag[--i];
592:       PREFETCH_NTA(aa+ai16+16);
593:       /*
594:          Scale the result by the diagonal 4x4 block,
595:          which was inverted as part of the factorization
596:       */
597:       SSE_INLINE_BEGIN_3(v,&t[idt],aa+ai16)
598:       /* First Column */
599:       SSE_COPY_PS(XMM0,XMM7)
600:       SSE_SHUFFLE(XMM0,XMM0,0x00)
601:       SSE_MULT_PS_M(XMM0,SSE_ARG_1,FLOAT_0)

603:       /* Second Column */
604:       SSE_COPY_PS(XMM1,XMM7)
605:       SSE_SHUFFLE(XMM1,XMM1,0x55)
606:       SSE_MULT_PS_M(XMM1,SSE_ARG_1,FLOAT_4)
607:       SSE_ADD_PS(XMM0,XMM1)

609:       SSE_PREFETCH_NTA(SSE_ARG_3,FLOAT_24)

611:       /* Third Column */
612:       SSE_COPY_PS(XMM2,XMM7)
613:       SSE_SHUFFLE(XMM2,XMM2,0xAA)
614:       SSE_MULT_PS_M(XMM2,SSE_ARG_1,FLOAT_8)
615:       SSE_ADD_PS(XMM0,XMM2)

617:       /* Fourth Column */
618:       SSE_COPY_PS(XMM3,XMM7)
619:       SSE_SHUFFLE(XMM3,XMM3,0xFF)
620:       SSE_MULT_PS_M(XMM3,SSE_ARG_1,FLOAT_12)
621:       SSE_ADD_PS(XMM0,XMM3)

623:       SSE_STORE_PS(SSE_ARG_2,FLOAT_0,XMM0)
624:       SSE_INLINE_END_3

626:       v    = aa + ai16 + 16;
627:       idt -= 4;
628:     }

630:     /* Convert t from single precision back to double precision (inplace)*/
631:     idt = 4*(n-1);
632:     for (i=n-1; i>=0; i--) {
633:       /*     CONVERT_FLOAT4_DOUBLE4(&x[idt],&t[idt]); */
634:       /* Unfortunately, CONVERT_ will count from 0 to 3 which doesn't work here. */
635:       PetscScalar *xtemp=&x[idt];
636:       MatScalar   *ttemp=&t[idt];
637:       xtemp[3] = (PetscScalar)ttemp[3];
638:       xtemp[2] = (PetscScalar)ttemp[2];
639:       xtemp[1] = (PetscScalar)ttemp[1];
640:       xtemp[0] = (PetscScalar)ttemp[0];
641:       idt     -= 4;
642:     }

644:   } /* End of artificial scope. */
645:   VecRestoreArray(bb,&b);
646:   VecRestoreArray(xx,&x);
647:   PetscLogFlops(2.0*16*(a->nz) - 4.0*A->cmap->n);
648:   SSE_SCOPE_END;
649:   return(0);
650: }

652: #endif