Actual source code: baijfact.c


  2: /*
  3:     Factorization code for BAIJ format.
  4: */
  5: #include <../src/mat/impls/baij/seq/baij.h>
  6: #include <petsc/private/kernels/blockinvert.h>

  8: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2(Mat B,Mat A,const MatFactorInfo *info)
  9: {
 10:   Mat            C     =B;
 11:   Mat_SeqBAIJ    *a    =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data;
 12:   IS             isrow = b->row,isicol = b->icol;
 14:   const PetscInt *r,*ic;
 15:   PetscInt       i,j,k,nz,nzL,row,*pj;
 16:   const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,bs2=a->bs2;
 17:   const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag;
 18:   MatScalar      *rtmp,*pc,*mwork,*pv;
 19:   MatScalar      *aa=a->a,*v;
 20:   PetscInt       flg;
 21:   PetscReal      shift = info->shiftamount;
 22:   PetscBool      allowzeropivot,zeropivotdetected;

 25:   ISGetIndices(isrow,&r);
 26:   ISGetIndices(isicol,&ic);
 27:   allowzeropivot = PetscNot(A->erroriffailure);

 29:   /* generate work space needed by the factorization */
 30:   PetscMalloc2(bs2*n,&rtmp,bs2,&mwork);
 31:   PetscArrayzero(rtmp,bs2*n);

 33:   for (i=0; i<n; i++) {
 34:     /* zero rtmp */
 35:     /* L part */
 36:     nz    = bi[i+1] - bi[i];
 37:     bjtmp = bj + bi[i];
 38:     for  (j=0; j<nz; j++) {
 39:       PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);
 40:     }

 42:     /* U part */
 43:     nz    = bdiag[i] - bdiag[i+1];
 44:     bjtmp = bj + bdiag[i+1]+1;
 45:     for  (j=0; j<nz; j++) {
 46:       PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);
 47:     }

 49:     /* load in initial (unfactored row) */
 50:     nz    = ai[r[i]+1] - ai[r[i]];
 51:     ajtmp = aj + ai[r[i]];
 52:     v     = aa + bs2*ai[r[i]];
 53:     for (j=0; j<nz; j++) {
 54:       PetscArraycpy(rtmp+bs2*ic[ajtmp[j]],v+bs2*j,bs2);
 55:     }

 57:     /* elimination */
 58:     bjtmp = bj + bi[i];
 59:     nzL   = bi[i+1] - bi[i];
 60:     for (k=0; k < nzL; k++) {
 61:       row = bjtmp[k];
 62:       pc  = rtmp + bs2*row;
 63:       for (flg=0,j=0; j<bs2; j++) {
 64:         if (pc[j] != (PetscScalar)0.0) {
 65:           flg = 1;
 66:           break;
 67:         }
 68:       }
 69:       if (flg) {
 70:         pv = b->a + bs2*bdiag[row];
 71:         /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
 72:         PetscKernel_A_gets_A_times_B_2(pc,pv,mwork);

 74:         pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
 75:         pv = b->a + bs2*(bdiag[row+1]+1);
 76:         nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */
 77:         for (j=0; j<nz; j++) {
 78:           /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
 79:           /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
 80:           v    = rtmp + 4*pj[j];
 81:           PetscKernel_A_gets_A_minus_B_times_C_2(v,pc,pv);
 82:           pv  += 4;
 83:         }
 84:         PetscLogFlops(16.0*nz+12); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
 85:       }
 86:     }

 88:     /* finished row so stick it into b->a */
 89:     /* L part */
 90:     pv = b->a + bs2*bi[i];
 91:     pj = b->j + bi[i];
 92:     nz = bi[i+1] - bi[i];
 93:     for (j=0; j<nz; j++) {
 94:       PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);
 95:     }

 97:     /* Mark diagonal and invert diagonal for simpler triangular solves */
 98:     pv   = b->a + bs2*bdiag[i];
 99:     pj   = b->j + bdiag[i];
100:     PetscArraycpy(pv,rtmp+bs2*pj[0],bs2);
101:     PetscKernel_A_gets_inverse_A_2(pv,shift,allowzeropivot,&zeropivotdetected);
102:     if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

104:     /* U part */
105:     pv = b->a + bs2*(bdiag[i+1]+1);
106:     pj = b->j + bdiag[i+1]+1;
107:     nz = bdiag[i] - bdiag[i+1] - 1;
108:     for (j=0; j<nz; j++) {
109:       PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);
110:     }
111:   }

113:   PetscFree2(rtmp,mwork);
114:   ISRestoreIndices(isicol,&ic);
115:   ISRestoreIndices(isrow,&r);

117:   C->ops->solve          = MatSolve_SeqBAIJ_2;
118:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2;
119:   C->assembled           = PETSC_TRUE;

121:   PetscLogFlops(1.333333333333*2*2*2*n); /* from inverting diagonal blocks */
122:   return(0);
123: }

125: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info)
126: {
127:   Mat            C =B;
128:   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data;
130:   PetscInt       i,j,k,nz,nzL,row,*pj;
131:   const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,bs2=a->bs2;
132:   const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag;
133:   MatScalar      *rtmp,*pc,*mwork,*pv;
134:   MatScalar      *aa=a->a,*v;
135:   PetscInt       flg;
136:   PetscReal      shift = info->shiftamount;
137:   PetscBool      allowzeropivot,zeropivotdetected;

140:   allowzeropivot = PetscNot(A->erroriffailure);

142:   /* generate work space needed by the factorization */
143:   PetscMalloc2(bs2*n,&rtmp,bs2,&mwork);
144:   PetscArrayzero(rtmp,bs2*n);

146:   for (i=0; i<n; i++) {
147:     /* zero rtmp */
148:     /* L part */
149:     nz    = bi[i+1] - bi[i];
150:     bjtmp = bj + bi[i];
151:     for  (j=0; j<nz; j++) {
152:       PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);
153:     }

155:     /* U part */
156:     nz    = bdiag[i] - bdiag[i+1];
157:     bjtmp = bj + bdiag[i+1]+1;
158:     for  (j=0; j<nz; j++) {
159:       PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);
160:     }

162:     /* load in initial (unfactored row) */
163:     nz    = ai[i+1] - ai[i];
164:     ajtmp = aj + ai[i];
165:     v     = aa + bs2*ai[i];
166:     for (j=0; j<nz; j++) {
167:       PetscArraycpy(rtmp+bs2*ajtmp[j],v+bs2*j,bs2);
168:     }

170:     /* elimination */
171:     bjtmp = bj + bi[i];
172:     nzL   = bi[i+1] - bi[i];
173:     for (k=0; k < nzL; k++) {
174:       row = bjtmp[k];
175:       pc  = rtmp + bs2*row;
176:       for (flg=0,j=0; j<bs2; j++) {
177:         if (pc[j]!=(PetscScalar)0.0) {
178:           flg = 1;
179:           break;
180:         }
181:       }
182:       if (flg) {
183:         pv = b->a + bs2*bdiag[row];
184:         /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
185:         PetscKernel_A_gets_A_times_B_2(pc,pv,mwork);

187:         pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
188:         pv = b->a + bs2*(bdiag[row+1]+1);
189:         nz = bdiag[row]-bdiag[row+1] - 1; /* num of entries in U(row,:) excluding diag */
190:         for (j=0; j<nz; j++) {
191:           /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
192:           /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
193:           v    = rtmp + 4*pj[j];
194:           PetscKernel_A_gets_A_minus_B_times_C_2(v,pc,pv);
195:           pv  += 4;
196:         }
197:         PetscLogFlops(16.0*nz+12); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
198:       }
199:     }

201:     /* finished row so stick it into b->a */
202:     /* L part */
203:     pv = b->a + bs2*bi[i];
204:     pj = b->j + bi[i];
205:     nz = bi[i+1] - bi[i];
206:     for (j=0; j<nz; j++) {
207:       PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);
208:     }

210:     /* Mark diagonal and invert diagonal for simpler triangular solves */
211:     pv   = b->a + bs2*bdiag[i];
212:     pj   = b->j + bdiag[i];
213:     PetscArraycpy(pv,rtmp+bs2*pj[0],bs2);
214:     PetscKernel_A_gets_inverse_A_2(pv,shift,allowzeropivot,&zeropivotdetected);
215:     if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

217:     /* U part */
218:     /*
219:     pv = b->a + bs2*bi[2*n-i];
220:     pj = b->j + bi[2*n-i];
221:     nz = bi[2*n-i+1] - bi[2*n-i] - 1;
222:     */
223:     pv = b->a + bs2*(bdiag[i+1]+1);
224:     pj = b->j + bdiag[i+1]+1;
225:     nz = bdiag[i] - bdiag[i+1] - 1;
226:     for (j=0; j<nz; j++) {
227:       PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);
228:     }
229:   }
230:   PetscFree2(rtmp,mwork);

232:   C->ops->solve          = MatSolve_SeqBAIJ_2_NaturalOrdering;
233:   C->ops->forwardsolve   = MatForwardSolve_SeqBAIJ_2_NaturalOrdering;
234:   C->ops->backwardsolve  = MatBackwardSolve_SeqBAIJ_2_NaturalOrdering;
235:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering;
236:   C->assembled           = PETSC_TRUE;

238:   PetscLogFlops(1.333333333333*2*2*2*n); /* from inverting diagonal blocks */
239:   return(0);
240: }

242: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_inplace(Mat B,Mat A,const MatFactorInfo *info)
243: {
244:   Mat            C     = B;
245:   Mat_SeqBAIJ    *a    = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data;
246:   IS             isrow = b->row,isicol = b->icol;
248:   const PetscInt *r,*ic;
249:   PetscInt       i,j,n = a->mbs,*bi = b->i,*bj = b->j;
250:   PetscInt       *ajtmpold,*ajtmp,nz,row;
251:   PetscInt       *diag_offset=b->diag,idx,*ai=a->i,*aj=a->j,*pj;
252:   MatScalar      *pv,*v,*rtmp,m1,m2,m3,m4,*pc,*w,*x,x1,x2,x3,x4;
253:   MatScalar      p1,p2,p3,p4;
254:   MatScalar      *ba   = b->a,*aa = a->a;
255:   PetscReal      shift = info->shiftamount;
256:   PetscBool      allowzeropivot,zeropivotdetected;

259:   allowzeropivot = PetscNot(A->erroriffailure);
260:   ISGetIndices(isrow,&r);
261:   ISGetIndices(isicol,&ic);
262:   PetscMalloc1(4*(n+1),&rtmp);

264:   for (i=0; i<n; i++) {
265:     nz    = bi[i+1] - bi[i];
266:     ajtmp = bj + bi[i];
267:     for  (j=0; j<nz; j++) {
268:       x = rtmp+4*ajtmp[j]; x[0] = x[1] = x[2] = x[3] = 0.0;
269:     }
270:     /* load in initial (unfactored row) */
271:     idx      = r[i];
272:     nz       = ai[idx+1] - ai[idx];
273:     ajtmpold = aj + ai[idx];
274:     v        = aa + 4*ai[idx];
275:     for (j=0; j<nz; j++) {
276:       x    = rtmp+4*ic[ajtmpold[j]];
277:       x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3];
278:       v   += 4;
279:     }
280:     row = *ajtmp++;
281:     while (row < i) {
282:       pc = rtmp + 4*row;
283:       p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3];
284:       if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) {
285:         pv    = ba + 4*diag_offset[row];
286:         pj    = bj + diag_offset[row] + 1;
287:         x1    = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3];
288:         pc[0] = m1 = p1*x1 + p3*x2;
289:         pc[1] = m2 = p2*x1 + p4*x2;
290:         pc[2] = m3 = p1*x3 + p3*x4;
291:         pc[3] = m4 = p2*x3 + p4*x4;
292:         nz    = bi[row+1] - diag_offset[row] - 1;
293:         pv   += 4;
294:         for (j=0; j<nz; j++) {
295:           x1    = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3];
296:           x     = rtmp + 4*pj[j];
297:           x[0] -= m1*x1 + m3*x2;
298:           x[1] -= m2*x1 + m4*x2;
299:           x[2] -= m1*x3 + m3*x4;
300:           x[3] -= m2*x3 + m4*x4;
301:           pv   += 4;
302:         }
303:         PetscLogFlops(16.0*nz+12.0);
304:       }
305:       row = *ajtmp++;
306:     }
307:     /* finished row so stick it into b->a */
308:     pv = ba + 4*bi[i];
309:     pj = bj + bi[i];
310:     nz = bi[i+1] - bi[i];
311:     for (j=0; j<nz; j++) {
312:       x     = rtmp+4*pj[j];
313:       pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3];
314:       pv   += 4;
315:     }
316:     /* invert diagonal block */
317:     w    = ba + 4*diag_offset[i];
318:     PetscKernel_A_gets_inverse_A_2(w,shift,allowzeropivot,&zeropivotdetected);
319:     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
320:   }

322:   PetscFree(rtmp);
323:   ISRestoreIndices(isicol,&ic);
324:   ISRestoreIndices(isrow,&r);

326:   C->ops->solve          = MatSolve_SeqBAIJ_2_inplace;
327:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_inplace;
328:   C->assembled           = PETSC_TRUE;

330:   PetscLogFlops(1.333333333333*8*b->mbs); /* from inverting diagonal blocks */
331:   return(0);
332: }
333: /*
334:       Version for when blocks are 2 by 2 Using natural ordering
335: */
336: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo *info)
337: {
338:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data;
340:   PetscInt       i,j,n = a->mbs,*bi = b->i,*bj = b->j;
341:   PetscInt       *ajtmpold,*ajtmp,nz,row;
342:   PetscInt       *diag_offset = b->diag,*ai=a->i,*aj=a->j,*pj;
343:   MatScalar      *pv,*v,*rtmp,*pc,*w,*x;
344:   MatScalar      p1,p2,p3,p4,m1,m2,m3,m4,x1,x2,x3,x4;
345:   MatScalar      *ba   = b->a,*aa = a->a;
346:   PetscReal      shift = info->shiftamount;
347:   PetscBool      allowzeropivot,zeropivotdetected;

350:   allowzeropivot = PetscNot(A->erroriffailure);
351:   PetscMalloc1(4*(n+1),&rtmp);
352:   for (i=0; i<n; i++) {
353:     nz    = bi[i+1] - bi[i];
354:     ajtmp = bj + bi[i];
355:     for  (j=0; j<nz; j++) {
356:       x    = rtmp+4*ajtmp[j];
357:       x[0] = x[1]  = x[2]  = x[3]  = 0.0;
358:     }
359:     /* load in initial (unfactored row) */
360:     nz       = ai[i+1] - ai[i];
361:     ajtmpold = aj + ai[i];
362:     v        = aa + 4*ai[i];
363:     for (j=0; j<nz; j++) {
364:       x    = rtmp+4*ajtmpold[j];
365:       x[0] = v[0];  x[1]  = v[1];  x[2]  = v[2];  x[3]  = v[3];
366:       v   += 4;
367:     }
368:     row = *ajtmp++;
369:     while (row < i) {
370:       pc = rtmp + 4*row;
371:       p1 = pc[0];  p2  = pc[1];  p3  = pc[2];  p4  = pc[3];
372:       if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) {
373:         pv    = ba + 4*diag_offset[row];
374:         pj    = bj + diag_offset[row] + 1;
375:         x1    = pv[0];  x2  = pv[1];  x3  = pv[2];  x4  = pv[3];
376:         pc[0] = m1 = p1*x1 + p3*x2;
377:         pc[1] = m2 = p2*x1 + p4*x2;
378:         pc[2] = m3 = p1*x3 + p3*x4;
379:         pc[3] = m4 = p2*x3 + p4*x4;
380:         nz    = bi[row+1] - diag_offset[row] - 1;
381:         pv   += 4;
382:         for (j=0; j<nz; j++) {
383:           x1    = pv[0];  x2  = pv[1];   x3 = pv[2];  x4  = pv[3];
384:           x     = rtmp + 4*pj[j];
385:           x[0] -= m1*x1 + m3*x2;
386:           x[1] -= m2*x1 + m4*x2;
387:           x[2] -= m1*x3 + m3*x4;
388:           x[3] -= m2*x3 + m4*x4;
389:           pv   += 4;
390:         }
391:         PetscLogFlops(16.0*nz+12.0);
392:       }
393:       row = *ajtmp++;
394:     }
395:     /* finished row so stick it into b->a */
396:     pv = ba + 4*bi[i];
397:     pj = bj + bi[i];
398:     nz = bi[i+1] - bi[i];
399:     for (j=0; j<nz; j++) {
400:       x     = rtmp+4*pj[j];
401:       pv[0] = x[0];  pv[1]  = x[1];  pv[2]  = x[2];  pv[3]  = x[3];
402:       /*
403:       printf(" col %d:",pj[j]);
404:       PetscInt j1;
405:       for (j1=0; j1<4; j1++) printf(" %g,",*(pv+j1));
406:       printf("\n");
407:       */
408:       pv += 4;
409:     }
410:     /* invert diagonal block */
411:     w = ba + 4*diag_offset[i];
412:     PetscKernel_A_gets_inverse_A_2(w,shift, allowzeropivot,&zeropivotdetected);
413:     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
414:   }

416:   PetscFree(rtmp);

418:   C->ops->solve          = MatSolve_SeqBAIJ_2_NaturalOrdering_inplace;
419:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering_inplace;
420:   C->assembled           = PETSC_TRUE;

422:   PetscLogFlops(1.333333333333*8*b->mbs); /* from inverting diagonal blocks */
423:   return(0);
424: }

426: /* ----------------------------------------------------------- */
427: /*
428:      Version for when blocks are 1 by 1.
429: */
430: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1(Mat B,Mat A,const MatFactorInfo *info)
431: {
432:   Mat             C     =B;
433:   Mat_SeqBAIJ     *a    =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data;
434:   IS              isrow = b->row,isicol = b->icol;
435:   PetscErrorCode  ierr;
436:   const PetscInt  *r,*ic,*ics;
437:   const PetscInt  n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bdiag=b->diag;
438:   PetscInt        i,j,k,nz,nzL,row,*pj;
439:   const PetscInt  *ajtmp,*bjtmp;
440:   MatScalar       *rtmp,*pc,multiplier,*pv;
441:   const MatScalar *aa=a->a,*v;
442:   PetscBool       row_identity,col_identity;
443:   FactorShiftCtx  sctx;
444:   const PetscInt  *ddiag;
445:   PetscReal       rs;
446:   MatScalar       d;

449:   /* MatPivotSetUp(): initialize shift context sctx */
450:   PetscMemzero(&sctx,sizeof(FactorShiftCtx));

452:   if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
453:     ddiag          = a->diag;
454:     sctx.shift_top = info->zeropivot;
455:     for (i=0; i<n; i++) {
456:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
457:       d  = (aa)[ddiag[i]];
458:       rs = -PetscAbsScalar(d) - PetscRealPart(d);
459:       v  = aa+ai[i];
460:       nz = ai[i+1] - ai[i];
461:       for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
462:       if (rs>sctx.shift_top) sctx.shift_top = rs;
463:     }
464:     sctx.shift_top *= 1.1;
465:     sctx.nshift_max = 5;
466:     sctx.shift_lo   = 0.;
467:     sctx.shift_hi   = 1.;
468:   }

470:   ISGetIndices(isrow,&r);
471:   ISGetIndices(isicol,&ic);
472:   PetscMalloc1(n+1,&rtmp);
473:   ics  = ic;

475:   do {
476:     sctx.newshift = PETSC_FALSE;
477:     for (i=0; i<n; i++) {
478:       /* zero rtmp */
479:       /* L part */
480:       nz    = bi[i+1] - bi[i];
481:       bjtmp = bj + bi[i];
482:       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

484:       /* U part */
485:       nz    = bdiag[i]-bdiag[i+1];
486:       bjtmp = bj + bdiag[i+1]+1;
487:       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

489:       /* load in initial (unfactored row) */
490:       nz    = ai[r[i]+1] - ai[r[i]];
491:       ajtmp = aj + ai[r[i]];
492:       v     = aa + ai[r[i]];
493:       for (j=0; j<nz; j++) rtmp[ics[ajtmp[j]]] = v[j];

495:       /* ZeropivotApply() */
496:       rtmp[i] += sctx.shift_amount;  /* shift the diagonal of the matrix */

498:       /* elimination */
499:       bjtmp = bj + bi[i];
500:       row   = *bjtmp++;
501:       nzL   = bi[i+1] - bi[i];
502:       for (k=0; k < nzL; k++) {
503:         pc = rtmp + row;
504:         if (*pc != (PetscScalar)0.0) {
505:           pv         = b->a + bdiag[row];
506:           multiplier = *pc * (*pv);
507:           *pc        = multiplier;

509:           pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
510:           pv = b->a + bdiag[row+1]+1;
511:           nz = bdiag[row]-bdiag[row+1]-1; /* num of entries in U(row,:) excluding diag */
512:           for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
513:           PetscLogFlops(2.0*nz);
514:         }
515:         row = *bjtmp++;
516:       }

518:       /* finished row so stick it into b->a */
519:       rs = 0.0;
520:       /* L part */
521:       pv = b->a + bi[i];
522:       pj = b->j + bi[i];
523:       nz = bi[i+1] - bi[i];
524:       for (j=0; j<nz; j++) {
525:         pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
526:       }

528:       /* U part */
529:       pv = b->a + bdiag[i+1]+1;
530:       pj = b->j + bdiag[i+1]+1;
531:       nz = bdiag[i] - bdiag[i+1]-1;
532:       for (j=0; j<nz; j++) {
533:         pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
534:       }

536:       sctx.rs = rs;
537:       sctx.pv = rtmp[i];
538:       MatPivotCheck(B,A,info,&sctx,i);
539:       if (sctx.newshift) break; /* break for-loop */
540:       rtmp[i] = sctx.pv; /* sctx.pv might be updated in the case of MAT_SHIFT_INBLOCKS */

542:       /* Mark diagonal and invert diagonal for simpler triangular solves */
543:       pv  = b->a + bdiag[i];
544:       *pv = (PetscScalar)1.0/rtmp[i];

546:     } /* endof for (i=0; i<n; i++) { */

548:     /* MatPivotRefine() */
549:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
550:       /*
551:        * if no shift in this attempt & shifting & started shifting & can refine,
552:        * then try lower shift
553:        */
554:       sctx.shift_hi       = sctx.shift_fraction;
555:       sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
556:       sctx.shift_amount   = sctx.shift_fraction * sctx.shift_top;
557:       sctx.newshift       = PETSC_TRUE;
558:       sctx.nshift++;
559:     }
560:   } while (sctx.newshift);

562:   PetscFree(rtmp);
563:   ISRestoreIndices(isicol,&ic);
564:   ISRestoreIndices(isrow,&r);

566:   ISIdentity(isrow,&row_identity);
567:   ISIdentity(isicol,&col_identity);
568:   if (row_identity && col_identity) {
569:     C->ops->solve          = MatSolve_SeqBAIJ_1_NaturalOrdering;
570:     C->ops->forwardsolve   = MatForwardSolve_SeqBAIJ_1_NaturalOrdering;
571:     C->ops->backwardsolve  = MatBackwardSolve_SeqBAIJ_1_NaturalOrdering;
572:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering;
573:   } else {
574:     C->ops->solve          = MatSolve_SeqBAIJ_1;
575:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1;
576:   }
577:   C->assembled = PETSC_TRUE;
578:   PetscLogFlops(C->cmap->n);

580:   /* MatShiftView(A,info,&sctx) */
581:   if (sctx.nshift) {
582:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
583:       PetscInfo4(A,"number of shift_pd tries %D, shift_amount %g, diagonal shifted up by %e fraction top_value %e\n",sctx.nshift,(double)sctx.shift_amount,(double)sctx.shift_fraction,(double)sctx.shift_top);
584:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
585:       PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
586:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
587:       PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %g\n",sctx.nshift,(double)info->shiftamount);
588:     }
589:   }
590:   return(0);
591: }

593: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1_inplace(Mat C,Mat A,const MatFactorInfo *info)
594: {
595:   Mat_SeqBAIJ    *a    = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data;
596:   IS             isrow = b->row,isicol = b->icol;
598:   const PetscInt *r,*ic;
599:   PetscInt       i,j,n = a->mbs,*bi = b->i,*bj = b->j;
600:   PetscInt       *ajtmpold,*ajtmp,nz,row,*ai = a->i,*aj = a->j;
601:   PetscInt       *diag_offset = b->diag,diag,*pj;
602:   MatScalar      *pv,*v,*rtmp,multiplier,*pc;
603:   MatScalar      *ba = b->a,*aa = a->a;
604:   PetscBool      row_identity, col_identity;

607:   ISGetIndices(isrow,&r);
608:   ISGetIndices(isicol,&ic);
609:   PetscMalloc1(n+1,&rtmp);

611:   for (i=0; i<n; i++) {
612:     nz    = bi[i+1] - bi[i];
613:     ajtmp = bj + bi[i];
614:     for  (j=0; j<nz; j++) rtmp[ajtmp[j]] = 0.0;

616:     /* load in initial (unfactored row) */
617:     nz       = ai[r[i]+1] - ai[r[i]];
618:     ajtmpold = aj + ai[r[i]];
619:     v        = aa + ai[r[i]];
620:     for (j=0; j<nz; j++) rtmp[ic[ajtmpold[j]]] =  v[j];

622:     row = *ajtmp++;
623:     while (row < i) {
624:       pc = rtmp + row;
625:       if (*pc != 0.0) {
626:         pv         = ba + diag_offset[row];
627:         pj         = bj + diag_offset[row] + 1;
628:         multiplier = *pc * *pv++;
629:         *pc        = multiplier;
630:         nz         = bi[row+1] - diag_offset[row] - 1;
631:         for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
632:         PetscLogFlops(1.0+2.0*nz);
633:       }
634:       row = *ajtmp++;
635:     }
636:     /* finished row so stick it into b->a */
637:     pv = ba + bi[i];
638:     pj = bj + bi[i];
639:     nz = bi[i+1] - bi[i];
640:     for (j=0; j<nz; j++) pv[j] = rtmp[pj[j]];
641:     diag = diag_offset[i] - bi[i];
642:     /* check pivot entry for current row */
643:     if (pv[diag] == 0.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot: row in original ordering %D in permuted ordering %D",r[i],i);
644:     pv[diag] = 1.0/pv[diag];
645:   }

647:   PetscFree(rtmp);
648:   ISRestoreIndices(isicol,&ic);
649:   ISRestoreIndices(isrow,&r);
650:   ISIdentity(isrow,&row_identity);
651:   ISIdentity(isicol,&col_identity);
652:   if (row_identity && col_identity) {
653:     C->ops->solve          = MatSolve_SeqBAIJ_1_NaturalOrdering_inplace;
654:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering_inplace;
655:   } else {
656:     C->ops->solve          = MatSolve_SeqBAIJ_1_inplace;
657:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_inplace;
658:   }
659:   C->assembled = PETSC_TRUE;
660:   PetscLogFlops(C->cmap->n);
661:   return(0);
662: }

664: static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type)
665: {
667:   *type = MATSOLVERPETSC;
668:   return(0);
669: }

671: PETSC_INTERN PetscErrorCode MatGetFactor_seqbaij_petsc(Mat A,MatFactorType ftype,Mat *B)
672: {
673:   PetscInt       n = A->rmap->n;

677: #if defined(PETSC_USE_COMPLEX)
678:   if (A->hermitian && (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported");
679: #endif
680:   MatCreate(PetscObjectComm((PetscObject)A),B);
681:   MatSetSizes(*B,n,n,n,n);
682:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU || ftype == MAT_FACTOR_ILUDT) {
683:     MatSetType(*B,MATSEQBAIJ);

685:     (*B)->ops->lufactorsymbolic  = MatLUFactorSymbolic_SeqBAIJ;
686:     (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqBAIJ;
687:     PetscStrallocpy(MATORDERINGND,(char**)&(*B)->preferredordering[MAT_FACTOR_LU]);
688:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ILU]);
689:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ILUDT]);
690:   } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
691:     MatSetType(*B,MATSEQSBAIJ);
692:     MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);

694:     (*B)->ops->iccfactorsymbolic      = MatICCFactorSymbolic_SeqBAIJ;
695:     (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqBAIJ;
696:     /*  Future optimization would be direct symbolic and numerical factorization for BAIJ to support orderings and Cholesky, instead of first converting to SBAIJ */
697:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]);
698:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC]);
699:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
700:   (*B)->factortype = ftype;
701:   (*B)->canuseordering = PETSC_TRUE;

703:   PetscFree((*B)->solvertype);
704:   PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);
705:   PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc);
706:   return(0);
707: }

709: /* ----------------------------------------------------------- */
710: PetscErrorCode MatLUFactor_SeqBAIJ(Mat A,IS row,IS col,const MatFactorInfo *info)
711: {
713:   Mat            C;

716:   MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&C);
717:   MatLUFactorSymbolic(C,A,row,col,info);
718:   MatLUFactorNumeric(C,A,info);

720:   A->ops->solve          = C->ops->solve;
721:   A->ops->solvetranspose = C->ops->solvetranspose;

723:   MatHeaderMerge(A,&C);
724:   PetscLogObjectParent((PetscObject)A,(PetscObject)((Mat_SeqBAIJ*)(A->data))->icol);
725:   return(0);
726: }

728: #include <../src/mat/impls/sbaij/seq/sbaij.h>
729: PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N(Mat C,Mat A,const MatFactorInfo *info)
730: {
732:   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data;
733:   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
734:   IS             ip=b->row;
735:   const PetscInt *rip;
736:   PetscInt       i,j,mbs=a->mbs,bs=A->rmap->bs,*bi=b->i,*bj=b->j,*bcol;
737:   PetscInt       *ai=a->i,*aj=a->j;
738:   PetscInt       k,jmin,jmax,*jl,*il,col,nexti,ili,nz;
739:   MatScalar      *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
740:   PetscReal      rs;
741:   FactorShiftCtx sctx;

744:   if (bs > 1) { /* convert A to a SBAIJ matrix and apply Cholesky factorization from it */
745:     if (!a->sbaijMat) {
746:       MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
747:     }
748:     (a->sbaijMat)->ops->choleskyfactornumeric(C,a->sbaijMat,info);
749:     MatDestroy(&a->sbaijMat);
750:     return(0);
751:   }

753:   /* MatPivotSetUp(): initialize shift context sctx */
754:   PetscMemzero(&sctx,sizeof(FactorShiftCtx));

756:   ISGetIndices(ip,&rip);
757:   PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&jl);

759:   sctx.shift_amount = 0.;
760:   sctx.nshift       = 0;
761:   do {
762:     sctx.newshift = PETSC_FALSE;
763:     for (i=0; i<mbs; i++) {
764:       rtmp[i] = 0.0; jl[i] = mbs; il[0] = 0;
765:     }

767:     for (k = 0; k<mbs; k++) {
768:       bval = ba + bi[k];
769:       /* initialize k-th row by the perm[k]-th row of A */
770:       jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
771:       for (j = jmin; j < jmax; j++) {
772:         col = rip[aj[j]];
773:         if (col >= k) { /* only take upper triangular entry */
774:           rtmp[col] = aa[j];
775:           *bval++   = 0.0; /* for in-place factorization */
776:         }
777:       }

779:       /* shift the diagonal of the matrix */
780:       if (sctx.nshift) rtmp[k] += sctx.shift_amount;

782:       /* modify k-th row by adding in those rows i with U(i,k)!=0 */
783:       dk = rtmp[k];
784:       i  = jl[k]; /* first row to be added to k_th row  */

786:       while (i < k) {
787:         nexti = jl[i]; /* next row to be added to k_th row */

789:         /* compute multiplier, update diag(k) and U(i,k) */
790:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
791:         uikdi   = -ba[ili]*ba[bi[i]]; /* diagonal(k) */
792:         dk     += uikdi*ba[ili];
793:         ba[ili] = uikdi; /* -U(i,k) */

795:         /* add multiple of row i to k-th row */
796:         jmin = ili + 1; jmax = bi[i+1];
797:         if (jmin < jmax) {
798:           for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
799:           /* update il and jl for row i */
800:           il[i] = jmin;
801:           j     = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
802:         }
803:         i = nexti;
804:       }

806:       /* shift the diagonals when zero pivot is detected */
807:       /* compute rs=sum of abs(off-diagonal) */
808:       rs   = 0.0;
809:       jmin = bi[k]+1;
810:       nz   = bi[k+1] - jmin;
811:       if (nz) {
812:         bcol = bj + jmin;
813:         while (nz--) {
814:           rs += PetscAbsScalar(rtmp[*bcol]);
815:           bcol++;
816:         }
817:       }

819:       sctx.rs = rs;
820:       sctx.pv = dk;
821:       MatPivotCheck(C,A,info,&sctx,k);
822:       if (sctx.newshift) break;
823:       dk = sctx.pv;

825:       /* copy data into U(k,:) */
826:       ba[bi[k]] = 1.0/dk; /* U(k,k) */
827:       jmin      = bi[k]+1; jmax = bi[k+1];
828:       if (jmin < jmax) {
829:         for (j=jmin; j<jmax; j++) {
830:           col = bj[j]; ba[j] = rtmp[col]; rtmp[col] = 0.0;
831:         }
832:         /* add the k-th row into il and jl */
833:         il[k] = jmin;
834:         i     = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
835:       }
836:     }
837:   } while (sctx.newshift);
838:   PetscFree3(rtmp,il,jl);

840:   ISRestoreIndices(ip,&rip);

842:   C->assembled    = PETSC_TRUE;
843:   C->preallocated = PETSC_TRUE;

845:   PetscLogFlops(C->rmap->N);
846:   if (sctx.nshift) {
847:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
848:       PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
849:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
850:       PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
851:     }
852:   }
853:   return(0);
854: }

856: PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info)
857: {
858:   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data;
859:   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
861:   PetscInt       i,j,am=a->mbs;
862:   PetscInt       *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
863:   PetscInt       k,jmin,*jl,*il,nexti,ili,*acol,*bcol,nz;
864:   MatScalar      *rtmp,*ba=b->a,*aa=a->a,dk,uikdi,*aval,*bval;
865:   PetscReal      rs;
866:   FactorShiftCtx sctx;

869:   /* MatPivotSetUp(): initialize shift context sctx */
870:   PetscMemzero(&sctx,sizeof(FactorShiftCtx));

872:   PetscMalloc3(am,&rtmp,am,&il,am,&jl);

874:   do {
875:     sctx.newshift = PETSC_FALSE;
876:     for (i=0; i<am; i++) {
877:       rtmp[i] = 0.0; jl[i] = am; il[0] = 0;
878:     }

880:     for (k = 0; k<am; k++) {
881:       /* initialize k-th row with elements nonzero in row perm(k) of A */
882:       nz   = ai[k+1] - ai[k];
883:       acol = aj + ai[k];
884:       aval = aa + ai[k];
885:       bval = ba + bi[k];
886:       while (nz--) {
887:         if (*acol < k) { /* skip lower triangular entries */
888:           acol++; aval++;
889:         } else {
890:           rtmp[*acol++] = *aval++;
891:           *bval++       = 0.0; /* for in-place factorization */
892:         }
893:       }

895:       /* shift the diagonal of the matrix */
896:       if (sctx.nshift) rtmp[k] += sctx.shift_amount;

898:       /* modify k-th row by adding in those rows i with U(i,k)!=0 */
899:       dk = rtmp[k];
900:       i  = jl[k]; /* first row to be added to k_th row  */

902:       while (i < k) {
903:         nexti = jl[i]; /* next row to be added to k_th row */
904:         /* compute multiplier, update D(k) and U(i,k) */
905:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
906:         uikdi   = -ba[ili]*ba[bi[i]];
907:         dk     += uikdi*ba[ili];
908:         ba[ili] = uikdi; /* -U(i,k) */

910:         /* add multiple of row i to k-th row ... */
911:         jmin = ili + 1;
912:         nz   = bi[i+1] - jmin;
913:         if (nz > 0) {
914:           bcol = bj + jmin;
915:           bval = ba + jmin;
916:           while (nz--) rtmp[*bcol++] += uikdi*(*bval++);
917:           /* update il and jl for i-th row */
918:           il[i] = jmin;
919:           j     = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
920:         }
921:         i = nexti;
922:       }

924:       /* shift the diagonals when zero pivot is detected */
925:       /* compute rs=sum of abs(off-diagonal) */
926:       rs   = 0.0;
927:       jmin = bi[k]+1;
928:       nz   = bi[k+1] - jmin;
929:       if (nz) {
930:         bcol = bj + jmin;
931:         while (nz--) {
932:           rs += PetscAbsScalar(rtmp[*bcol]);
933:           bcol++;
934:         }
935:       }

937:       sctx.rs = rs;
938:       sctx.pv = dk;
939:       MatPivotCheck(C,A,info,&sctx,k);
940:       if (sctx.newshift) break;    /* sctx.shift_amount is updated */
941:       dk = sctx.pv;

943:       /* copy data into U(k,:) */
944:       ba[bi[k]] = 1.0/dk;
945:       jmin      = bi[k]+1;
946:       nz        = bi[k+1] - jmin;
947:       if (nz) {
948:         bcol = bj + jmin;
949:         bval = ba + jmin;
950:         while (nz--) {
951:           *bval++       = rtmp[*bcol];
952:           rtmp[*bcol++] = 0.0;
953:         }
954:         /* add k-th row into il and jl */
955:         il[k] = jmin;
956:         i     = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
957:       }
958:     }
959:   } while (sctx.newshift);
960:   PetscFree3(rtmp,il,jl);

962:   C->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
963:   C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
964:   C->assembled           = PETSC_TRUE;
965:   C->preallocated        = PETSC_TRUE;

967:   PetscLogFlops(C->rmap->N);
968:   if (sctx.nshift) {
969:     if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
970:       PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
971:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
972:       PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
973:     }
974:   }
975:   return(0);
976: }

978: #include <petscbt.h>
979: #include <../src/mat/utils/freespace.h>
980: PetscErrorCode MatICCFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
981: {
982:   Mat_SeqBAIJ        *a = (Mat_SeqBAIJ*)A->data;
983:   Mat_SeqSBAIJ       *b;
984:   Mat                B;
985:   PetscErrorCode     ierr;
986:   PetscBool          perm_identity,missing;
987:   PetscInt           reallocs=0,i,*ai=a->i,*aj=a->j,am=a->mbs,bs=A->rmap->bs,*ui;
988:   const PetscInt     *rip;
989:   PetscInt           jmin,jmax,nzk,k,j,*jl,prow,*il,nextprow;
990:   PetscInt           nlnk,*lnk,*lnk_lvl=NULL,ncols,ncols_upper,*cols,*cols_lvl,*uj,**uj_ptr,**uj_lvl_ptr;
991:   PetscReal          fill          =info->fill,levels=info->levels;
992:   PetscFreeSpaceList free_space    =NULL,current_space=NULL;
993:   PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL;
994:   PetscBT            lnkbt;

997:   MatMissingDiagonal(A,&missing,&i);
998:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);

1000:   if (bs > 1) {
1001:     if (!a->sbaijMat) {
1002:       MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
1003:     }
1004:     (fact)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ;  /* undue the change made in MatGetFactor_seqbaij_petsc */

1006:     MatICCFactorSymbolic(fact,a->sbaijMat,perm,info);
1007:     return(0);
1008:   }

1010:   ISIdentity(perm,&perm_identity);
1011:   ISGetIndices(perm,&rip);

1013:   /* special case that simply copies fill pattern */
1014:   if (!levels && perm_identity) {
1015:     PetscMalloc1(am+1,&ui);
1016:     for (i=0; i<am; i++) ui[i] = ai[i+1] - a->diag[i]; /* ui: rowlengths - changes when !perm_identity */
1017:     B    = fact;
1018:     MatSeqSBAIJSetPreallocation(B,1,0,ui);

1020:     b  = (Mat_SeqSBAIJ*)B->data;
1021:     uj = b->j;
1022:     for (i=0; i<am; i++) {
1023:       aj = a->j + a->diag[i];
1024:       for (j=0; j<ui[i]; j++) *uj++ = *aj++;
1025:       b->ilen[i] = ui[i];
1026:     }
1027:     PetscFree(ui);

1029:     B->factortype = MAT_FACTOR_NONE;

1031:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1032:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1033:     B->factortype = MAT_FACTOR_ICC;

1035:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1036:     return(0);
1037:   }

1039:   /* initialization */
1040:   PetscMalloc1(am+1,&ui);
1041:   ui[0] = 0;
1042:   PetscMalloc1(2*am+1,&cols_lvl);

1044:   /* jl: linked list for storing indices of the pivot rows
1045:      il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
1046:   PetscMalloc4(am,&uj_ptr,am,&uj_lvl_ptr,am,&il,am,&jl);
1047:   for (i=0; i<am; i++) {
1048:     jl[i] = am; il[i] = 0;
1049:   }

1051:   /* create and initialize a linked list for storing column indices of the active row k */
1052:   nlnk = am + 1;
1053:   PetscIncompleteLLCreate(am,am,nlnk,lnk,lnk_lvl,lnkbt);

1055:   /* initial FreeSpace size is fill*(ai[am]+am)/2 */
1056:   PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am]/2,am/2)),&free_space);

1058:   current_space = free_space;

1060:   PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am]/2,am/2)),&free_space_lvl);
1061:   current_space_lvl = free_space_lvl;

1063:   for (k=0; k<am; k++) {  /* for each active row k */
1064:     /* initialize lnk by the column indices of row rip[k] of A */
1065:     nzk         = 0;
1066:     ncols       = ai[rip[k]+1] - ai[rip[k]];
1067:     ncols_upper = 0;
1068:     cols        = cols_lvl + am;
1069:     for (j=0; j<ncols; j++) {
1070:       i = rip[*(aj + ai[rip[k]] + j)];
1071:       if (i >= k) { /* only take upper triangular entry */
1072:         cols[ncols_upper]     = i;
1073:         cols_lvl[ncols_upper] = -1;  /* initialize level for nonzero entries */
1074:         ncols_upper++;
1075:       }
1076:     }
1077:     PetscIncompleteLLAdd(ncols_upper,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);
1078:     nzk += nlnk;

1080:     /* update lnk by computing fill-in for each pivot row to be merged in */
1081:     prow = jl[k]; /* 1st pivot row */

1083:     while (prow < k) {
1084:       nextprow = jl[prow];

1086:       /* merge prow into k-th row */
1087:       jmin  = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
1088:       jmax  = ui[prow+1];
1089:       ncols = jmax-jmin;
1090:       i     = jmin - ui[prow];
1091:       cols  = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
1092:       for (j=0; j<ncols; j++) cols_lvl[j] = *(uj_lvl_ptr[prow] + i + j);
1093:       PetscIncompleteLLAddSorted(ncols,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);
1094:       nzk += nlnk;

1096:       /* update il and jl for prow */
1097:       if (jmin < jmax) {
1098:         il[prow] = jmin;

1100:         j = *cols; jl[prow] = jl[j]; jl[j] = prow;
1101:       }
1102:       prow = nextprow;
1103:     }

1105:     /* if free space is not available, make more free space */
1106:     if (current_space->local_remaining<nzk) {
1107:       i    = am - k + 1; /* num of unfactored rows */
1108:       i    = PetscMin(PetscIntMultTruncate(i,nzk), PetscIntMultTruncate(i,i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1109:       PetscFreeSpaceGet(i,&current_space);
1110:       PetscFreeSpaceGet(i,&current_space_lvl);
1111:       reallocs++;
1112:     }

1114:     /* copy data into free_space and free_space_lvl, then initialize lnk */
1115:     PetscIncompleteLLClean(am,am,nzk,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);

1117:     /* add the k-th row into il and jl */
1118:     if (nzk-1 > 0) {
1119:       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
1120:       jl[k] = jl[i]; jl[i] = k;
1121:       il[k] = ui[k] + 1;
1122:     }
1123:     uj_ptr[k]     = current_space->array;
1124:     uj_lvl_ptr[k] = current_space_lvl->array;

1126:     current_space->array           += nzk;
1127:     current_space->local_used      += nzk;
1128:     current_space->local_remaining -= nzk;

1130:     current_space_lvl->array           += nzk;
1131:     current_space_lvl->local_used      += nzk;
1132:     current_space_lvl->local_remaining -= nzk;

1134:     ui[k+1] = ui[k] + nzk;
1135:   }

1137:   ISRestoreIndices(perm,&rip);
1138:   PetscFree4(uj_ptr,uj_lvl_ptr,il,jl);
1139:   PetscFree(cols_lvl);

1141:   /* copy free_space into uj and free free_space; set uj in new datastructure; */
1142:   PetscMalloc1(ui[am]+1,&uj);
1143:   PetscFreeSpaceContiguous(&free_space,uj);
1144:   PetscIncompleteLLDestroy(lnk,lnkbt);
1145:   PetscFreeSpaceDestroy(free_space_lvl);

1147:   /* put together the new matrix in MATSEQSBAIJ format */
1148:   B    = fact;
1149:   MatSeqSBAIJSetPreallocation(B,1,MAT_SKIP_ALLOCATION,NULL);

1151:   b                = (Mat_SeqSBAIJ*)B->data;
1152:   b->singlemalloc  = PETSC_FALSE;
1153:   b->free_a        = PETSC_TRUE;
1154:   b->free_ij       = PETSC_TRUE;

1156:   PetscMalloc1(ui[am]+1,&b->a);

1158:   b->j             = uj;
1159:   b->i             = ui;
1160:   b->diag          = NULL;
1161:   b->ilen          = NULL;
1162:   b->imax          = NULL;
1163:   b->row           = perm;
1164:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */

1166:   PetscObjectReference((PetscObject)perm);

1168:   b->icol = perm;

1170:   PetscObjectReference((PetscObject)perm);
1171:   PetscMalloc1(am+1,&b->solve_work);
1172:   PetscLogObjectMemory((PetscObject)B,(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));

1174:   b->maxnz = b->nz = ui[am];

1176:   B->info.factor_mallocs   = reallocs;
1177:   B->info.fill_ratio_given = fill;
1178:   if (ai[am] != 0.) {
1179:     /* nonzeros in lower triangular part of A (includign diagonals)= (ai[am]+am)/2 */
1180:     B->info.fill_ratio_needed = ((PetscReal)2*ui[am])/(ai[am]+am);
1181:   } else {
1182:     B->info.fill_ratio_needed = 0.0;
1183:   }
1184: #if defined(PETSC_USE_INFO)
1185:   if (ai[am] != 0) {
1186:     PetscReal af = B->info.fill_ratio_needed;
1187:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
1188:     PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
1189:     PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
1190:   } else {
1191:     PetscInfo(A,"Empty matrix\n");
1192:   }
1193: #endif
1194:   if (perm_identity) {
1195:     B->ops->solve                 = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1196:     B->ops->solvetranspose        = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1197:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1198:   } else {
1199:     (fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
1200:   }
1201:   return(0);
1202: }

1204: PetscErrorCode MatCholeskyFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
1205: {
1206:   Mat_SeqBAIJ        *a = (Mat_SeqBAIJ*)A->data;
1207:   Mat_SeqSBAIJ       *b;
1208:   Mat                B;
1209:   PetscErrorCode     ierr;
1210:   PetscBool          perm_identity,missing;
1211:   PetscReal          fill = info->fill;
1212:   const PetscInt     *rip;
1213:   PetscInt           i,mbs=a->mbs,bs=A->rmap->bs,*ai=a->i,*aj=a->j,reallocs=0,prow;
1214:   PetscInt           *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
1215:   PetscInt           nlnk,*lnk,ncols,ncols_upper,*cols,*uj,**ui_ptr,*uj_ptr;
1216:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
1217:   PetscBT            lnkbt;

1220:   if (bs > 1) { /* convert to seqsbaij */
1221:     if (!a->sbaijMat) {
1222:       MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
1223:     }
1224:     (fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */

1226:     MatCholeskyFactorSymbolic(fact,a->sbaijMat,perm,info);
1227:     return(0);
1228:   }

1230:   MatMissingDiagonal(A,&missing,&i);
1231:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);

1233:   /* check whether perm is the identity mapping */
1234:   ISIdentity(perm,&perm_identity);
1235:   if (!perm_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported");
1236:   ISGetIndices(perm,&rip);

1238:   /* initialization */
1239:   PetscMalloc1(mbs+1,&ui);
1240:   ui[0] = 0;

1242:   /* jl: linked list for storing indices of the pivot rows
1243:      il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
1244:   PetscMalloc4(mbs,&ui_ptr,mbs,&il,mbs,&jl,mbs,&cols);
1245:   for (i=0; i<mbs; i++) {
1246:     jl[i] = mbs; il[i] = 0;
1247:   }

1249:   /* create and initialize a linked list for storing column indices of the active row k */
1250:   nlnk = mbs + 1;
1251:   PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);

1253:   /* initial FreeSpace size is fill* (ai[mbs]+mbs)/2 */
1254:   PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[mbs]/2,mbs/2)),&free_space);

1256:   current_space = free_space;

1258:   for (k=0; k<mbs; k++) {  /* for each active row k */
1259:     /* initialize lnk by the column indices of row rip[k] of A */
1260:     nzk         = 0;
1261:     ncols       = ai[rip[k]+1] - ai[rip[k]];
1262:     ncols_upper = 0;
1263:     for (j=0; j<ncols; j++) {
1264:       i = rip[*(aj + ai[rip[k]] + j)];
1265:       if (i >= k) { /* only take upper triangular entry */
1266:         cols[ncols_upper] = i;
1267:         ncols_upper++;
1268:       }
1269:     }
1270:     PetscLLAdd(ncols_upper,cols,mbs,nlnk,lnk,lnkbt);
1271:     nzk += nlnk;

1273:     /* update lnk by computing fill-in for each pivot row to be merged in */
1274:     prow = jl[k]; /* 1st pivot row */

1276:     while (prow < k) {
1277:       nextprow = jl[prow];
1278:       /* merge prow into k-th row */
1279:       jmin   = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
1280:       jmax   = ui[prow+1];
1281:       ncols  = jmax-jmin;
1282:       uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
1283:       PetscLLAddSorted(ncols,uj_ptr,mbs,nlnk,lnk,lnkbt);
1284:       nzk   += nlnk;

1286:       /* update il and jl for prow */
1287:       if (jmin < jmax) {
1288:         il[prow] = jmin;
1289:         j        = *uj_ptr;
1290:         jl[prow] = jl[j];
1291:         jl[j]    = prow;
1292:       }
1293:       prow = nextprow;
1294:     }

1296:     /* if free space is not available, make more free space */
1297:     if (current_space->local_remaining<nzk) {
1298:       i    = mbs - k + 1; /* num of unfactored rows */
1299:       i    = PetscMin(PetscIntMultTruncate(i,nzk), PetscIntMultTruncate(i,i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1300:       PetscFreeSpaceGet(i,&current_space);
1301:       reallocs++;
1302:     }

1304:     /* copy data into free space, then initialize lnk */
1305:     PetscLLClean(mbs,mbs,nzk,lnk,current_space->array,lnkbt);

1307:     /* add the k-th row into il and jl */
1308:     if (nzk-1 > 0) {
1309:       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
1310:       jl[k] = jl[i]; jl[i] = k;
1311:       il[k] = ui[k] + 1;
1312:     }
1313:     ui_ptr[k]                       = current_space->array;
1314:     current_space->array           += nzk;
1315:     current_space->local_used      += nzk;
1316:     current_space->local_remaining -= nzk;

1318:     ui[k+1] = ui[k] + nzk;
1319:   }

1321:   ISRestoreIndices(perm,&rip);
1322:   PetscFree4(ui_ptr,il,jl,cols);

1324:   /* copy free_space into uj and free free_space; set uj in new datastructure; */
1325:   PetscMalloc1(ui[mbs]+1,&uj);
1326:   PetscFreeSpaceContiguous(&free_space,uj);
1327:   PetscLLDestroy(lnk,lnkbt);

1329:   /* put together the new matrix in MATSEQSBAIJ format */
1330:   B    = fact;
1331:   MatSeqSBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,NULL);

1333:   b               = (Mat_SeqSBAIJ*)B->data;
1334:   b->singlemalloc = PETSC_FALSE;
1335:   b->free_a       = PETSC_TRUE;
1336:   b->free_ij      = PETSC_TRUE;

1338:   PetscMalloc1(ui[mbs]+1,&b->a);

1340:   b->j             = uj;
1341:   b->i             = ui;
1342:   b->diag          = NULL;
1343:   b->ilen          = NULL;
1344:   b->imax          = NULL;
1345:   b->row           = perm;
1346:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */

1348:   PetscObjectReference((PetscObject)perm);
1349:   b->icol  = perm;
1350:   PetscObjectReference((PetscObject)perm);
1351:   PetscMalloc1(mbs+1,&b->solve_work);
1352:   PetscLogObjectMemory((PetscObject)B,(ui[mbs]-mbs)*(sizeof(PetscInt)+sizeof(MatScalar)));
1353:   b->maxnz = b->nz = ui[mbs];

1355:   B->info.factor_mallocs   = reallocs;
1356:   B->info.fill_ratio_given = fill;
1357:   if (ai[mbs] != 0.) {
1358:     /* nonzeros in lower triangular part of A = (ai[mbs]+mbs)/2 */
1359:     B->info.fill_ratio_needed = ((PetscReal)2*ui[mbs])/(ai[mbs]+mbs);
1360:   } else {
1361:     B->info.fill_ratio_needed = 0.0;
1362:   }
1363: #if defined(PETSC_USE_INFO)
1364:   if (ai[mbs] != 0.) {
1365:     PetscReal af = B->info.fill_ratio_needed;
1366:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
1367:     PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
1368:     PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
1369:   } else {
1370:     PetscInfo(A,"Empty matrix\n");
1371:   }
1372: #endif
1373:   if (perm_identity) {
1374:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1375:   } else {
1376:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
1377:   }
1378:   return(0);
1379: }

1381: PetscErrorCode MatSolve_SeqBAIJ_N_NaturalOrdering(Mat A,Vec bb,Vec xx)
1382: {
1383:   Mat_SeqBAIJ       *a=(Mat_SeqBAIJ*)A->data;
1384:   PetscErrorCode    ierr;
1385:   const PetscInt    *ai=a->i,*aj=a->j,*adiag=a->diag,*vi;
1386:   PetscInt          i,k,n=a->mbs;
1387:   PetscInt          nz,bs=A->rmap->bs,bs2=a->bs2;
1388:   const MatScalar   *aa=a->a,*v;
1389:   PetscScalar       *x,*s,*t,*ls;
1390:   const PetscScalar *b;

1393:   VecGetArrayRead(bb,&b);
1394:   VecGetArray(xx,&x);
1395:   t    = a->solve_work;

1397:   /* forward solve the lower triangular */
1398:   PetscArraycpy(t,b,bs); /* copy 1st block of b to t */

1400:   for (i=1; i<n; i++) {
1401:     v    = aa + bs2*ai[i];
1402:     vi   = aj + ai[i];
1403:     nz   = ai[i+1] - ai[i];
1404:     s    = t + bs*i;
1405:     PetscArraycpy(s,b+bs*i,bs); /* copy i_th block of b to t */
1406:     for (k=0;k<nz;k++) {
1407:       PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[k]);
1408:       v += bs2;
1409:     }
1410:   }

1412:   /* backward solve the upper triangular */
1413:   ls = a->solve_work + A->cmap->n;
1414:   for (i=n-1; i>=0; i--) {
1415:     v    = aa + bs2*(adiag[i+1]+1);
1416:     vi   = aj + adiag[i+1]+1;
1417:     nz   = adiag[i] - adiag[i+1]-1;
1418:     PetscArraycpy(ls,t+i*bs,bs);
1419:     for (k=0; k<nz; k++) {
1420:       PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[k]);
1421:       v += bs2;
1422:     }
1423:     PetscKernel_w_gets_A_times_v(bs,ls,aa+bs2*adiag[i],t+i*bs); /* *inv(diagonal[i]) */
1424:     PetscArraycpy(x+i*bs,t+i*bs,bs);
1425:   }

1427:   VecRestoreArrayRead(bb,&b);
1428:   VecRestoreArray(xx,&x);
1429:   PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);
1430:   return(0);
1431: }

1433: PetscErrorCode MatSolve_SeqBAIJ_N(Mat A,Vec bb,Vec xx)
1434: {
1435:   Mat_SeqBAIJ        *a   =(Mat_SeqBAIJ*)A->data;
1436:   IS                 iscol=a->col,isrow=a->row;
1437:   PetscErrorCode     ierr;
1438:   const PetscInt     *r,*c,*rout,*cout,*ai=a->i,*aj=a->j,*adiag=a->diag,*vi;
1439:   PetscInt           i,m,n=a->mbs;
1440:   PetscInt           nz,bs=A->rmap->bs,bs2=a->bs2;
1441:   const MatScalar    *aa=a->a,*v;
1442:   PetscScalar        *x,*s,*t,*ls;
1443:   const PetscScalar  *b;

1446:   VecGetArrayRead(bb,&b);
1447:   VecGetArray(xx,&x);
1448:   t    = a->solve_work;

1450:   ISGetIndices(isrow,&rout); r = rout;
1451:   ISGetIndices(iscol,&cout); c = cout;

1453:   /* forward solve the lower triangular */
1454:   PetscArraycpy(t,b+bs*r[0],bs);
1455:   for (i=1; i<n; i++) {
1456:     v    = aa + bs2*ai[i];
1457:     vi   = aj + ai[i];
1458:     nz   = ai[i+1] - ai[i];
1459:     s    = t + bs*i;
1460:     PetscArraycpy(s,b+bs*r[i],bs);
1461:     for (m=0; m<nz; m++) {
1462:       PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[m]);
1463:       v += bs2;
1464:     }
1465:   }

1467:   /* backward solve the upper triangular */
1468:   ls = a->solve_work + A->cmap->n;
1469:   for (i=n-1; i>=0; i--) {
1470:     v    = aa + bs2*(adiag[i+1]+1);
1471:     vi   = aj + adiag[i+1]+1;
1472:     nz   = adiag[i] - adiag[i+1] - 1;
1473:     PetscArraycpy(ls,t+i*bs,bs);
1474:     for (m=0; m<nz; m++) {
1475:       PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[m]);
1476:       v += bs2;
1477:     }
1478:     PetscKernel_w_gets_A_times_v(bs,ls,v,t+i*bs); /* *inv(diagonal[i]) */
1479:     PetscArraycpy(x + bs*c[i],t+i*bs,bs);
1480:   }
1481:   ISRestoreIndices(isrow,&rout);
1482:   ISRestoreIndices(iscol,&cout);
1483:   VecRestoreArrayRead(bb,&b);
1484:   VecRestoreArray(xx,&x);
1485:   PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);
1486:   return(0);
1487: }

1489: /*
1490:     For each block in an block array saves the largest absolute value in the block into another array
1491: */
1492: static PetscErrorCode MatBlockAbs_private(PetscInt nbs,PetscInt bs2,PetscScalar *blockarray,PetscReal *absarray)
1493: {
1495:   PetscInt       i,j;

1498:   PetscArrayzero(absarray,nbs+1);
1499:   for (i=0; i<nbs; i++) {
1500:     for (j=0; j<bs2; j++) {
1501:       if (absarray[i] < PetscAbsScalar(blockarray[i*nbs+j])) absarray[i] = PetscAbsScalar(blockarray[i*nbs+j]);
1502:     }
1503:   }
1504:   return(0);
1505: }

1507: /*
1508:      This needs to be renamed and called by the regular MatILUFactor_SeqBAIJ when drop tolerance is used
1509: */
1510: PetscErrorCode MatILUDTFactor_SeqBAIJ(Mat A,IS isrow,IS iscol,const MatFactorInfo *info,Mat *fact)
1511: {
1512:   Mat            B = *fact;
1513:   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data,*b;
1514:   IS             isicol;
1516:   const PetscInt *r,*ic;
1517:   PetscInt       i,mbs=a->mbs,bs=A->rmap->bs,bs2=a->bs2,*ai=a->i,*aj=a->j,*ajtmp,*adiag;
1518:   PetscInt       *bi,*bj,*bdiag;

1520:   PetscInt  row,nzi,nzi_bl,nzi_bu,*im,dtcount,nzi_al,nzi_au;
1521:   PetscInt  nlnk,*lnk;
1522:   PetscBT   lnkbt;
1523:   PetscBool row_identity,icol_identity;
1524:   MatScalar *aatmp,*pv,*batmp,*ba,*rtmp,*pc,*multiplier,*vtmp;
1525:   PetscInt  j,nz,*pj,*bjtmp,k,ncut,*jtmp;

1527:   PetscReal dt=info->dt;          /* shift=info->shiftamount; */
1528:   PetscInt  nnz_max;
1529:   PetscBool missing;
1530:   PetscReal *vtmp_abs;
1531:   MatScalar *v_work;
1532:   PetscInt  *v_pivots;
1533:   PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE;

1536:   /* ------- symbolic factorization, can be reused ---------*/
1537:   MatMissingDiagonal(A,&missing,&i);
1538:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);
1539:   adiag=a->diag;

1541:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);

1543:   /* bdiag is location of diagonal in factor */
1544:   PetscMalloc1(mbs+1,&bdiag);

1546:   /* allocate row pointers bi */
1547:   PetscMalloc1(2*mbs+2,&bi);

1549:   /* allocate bj and ba; max num of nonzero entries is (ai[n]+2*n*dtcount+2) */
1550:   dtcount = (PetscInt)info->dtcount;
1551:   if (dtcount > mbs-1) dtcount = mbs-1;
1552:   nnz_max = ai[mbs]+2*mbs*dtcount +2;
1553:   /* printf("MatILUDTFactor_SeqBAIJ, bs %d, ai[mbs] %d, nnz_max  %d, dtcount %d\n",bs,ai[mbs],nnz_max,dtcount); */
1554:   PetscMalloc1(nnz_max,&bj);
1555:   nnz_max = nnz_max*bs2;
1556:   PetscMalloc1(nnz_max,&ba);

1558:   /* put together the new matrix */
1559:   MatSeqBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,NULL);
1560:   PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);

1562:   b               = (Mat_SeqBAIJ*)(B)->data;
1563:   b->free_a       = PETSC_TRUE;
1564:   b->free_ij      = PETSC_TRUE;
1565:   b->singlemalloc = PETSC_FALSE;

1567:   b->a    = ba;
1568:   b->j    = bj;
1569:   b->i    = bi;
1570:   b->diag = bdiag;
1571:   b->ilen = NULL;
1572:   b->imax = NULL;
1573:   b->row  = isrow;
1574:   b->col  = iscol;

1576:   PetscObjectReference((PetscObject)isrow);
1577:   PetscObjectReference((PetscObject)iscol);

1579:   b->icol  = isicol;
1580:   PetscMalloc1(bs*(mbs+1),&b->solve_work);
1581:   PetscLogObjectMemory((PetscObject)B,nnz_max*(sizeof(PetscInt)+sizeof(MatScalar)));
1582:   b->maxnz = nnz_max/bs2;

1584:   (B)->factortype            = MAT_FACTOR_ILUDT;
1585:   (B)->info.factor_mallocs   = 0;
1586:   (B)->info.fill_ratio_given = ((PetscReal)nnz_max)/((PetscReal)(ai[mbs]*bs2));
1587:   /* ------- end of symbolic factorization ---------*/
1588:   ISGetIndices(isrow,&r);
1589:   ISGetIndices(isicol,&ic);

1591:   /* linked list for storing column indices of the active row */
1592:   nlnk = mbs + 1;
1593:   PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);

1595:   /* im: used by PetscLLAddSortedLU(); jtmp: working array for column indices of active row */
1596:   PetscMalloc2(mbs,&im,mbs,&jtmp);
1597:   /* rtmp, vtmp: working arrays for sparse and contiguous row entries of active row */
1598:   PetscMalloc2(mbs*bs2,&rtmp,mbs*bs2,&vtmp);
1599:   PetscMalloc1(mbs+1,&vtmp_abs);
1600:   PetscMalloc3(bs,&v_work,bs2,&multiplier,bs,&v_pivots);

1602:   allowzeropivot = PetscNot(A->erroriffailure);
1603:   bi[0]       = 0;
1604:   bdiag[0]    = (nnz_max/bs2)-1; /* location of diagonal in factor B */
1605:   bi[2*mbs+1] = bdiag[0]+1; /* endof bj and ba array */
1606:   for (i=0; i<mbs; i++) {
1607:     /* copy initial fill into linked list */
1608:     nzi = ai[r[i]+1] - ai[r[i]];
1609:     if (!nzi) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",r[i],i);
1610:     nzi_al = adiag[r[i]] - ai[r[i]];
1611:     nzi_au = ai[r[i]+1] - adiag[r[i]] -1;

1613:     /* load in initial unfactored row */
1614:     ajtmp = aj + ai[r[i]];
1615:     PetscLLAddPerm(nzi,ajtmp,ic,mbs,nlnk,lnk,lnkbt);
1616:     PetscArrayzero(rtmp,mbs*bs2);
1617:     aatmp = a->a + bs2*ai[r[i]];
1618:     for (j=0; j<nzi; j++) {
1619:       PetscArraycpy(rtmp+bs2*ic[ajtmp[j]],aatmp+bs2*j,bs2);
1620:     }

1622:     /* add pivot rows into linked list */
1623:     row = lnk[mbs];
1624:     while (row < i) {
1625:       nzi_bl = bi[row+1] - bi[row] + 1;
1626:       bjtmp  = bj + bdiag[row+1]+1; /* points to 1st column next to the diagonal in U */
1627:       PetscLLAddSortedLU(bjtmp,row,nlnk,lnk,lnkbt,i,nzi_bl,im);
1628:       nzi   += nlnk;
1629:       row    = lnk[row];
1630:     }

1632:     /* copy data from lnk into jtmp, then initialize lnk */
1633:     PetscLLClean(mbs,mbs,nzi,lnk,jtmp,lnkbt);

1635:     /* numerical factorization */
1636:     bjtmp = jtmp;
1637:     row   = *bjtmp++; /* 1st pivot row */

1639:     while  (row < i) {
1640:       pc = rtmp + bs2*row;
1641:       pv = ba + bs2*bdiag[row]; /* inv(diag) of the pivot row */
1642:       PetscKernel_A_gets_A_times_B(bs,pc,pv,multiplier); /* pc= multiplier = pc*inv(diag[row]) */
1643:       MatBlockAbs_private(1,bs2,pc,vtmp_abs);
1644:       if (vtmp_abs[0] > dt) { /* apply tolerance dropping rule */
1645:         pj = bj + bdiag[row+1] + 1;         /* point to 1st entry of U(row,:) */
1646:         pv = ba + bs2*(bdiag[row+1] + 1);
1647:         nz = bdiag[row] - bdiag[row+1] - 1;         /* num of entries in U(row,:), excluding diagonal */
1648:         for (j=0; j<nz; j++) {
1649:           PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j);
1650:         }
1651:         /* PetscLogFlops(bslog*(nz+1.0)-bs); */
1652:       }
1653:       row = *bjtmp++;
1654:     }

1656:     /* copy sparse rtmp into contiguous vtmp; separate L and U part */
1657:     nzi_bl = 0; j = 0;
1658:     while (jtmp[j] < i) { /* L-part. Note: jtmp is sorted */
1659:       PetscArraycpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2);
1660:       nzi_bl++; j++;
1661:     }
1662:     nzi_bu = nzi - nzi_bl -1;

1664:     while (j < nzi) { /* U-part */
1665:       PetscArraycpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2);
1666:       j++;
1667:     }

1669:     MatBlockAbs_private(nzi,bs2,vtmp,vtmp_abs);

1671:     bjtmp = bj + bi[i];
1672:     batmp = ba + bs2*bi[i];
1673:     /* apply level dropping rule to L part */
1674:     ncut = nzi_al + dtcount;
1675:     if (ncut < nzi_bl) {
1676:       PetscSortSplitReal(ncut,nzi_bl,vtmp_abs,jtmp);
1677:       PetscSortIntWithScalarArray(ncut,jtmp,vtmp);
1678:     } else {
1679:       ncut = nzi_bl;
1680:     }
1681:     for (j=0; j<ncut; j++) {
1682:       bjtmp[j] = jtmp[j];
1683:       PetscArraycpy(batmp+bs2*j,rtmp+bs2*bjtmp[j],bs2);
1684:     }
1685:     bi[i+1] = bi[i] + ncut;
1686:     nzi     = ncut + 1;

1688:     /* apply level dropping rule to U part */
1689:     ncut = nzi_au + dtcount;
1690:     if (ncut < nzi_bu) {
1691:       PetscSortSplitReal(ncut,nzi_bu,vtmp_abs+nzi_bl+1,jtmp+nzi_bl+1);
1692:       PetscSortIntWithScalarArray(ncut,jtmp+nzi_bl+1,vtmp+nzi_bl+1);
1693:     } else {
1694:       ncut = nzi_bu;
1695:     }
1696:     nzi += ncut;

1698:     /* mark bdiagonal */
1699:     bdiag[i+1]    = bdiag[i] - (ncut + 1);
1700:     bi[2*mbs - i] = bi[2*mbs - i +1] - (ncut + 1);

1702:     bjtmp  = bj + bdiag[i];
1703:     batmp  = ba + bs2*bdiag[i];
1704:     PetscArraycpy(batmp,rtmp+bs2*i,bs2);
1705:     *bjtmp = i;

1707:     bjtmp = bj + bdiag[i+1]+1;
1708:     batmp = ba + (bdiag[i+1]+1)*bs2;

1710:     for (k=0; k<ncut; k++) {
1711:       bjtmp[k] = jtmp[nzi_bl+1+k];
1712:       PetscArraycpy(batmp+bs2*k,rtmp+bs2*bjtmp[k],bs2);
1713:     }

1715:     im[i] = nzi; /* used by PetscLLAddSortedLU() */

1717:     /* invert diagonal block for simpler triangular solves - add shift??? */
1718:     batmp = ba + bs2*bdiag[i];

1720:     PetscKernel_A_gets_inverse_A(bs,batmp,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
1721:     if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1722:   } /* for (i=0; i<mbs; i++) */
1723:   PetscFree3(v_work,multiplier,v_pivots);

1725:   /* printf("end of L %d, beginning of U %d\n",bi[mbs],bdiag[mbs]); */
1726:   if (bi[mbs] >= bdiag[mbs]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"end of L array %d cannot >= the beginning of U array %d",bi[mbs],bdiag[mbs]);

1728:   ISRestoreIndices(isrow,&r);
1729:   ISRestoreIndices(isicol,&ic);

1731:   PetscLLDestroy(lnk,lnkbt);

1733:   PetscFree2(im,jtmp);
1734:   PetscFree2(rtmp,vtmp);

1736:   PetscLogFlops(bs2*B->cmap->n);
1737:   b->maxnz = b->nz = bi[mbs] + bdiag[0] - bdiag[mbs];

1739:   ISIdentity(isrow,&row_identity);
1740:   ISIdentity(isicol,&icol_identity);
1741:   if (row_identity && icol_identity) {
1742:     B->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
1743:   } else {
1744:     B->ops->solve = MatSolve_SeqBAIJ_N;
1745:   }

1747:   B->ops->solveadd          = NULL;
1748:   B->ops->solvetranspose    = NULL;
1749:   B->ops->solvetransposeadd = NULL;
1750:   B->ops->matsolve          = NULL;
1751:   B->assembled              = PETSC_TRUE;
1752:   B->preallocated           = PETSC_TRUE;
1753:   return(0);
1754: }