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;
 13:   const PetscInt *r,*ic;
 14:   PetscInt       i,j,k,nz,nzL,row,*pj;
 15:   const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,bs2=a->bs2;
 16:   const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag;
 17:   MatScalar      *rtmp,*pc,*mwork,*pv;
 18:   MatScalar      *aa=a->a,*v;
 19:   PetscInt       flg;
 20:   PetscReal      shift = info->shiftamount;
 21:   PetscBool      allowzeropivot,zeropivotdetected;

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

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

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

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

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

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

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

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

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

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

111:   PetscFree2(rtmp,mwork);
112:   ISRestoreIndices(isicol,&ic);
113:   ISRestoreIndices(isrow,&r);

115:   C->ops->solve          = MatSolve_SeqBAIJ_2;
116:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2;
117:   C->assembled           = PETSC_TRUE;

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

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

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

138:   /* generate work space needed by the factorization */
139:   PetscMalloc2(bs2*n,&rtmp,bs2,&mwork);
140:   PetscArrayzero(rtmp,bs2*n);

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

151:     /* U part */
152:     nz    = bdiag[i] - bdiag[i+1];
153:     bjtmp = bj + bdiag[i+1]+1;
154:     for  (j=0; j<nz; j++) {
155:       PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);
156:     }

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

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

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

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

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

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

228:   C->ops->solve          = MatSolve_SeqBAIJ_2_NaturalOrdering;
229:   C->ops->forwardsolve   = MatForwardSolve_SeqBAIJ_2_NaturalOrdering;
230:   C->ops->backwardsolve  = MatBackwardSolve_SeqBAIJ_2_NaturalOrdering;
231:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering;
232:   C->assembled           = PETSC_TRUE;

234:   PetscLogFlops(1.333333333333*2*2*2*n); /* from inverting diagonal blocks */
235:   return 0;
236: }

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

253:   allowzeropivot = PetscNot(A->erroriffailure);
254:   ISGetIndices(isrow,&r);
255:   ISGetIndices(isicol,&ic);
256:   PetscMalloc1(4*(n+1),&rtmp);

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

316:   PetscFree(rtmp);
317:   ISRestoreIndices(isicol,&ic);
318:   ISRestoreIndices(isrow,&r);

320:   C->ops->solve          = MatSolve_SeqBAIJ_2_inplace;
321:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_inplace;
322:   C->assembled           = PETSC_TRUE;

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

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

408:   PetscFree(rtmp);

410:   C->ops->solve          = MatSolve_SeqBAIJ_2_NaturalOrdering_inplace;
411:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering_inplace;
412:   C->assembled           = PETSC_TRUE;

414:   PetscLogFlops(1.333333333333*8*b->mbs); /* from inverting diagonal blocks */
415:   return 0;
416: }

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

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

442:   if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
443:     ddiag          = a->diag;
444:     sctx.shift_top = info->zeropivot;
445:     for (i=0; i<n; i++) {
446:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
447:       d  = (aa)[ddiag[i]];
448:       rs = -PetscAbsScalar(d) - PetscRealPart(d);
449:       v  = aa+ai[i];
450:       nz = ai[i+1] - ai[i];
451:       for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
452:       if (rs>sctx.shift_top) sctx.shift_top = rs;
453:     }
454:     sctx.shift_top *= 1.1;
455:     sctx.nshift_max = 5;
456:     sctx.shift_lo   = 0.;
457:     sctx.shift_hi   = 1.;
458:   }

460:   ISGetIndices(isrow,&r);
461:   ISGetIndices(isicol,&ic);
462:   PetscMalloc1(n+1,&rtmp);
463:   ics  = ic;

465:   do {
466:     sctx.newshift = PETSC_FALSE;
467:     for (i=0; i<n; i++) {
468:       /* zero rtmp */
469:       /* L part */
470:       nz    = bi[i+1] - bi[i];
471:       bjtmp = bj + bi[i];
472:       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

474:       /* U part */
475:       nz    = bdiag[i]-bdiag[i+1];
476:       bjtmp = bj + bdiag[i+1]+1;
477:       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

479:       /* load in initial (unfactored row) */
480:       nz    = ai[r[i]+1] - ai[r[i]];
481:       ajtmp = aj + ai[r[i]];
482:       v     = aa + ai[r[i]];
483:       for (j=0; j<nz; j++) rtmp[ics[ajtmp[j]]] = v[j];

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

488:       /* elimination */
489:       bjtmp = bj + bi[i];
490:       row   = *bjtmp++;
491:       nzL   = bi[i+1] - bi[i];
492:       for (k=0; k < nzL; k++) {
493:         pc = rtmp + row;
494:         if (*pc != (PetscScalar)0.0) {
495:           pv         = b->a + bdiag[row];
496:           multiplier = *pc * (*pv);
497:           *pc        = multiplier;

499:           pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
500:           pv = b->a + bdiag[row+1]+1;
501:           nz = bdiag[row]-bdiag[row+1]-1; /* num of entries in U(row,:) excluding diag */
502:           for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
503:           PetscLogFlops(2.0*nz);
504:         }
505:         row = *bjtmp++;
506:       }

508:       /* finished row so stick it into b->a */
509:       rs = 0.0;
510:       /* L part */
511:       pv = b->a + bi[i];
512:       pj = b->j + bi[i];
513:       nz = bi[i+1] - bi[i];
514:       for (j=0; j<nz; j++) {
515:         pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
516:       }

518:       /* U part */
519:       pv = b->a + bdiag[i+1]+1;
520:       pj = b->j + bdiag[i+1]+1;
521:       nz = bdiag[i] - bdiag[i+1]-1;
522:       for (j=0; j<nz; j++) {
523:         pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
524:       }

526:       sctx.rs = rs;
527:       sctx.pv = rtmp[i];
528:       MatPivotCheck(B,A,info,&sctx,i);
529:       if (sctx.newshift) break; /* break for-loop */
530:       rtmp[i] = sctx.pv; /* sctx.pv might be updated in the case of MAT_SHIFT_INBLOCKS */

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

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

538:     /* MatPivotRefine() */
539:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
540:       /*
541:        * if no shift in this attempt & shifting & started shifting & can refine,
542:        * then try lower shift
543:        */
544:       sctx.shift_hi       = sctx.shift_fraction;
545:       sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
546:       sctx.shift_amount   = sctx.shift_fraction * sctx.shift_top;
547:       sctx.newshift       = PETSC_TRUE;
548:       sctx.nshift++;
549:     }
550:   } while (sctx.newshift);

552:   PetscFree(rtmp);
553:   ISRestoreIndices(isicol,&ic);
554:   ISRestoreIndices(isrow,&r);

556:   ISIdentity(isrow,&row_identity);
557:   ISIdentity(isicol,&col_identity);
558:   if (row_identity && col_identity) {
559:     C->ops->solve          = MatSolve_SeqBAIJ_1_NaturalOrdering;
560:     C->ops->forwardsolve   = MatForwardSolve_SeqBAIJ_1_NaturalOrdering;
561:     C->ops->backwardsolve  = MatBackwardSolve_SeqBAIJ_1_NaturalOrdering;
562:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering;
563:   } else {
564:     C->ops->solve          = MatSolve_SeqBAIJ_1;
565:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1;
566:   }
567:   C->assembled = PETSC_TRUE;
568:   PetscLogFlops(C->cmap->n);

570:   /* MatShiftView(A,info,&sctx) */
571:   if (sctx.nshift) {
572:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
573:       PetscInfo(A,"number of shift_pd tries %" PetscInt_FMT ", 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);
574:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
575:       PetscInfo(A,"number of shift_nz tries %" PetscInt_FMT ", shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
576:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
577:       PetscInfo(A,"number of shift_inblocks applied %" PetscInt_FMT ", each shift_amount %g\n",sctx.nshift,(double)info->shiftamount);
578:     }
579:   }
580:   return 0;
581: }

583: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1_inplace(Mat C,Mat A,const MatFactorInfo *info)
584: {
585:   Mat_SeqBAIJ    *a    = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data;
586:   IS             isrow = b->row,isicol = b->icol;
587:   const PetscInt *r,*ic;
588:   PetscInt       i,j,n = a->mbs,*bi = b->i,*bj = b->j;
589:   PetscInt       *ajtmpold,*ajtmp,nz,row,*ai = a->i,*aj = a->j;
590:   PetscInt       *diag_offset = b->diag,diag,*pj;
591:   MatScalar      *pv,*v,*rtmp,multiplier,*pc;
592:   MatScalar      *ba = b->a,*aa = a->a;
593:   PetscBool      row_identity, col_identity;

595:   ISGetIndices(isrow,&r);
596:   ISGetIndices(isicol,&ic);
597:   PetscMalloc1(n+1,&rtmp);

599:   for (i=0; i<n; i++) {
600:     nz    = bi[i+1] - bi[i];
601:     ajtmp = bj + bi[i];
602:     for  (j=0; j<nz; j++) rtmp[ajtmp[j]] = 0.0;

604:     /* load in initial (unfactored row) */
605:     nz       = ai[r[i]+1] - ai[r[i]];
606:     ajtmpold = aj + ai[r[i]];
607:     v        = aa + ai[r[i]];
608:     for (j=0; j<nz; j++) rtmp[ic[ajtmpold[j]]] =  v[j];

610:     row = *ajtmp++;
611:     while (row < i) {
612:       pc = rtmp + row;
613:       if (*pc != 0.0) {
614:         pv         = ba + diag_offset[row];
615:         pj         = bj + diag_offset[row] + 1;
616:         multiplier = *pc * *pv++;
617:         *pc        = multiplier;
618:         nz         = bi[row+1] - diag_offset[row] - 1;
619:         for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
620:         PetscLogFlops(1.0+2.0*nz);
621:       }
622:       row = *ajtmp++;
623:     }
624:     /* finished row so stick it into b->a */
625:     pv = ba + bi[i];
626:     pj = bj + bi[i];
627:     nz = bi[i+1] - bi[i];
628:     for (j=0; j<nz; j++) pv[j] = rtmp[pj[j]];
629:     diag = diag_offset[i] - bi[i];
630:     /* check pivot entry for current row */
632:     pv[diag] = 1.0/pv[diag];
633:   }

635:   PetscFree(rtmp);
636:   ISRestoreIndices(isicol,&ic);
637:   ISRestoreIndices(isrow,&r);
638:   ISIdentity(isrow,&row_identity);
639:   ISIdentity(isicol,&col_identity);
640:   if (row_identity && col_identity) {
641:     C->ops->solve          = MatSolve_SeqBAIJ_1_NaturalOrdering_inplace;
642:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering_inplace;
643:   } else {
644:     C->ops->solve          = MatSolve_SeqBAIJ_1_inplace;
645:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_inplace;
646:   }
647:   C->assembled = PETSC_TRUE;
648:   PetscLogFlops(C->cmap->n);
649:   return 0;
650: }

652: static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type)
653: {
654:   *type = MATSOLVERPETSC;
655:   return 0;
656: }

658: PETSC_INTERN PetscErrorCode MatGetFactor_seqbaij_petsc(Mat A,MatFactorType ftype,Mat *B)
659: {
660:   PetscInt       n = A->rmap->n;

662: #if defined(PETSC_USE_COMPLEX)
664: #endif
665:   MatCreate(PetscObjectComm((PetscObject)A),B);
666:   MatSetSizes(*B,n,n,n,n);
667:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU || ftype == MAT_FACTOR_ILUDT) {
668:     MatSetType(*B,MATSEQBAIJ);

670:     (*B)->ops->lufactorsymbolic  = MatLUFactorSymbolic_SeqBAIJ;
671:     (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqBAIJ;
672:     PetscStrallocpy(MATORDERINGND,(char**)&(*B)->preferredordering[MAT_FACTOR_LU]);
673:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ILU]);
674:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ILUDT]);
675:   } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
676:     MatSetType(*B,MATSEQSBAIJ);
677:     MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);

679:     (*B)->ops->iccfactorsymbolic      = MatICCFactorSymbolic_SeqBAIJ;
680:     (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqBAIJ;
681:     /*  Future optimization would be direct symbolic and numerical factorization for BAIJ to support orderings and Cholesky, instead of first converting to SBAIJ */
682:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]);
683:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC]);
684:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
685:   (*B)->factortype = ftype;
686:   (*B)->canuseordering = PETSC_TRUE;

688:   PetscFree((*B)->solvertype);
689:   PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);
690:   PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc);
691:   return 0;
692: }

694: /* ----------------------------------------------------------- */
695: PetscErrorCode MatLUFactor_SeqBAIJ(Mat A,IS row,IS col,const MatFactorInfo *info)
696: {
697:   Mat            C;

699:   MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&C);
700:   MatLUFactorSymbolic(C,A,row,col,info);
701:   MatLUFactorNumeric(C,A,info);

703:   A->ops->solve          = C->ops->solve;
704:   A->ops->solvetranspose = C->ops->solvetranspose;

706:   MatHeaderMerge(A,&C);
707:   PetscLogObjectParent((PetscObject)A,(PetscObject)((Mat_SeqBAIJ*)(A->data))->icol);
708:   return 0;
709: }

711: #include <../src/mat/impls/sbaij/seq/sbaij.h>
712: PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N(Mat C,Mat A,const MatFactorInfo *info)
713: {
714:   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data;
715:   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
716:   IS             ip=b->row;
717:   const PetscInt *rip;
718:   PetscInt       i,j,mbs=a->mbs,bs=A->rmap->bs,*bi=b->i,*bj=b->j,*bcol;
719:   PetscInt       *ai=a->i,*aj=a->j;
720:   PetscInt       k,jmin,jmax,*jl,*il,col,nexti,ili,nz;
721:   MatScalar      *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
722:   PetscReal      rs;
723:   FactorShiftCtx sctx;

725:   if (bs > 1) { /* convert A to a SBAIJ matrix and apply Cholesky factorization from it */
726:     if (!a->sbaijMat) {
727:       MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
728:     }
729:     (a->sbaijMat)->ops->choleskyfactornumeric(C,a->sbaijMat,info);
730:     MatDestroy(&a->sbaijMat);
731:     return 0;
732:   }

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

737:   ISGetIndices(ip,&rip);
738:   PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&jl);

740:   sctx.shift_amount = 0.;
741:   sctx.nshift       = 0;
742:   do {
743:     sctx.newshift = PETSC_FALSE;
744:     for (i=0; i<mbs; i++) {
745:       rtmp[i] = 0.0; jl[i] = mbs; il[0] = 0;
746:     }

748:     for (k = 0; k<mbs; k++) {
749:       bval = ba + bi[k];
750:       /* initialize k-th row by the perm[k]-th row of A */
751:       jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
752:       for (j = jmin; j < jmax; j++) {
753:         col = rip[aj[j]];
754:         if (col >= k) { /* only take upper triangular entry */
755:           rtmp[col] = aa[j];
756:           *bval++   = 0.0; /* for in-place factorization */
757:         }
758:       }

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

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

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

770:         /* compute multiplier, update diag(k) and U(i,k) */
771:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
772:         uikdi   = -ba[ili]*ba[bi[i]]; /* diagonal(k) */
773:         dk     += uikdi*ba[ili];
774:         ba[ili] = uikdi; /* -U(i,k) */

776:         /* add multiple of row i to k-th row */
777:         jmin = ili + 1; jmax = bi[i+1];
778:         if (jmin < jmax) {
779:           for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
780:           /* update il and jl for row i */
781:           il[i] = jmin;
782:           j     = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
783:         }
784:         i = nexti;
785:       }

787:       /* shift the diagonals when zero pivot is detected */
788:       /* compute rs=sum of abs(off-diagonal) */
789:       rs   = 0.0;
790:       jmin = bi[k]+1;
791:       nz   = bi[k+1] - jmin;
792:       if (nz) {
793:         bcol = bj + jmin;
794:         while (nz--) {
795:           rs += PetscAbsScalar(rtmp[*bcol]);
796:           bcol++;
797:         }
798:       }

800:       sctx.rs = rs;
801:       sctx.pv = dk;
802:       MatPivotCheck(C,A,info,&sctx,k);
803:       if (sctx.newshift) break;
804:       dk = sctx.pv;

806:       /* copy data into U(k,:) */
807:       ba[bi[k]] = 1.0/dk; /* U(k,k) */
808:       jmin      = bi[k]+1; jmax = bi[k+1];
809:       if (jmin < jmax) {
810:         for (j=jmin; j<jmax; j++) {
811:           col = bj[j]; ba[j] = rtmp[col]; rtmp[col] = 0.0;
812:         }
813:         /* add the k-th row into il and jl */
814:         il[k] = jmin;
815:         i     = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
816:       }
817:     }
818:   } while (sctx.newshift);
819:   PetscFree3(rtmp,il,jl);

821:   ISRestoreIndices(ip,&rip);

823:   C->assembled    = PETSC_TRUE;
824:   C->preallocated = PETSC_TRUE;

826:   PetscLogFlops(C->rmap->N);
827:   if (sctx.nshift) {
828:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
829:       PetscInfo(A,"number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
830:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
831:       PetscInfo(A,"number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
832:     }
833:   }
834:   return 0;
835: }

837: PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info)
838: {
839:   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data;
840:   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
841:   PetscInt       i,j,am=a->mbs;
842:   PetscInt       *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
843:   PetscInt       k,jmin,*jl,*il,nexti,ili,*acol,*bcol,nz;
844:   MatScalar      *rtmp,*ba=b->a,*aa=a->a,dk,uikdi,*aval,*bval;
845:   PetscReal      rs;
846:   FactorShiftCtx sctx;

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

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

853:   do {
854:     sctx.newshift = PETSC_FALSE;
855:     for (i=0; i<am; i++) {
856:       rtmp[i] = 0.0; jl[i] = am; il[0] = 0;
857:     }

859:     for (k = 0; k<am; k++) {
860:       /* initialize k-th row with elements nonzero in row perm(k) of A */
861:       nz   = ai[k+1] - ai[k];
862:       acol = aj + ai[k];
863:       aval = aa + ai[k];
864:       bval = ba + bi[k];
865:       while (nz--) {
866:         if (*acol < k) { /* skip lower triangular entries */
867:           acol++; aval++;
868:         } else {
869:           rtmp[*acol++] = *aval++;
870:           *bval++       = 0.0; /* for in-place factorization */
871:         }
872:       }

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

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

881:       while (i < k) {
882:         nexti = jl[i]; /* next row to be added to k_th row */
883:         /* compute multiplier, update D(k) and U(i,k) */
884:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
885:         uikdi   = -ba[ili]*ba[bi[i]];
886:         dk     += uikdi*ba[ili];
887:         ba[ili] = uikdi; /* -U(i,k) */

889:         /* add multiple of row i to k-th row ... */
890:         jmin = ili + 1;
891:         nz   = bi[i+1] - jmin;
892:         if (nz > 0) {
893:           bcol = bj + jmin;
894:           bval = ba + jmin;
895:           while (nz--) rtmp[*bcol++] += uikdi*(*bval++);
896:           /* update il and jl for i-th row */
897:           il[i] = jmin;
898:           j     = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
899:         }
900:         i = nexti;
901:       }

903:       /* shift the diagonals when zero pivot is detected */
904:       /* compute rs=sum of abs(off-diagonal) */
905:       rs   = 0.0;
906:       jmin = bi[k]+1;
907:       nz   = bi[k+1] - jmin;
908:       if (nz) {
909:         bcol = bj + jmin;
910:         while (nz--) {
911:           rs += PetscAbsScalar(rtmp[*bcol]);
912:           bcol++;
913:         }
914:       }

916:       sctx.rs = rs;
917:       sctx.pv = dk;
918:       MatPivotCheck(C,A,info,&sctx,k);
919:       if (sctx.newshift) break;    /* sctx.shift_amount is updated */
920:       dk = sctx.pv;

922:       /* copy data into U(k,:) */
923:       ba[bi[k]] = 1.0/dk;
924:       jmin      = bi[k]+1;
925:       nz        = bi[k+1] - jmin;
926:       if (nz) {
927:         bcol = bj + jmin;
928:         bval = ba + jmin;
929:         while (nz--) {
930:           *bval++       = rtmp[*bcol];
931:           rtmp[*bcol++] = 0.0;
932:         }
933:         /* add k-th row into il and jl */
934:         il[k] = jmin;
935:         i     = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
936:       }
937:     }
938:   } while (sctx.newshift);
939:   PetscFree3(rtmp,il,jl);

941:   C->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
942:   C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
943:   C->assembled           = PETSC_TRUE;
944:   C->preallocated        = PETSC_TRUE;

946:   PetscLogFlops(C->rmap->N);
947:   if (sctx.nshift) {
948:     if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
949:       PetscInfo(A,"number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
950:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
951:       PetscInfo(A,"number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
952:     }
953:   }
954:   return 0;
955: }

957: #include <petscbt.h>
958: #include <../src/mat/utils/freespace.h>
959: PetscErrorCode MatICCFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
960: {
961:   Mat_SeqBAIJ        *a = (Mat_SeqBAIJ*)A->data;
962:   Mat_SeqSBAIJ       *b;
963:   Mat                B;
964:   PetscBool          perm_identity,missing;
965:   PetscInt           reallocs=0,i,*ai=a->i,*aj=a->j,am=a->mbs,bs=A->rmap->bs,*ui;
966:   const PetscInt     *rip;
967:   PetscInt           jmin,jmax,nzk,k,j,*jl,prow,*il,nextprow;
968:   PetscInt           nlnk,*lnk,*lnk_lvl=NULL,ncols,ncols_upper,*cols,*cols_lvl,*uj,**uj_ptr,**uj_lvl_ptr;
969:   PetscReal          fill          =info->fill,levels=info->levels;
970:   PetscFreeSpaceList free_space    =NULL,current_space=NULL;
971:   PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL;
972:   PetscBT            lnkbt;

974:   MatMissingDiagonal(A,&missing,&i);

977:   if (bs > 1) {
978:     if (!a->sbaijMat) {
979:       MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
980:     }
981:     (fact)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ;  /* undue the change made in MatGetFactor_seqbaij_petsc */

983:     MatICCFactorSymbolic(fact,a->sbaijMat,perm,info);
984:     return 0;
985:   }

987:   ISIdentity(perm,&perm_identity);
988:   ISGetIndices(perm,&rip);

990:   /* special case that simply copies fill pattern */
991:   if (!levels && perm_identity) {
992:     PetscMalloc1(am+1,&ui);
993:     for (i=0; i<am; i++) ui[i] = ai[i+1] - a->diag[i]; /* ui: rowlengths - changes when !perm_identity */
994:     B    = fact;
995:     MatSeqSBAIJSetPreallocation(B,1,0,ui);

997:     b  = (Mat_SeqSBAIJ*)B->data;
998:     uj = b->j;
999:     for (i=0; i<am; i++) {
1000:       aj = a->j + a->diag[i];
1001:       for (j=0; j<ui[i]; j++) *uj++ = *aj++;
1002:       b->ilen[i] = ui[i];
1003:     }
1004:     PetscFree(ui);

1006:     B->factortype = MAT_FACTOR_NONE;

1008:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1009:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1010:     B->factortype = MAT_FACTOR_ICC;

1012:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1013:     return 0;
1014:   }

1016:   /* initialization */
1017:   PetscMalloc1(am+1,&ui);
1018:   ui[0] = 0;
1019:   PetscMalloc1(2*am+1,&cols_lvl);

1021:   /* jl: linked list for storing indices of the pivot rows
1022:      il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
1023:   PetscMalloc4(am,&uj_ptr,am,&uj_lvl_ptr,am,&il,am,&jl);
1024:   for (i=0; i<am; i++) {
1025:     jl[i] = am; il[i] = 0;
1026:   }

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

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

1035:   current_space = free_space;

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

1040:   for (k=0; k<am; k++) {  /* for each active row k */
1041:     /* initialize lnk by the column indices of row rip[k] of A */
1042:     nzk         = 0;
1043:     ncols       = ai[rip[k]+1] - ai[rip[k]];
1044:     ncols_upper = 0;
1045:     cols        = cols_lvl + am;
1046:     for (j=0; j<ncols; j++) {
1047:       i = rip[*(aj + ai[rip[k]] + j)];
1048:       if (i >= k) { /* only take upper triangular entry */
1049:         cols[ncols_upper]     = i;
1050:         cols_lvl[ncols_upper] = -1;  /* initialize level for nonzero entries */
1051:         ncols_upper++;
1052:       }
1053:     }
1054:     PetscIncompleteLLAdd(ncols_upper,cols,levels,cols_lvl,am,&nlnk,lnk,lnk_lvl,lnkbt);
1055:     nzk += nlnk;

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

1060:     while (prow < k) {
1061:       nextprow = jl[prow];

1063:       /* merge prow into k-th row */
1064:       jmin  = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
1065:       jmax  = ui[prow+1];
1066:       ncols = jmax-jmin;
1067:       i     = jmin - ui[prow];
1068:       cols  = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
1069:       for (j=0; j<ncols; j++) cols_lvl[j] = *(uj_lvl_ptr[prow] + i + j);
1070:       PetscIncompleteLLAddSorted(ncols,cols,levels,cols_lvl,am,&nlnk,lnk,lnk_lvl,lnkbt);
1071:       nzk += nlnk;

1073:       /* update il and jl for prow */
1074:       if (jmin < jmax) {
1075:         il[prow] = jmin;

1077:         j = *cols; jl[prow] = jl[j]; jl[j] = prow;
1078:       }
1079:       prow = nextprow;
1080:     }

1082:     /* if free space is not available, make more free space */
1083:     if (current_space->local_remaining<nzk) {
1084:       i    = am - k + 1; /* num of unfactored rows */
1085:       i    = PetscMin(PetscIntMultTruncate(i,nzk), PetscIntMultTruncate(i,i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1086:       PetscFreeSpaceGet(i,&current_space);
1087:       PetscFreeSpaceGet(i,&current_space_lvl);
1088:       reallocs++;
1089:     }

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

1094:     /* add the k-th row into il and jl */
1095:     if (nzk-1 > 0) {
1096:       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
1097:       jl[k] = jl[i]; jl[i] = k;
1098:       il[k] = ui[k] + 1;
1099:     }
1100:     uj_ptr[k]     = current_space->array;
1101:     uj_lvl_ptr[k] = current_space_lvl->array;

1103:     current_space->array           += nzk;
1104:     current_space->local_used      += nzk;
1105:     current_space->local_remaining -= nzk;

1107:     current_space_lvl->array           += nzk;
1108:     current_space_lvl->local_used      += nzk;
1109:     current_space_lvl->local_remaining -= nzk;

1111:     ui[k+1] = ui[k] + nzk;
1112:   }

1114:   ISRestoreIndices(perm,&rip);
1115:   PetscFree4(uj_ptr,uj_lvl_ptr,il,jl);
1116:   PetscFree(cols_lvl);

1118:   /* copy free_space into uj and free free_space; set uj in new datastructure; */
1119:   PetscMalloc1(ui[am]+1,&uj);
1120:   PetscFreeSpaceContiguous(&free_space,uj);
1121:   PetscIncompleteLLDestroy(lnk,lnkbt);
1122:   PetscFreeSpaceDestroy(free_space_lvl);

1124:   /* put together the new matrix in MATSEQSBAIJ format */
1125:   B    = fact;
1126:   MatSeqSBAIJSetPreallocation(B,1,MAT_SKIP_ALLOCATION,NULL);

1128:   b                = (Mat_SeqSBAIJ*)B->data;
1129:   b->singlemalloc  = PETSC_FALSE;
1130:   b->free_a        = PETSC_TRUE;
1131:   b->free_ij       = PETSC_TRUE;

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

1135:   b->j             = uj;
1136:   b->i             = ui;
1137:   b->diag          = NULL;
1138:   b->ilen          = NULL;
1139:   b->imax          = NULL;
1140:   b->row           = perm;
1141:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */

1143:   PetscObjectReference((PetscObject)perm);

1145:   b->icol = perm;

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

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

1153:   B->info.factor_mallocs   = reallocs;
1154:   B->info.fill_ratio_given = fill;
1155:   if (ai[am] != 0.) {
1156:     /* nonzeros in lower triangular part of A (includign diagonals)= (ai[am]+am)/2 */
1157:     B->info.fill_ratio_needed = ((PetscReal)2*ui[am])/(ai[am]+am);
1158:   } else {
1159:     B->info.fill_ratio_needed = 0.0;
1160:   }
1161: #if defined(PETSC_USE_INFO)
1162:   if (ai[am] != 0) {
1163:     PetscReal af = B->info.fill_ratio_needed;
1164:     PetscInfo(A,"Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
1165:     PetscInfo(A,"Run with -pc_factor_fill %g or use \n",(double)af);
1166:     PetscInfo(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
1167:   } else {
1168:     PetscInfo(A,"Empty matrix\n");
1169:   }
1170: #endif
1171:   if (perm_identity) {
1172:     B->ops->solve                 = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1173:     B->ops->solvetranspose        = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1174:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1175:   } else {
1176:     (fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
1177:   }
1178:   return 0;
1179: }

1181: PetscErrorCode MatCholeskyFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
1182: {
1183:   Mat_SeqBAIJ        *a = (Mat_SeqBAIJ*)A->data;
1184:   Mat_SeqSBAIJ       *b;
1185:   Mat                B;
1186:   PetscBool          perm_identity,missing;
1187:   PetscReal          fill = info->fill;
1188:   const PetscInt     *rip;
1189:   PetscInt           i,mbs=a->mbs,bs=A->rmap->bs,*ai=a->i,*aj=a->j,reallocs=0,prow;
1190:   PetscInt           *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
1191:   PetscInt           nlnk,*lnk,ncols,ncols_upper,*cols,*uj,**ui_ptr,*uj_ptr;
1192:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
1193:   PetscBT            lnkbt;

1195:   if (bs > 1) { /* convert to seqsbaij */
1196:     if (!a->sbaijMat) {
1197:       MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
1198:     }
1199:     (fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */

1201:     MatCholeskyFactorSymbolic(fact,a->sbaijMat,perm,info);
1202:     return 0;
1203:   }

1205:   MatMissingDiagonal(A,&missing,&i);

1208:   /* check whether perm is the identity mapping */
1209:   ISIdentity(perm,&perm_identity);
1211:   ISGetIndices(perm,&rip);

1213:   /* initialization */
1214:   PetscMalloc1(mbs+1,&ui);
1215:   ui[0] = 0;

1217:   /* jl: linked list for storing indices of the pivot rows
1218:      il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
1219:   PetscMalloc4(mbs,&ui_ptr,mbs,&il,mbs,&jl,mbs,&cols);
1220:   for (i=0; i<mbs; i++) {
1221:     jl[i] = mbs; il[i] = 0;
1222:   }

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

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

1231:   current_space = free_space;

1233:   for (k=0; k<mbs; k++) {  /* for each active row k */
1234:     /* initialize lnk by the column indices of row rip[k] of A */
1235:     nzk         = 0;
1236:     ncols       = ai[rip[k]+1] - ai[rip[k]];
1237:     ncols_upper = 0;
1238:     for (j=0; j<ncols; j++) {
1239:       i = rip[*(aj + ai[rip[k]] + j)];
1240:       if (i >= k) { /* only take upper triangular entry */
1241:         cols[ncols_upper] = i;
1242:         ncols_upper++;
1243:       }
1244:     }
1245:     PetscLLAdd(ncols_upper,cols,mbs,&nlnk,lnk,lnkbt);
1246:     nzk += nlnk;

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

1251:     while (prow < k) {
1252:       nextprow = jl[prow];
1253:       /* merge prow into k-th row */
1254:       jmin   = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
1255:       jmax   = ui[prow+1];
1256:       ncols  = jmax-jmin;
1257:       uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
1258:       PetscLLAddSorted(ncols,uj_ptr,mbs,&nlnk,lnk,lnkbt);
1259:       nzk   += nlnk;

1261:       /* update il and jl for prow */
1262:       if (jmin < jmax) {
1263:         il[prow] = jmin;
1264:         j        = *uj_ptr;
1265:         jl[prow] = jl[j];
1266:         jl[j]    = prow;
1267:       }
1268:       prow = nextprow;
1269:     }

1271:     /* if free space is not available, make more free space */
1272:     if (current_space->local_remaining<nzk) {
1273:       i    = mbs - k + 1; /* num of unfactored rows */
1274:       i    = PetscMin(PetscIntMultTruncate(i,nzk), PetscIntMultTruncate(i,i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1275:       PetscFreeSpaceGet(i,&current_space);
1276:       reallocs++;
1277:     }

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

1282:     /* add the k-th row into il and jl */
1283:     if (nzk-1 > 0) {
1284:       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
1285:       jl[k] = jl[i]; jl[i] = k;
1286:       il[k] = ui[k] + 1;
1287:     }
1288:     ui_ptr[k]                       = current_space->array;
1289:     current_space->array           += nzk;
1290:     current_space->local_used      += nzk;
1291:     current_space->local_remaining -= nzk;

1293:     ui[k+1] = ui[k] + nzk;
1294:   }

1296:   ISRestoreIndices(perm,&rip);
1297:   PetscFree4(ui_ptr,il,jl,cols);

1299:   /* copy free_space into uj and free free_space; set uj in new datastructure; */
1300:   PetscMalloc1(ui[mbs]+1,&uj);
1301:   PetscFreeSpaceContiguous(&free_space,uj);
1302:   PetscLLDestroy(lnk,lnkbt);

1304:   /* put together the new matrix in MATSEQSBAIJ format */
1305:   B    = fact;
1306:   MatSeqSBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,NULL);

1308:   b               = (Mat_SeqSBAIJ*)B->data;
1309:   b->singlemalloc = PETSC_FALSE;
1310:   b->free_a       = PETSC_TRUE;
1311:   b->free_ij      = PETSC_TRUE;

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

1315:   b->j             = uj;
1316:   b->i             = ui;
1317:   b->diag          = NULL;
1318:   b->ilen          = NULL;
1319:   b->imax          = NULL;
1320:   b->row           = perm;
1321:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */

1323:   PetscObjectReference((PetscObject)perm);
1324:   b->icol  = perm;
1325:   PetscObjectReference((PetscObject)perm);
1326:   PetscMalloc1(mbs+1,&b->solve_work);
1327:   PetscLogObjectMemory((PetscObject)B,(ui[mbs]-mbs)*(sizeof(PetscInt)+sizeof(MatScalar)));
1328:   b->maxnz = b->nz = ui[mbs];

1330:   B->info.factor_mallocs   = reallocs;
1331:   B->info.fill_ratio_given = fill;
1332:   if (ai[mbs] != 0.) {
1333:     /* nonzeros in lower triangular part of A = (ai[mbs]+mbs)/2 */
1334:     B->info.fill_ratio_needed = ((PetscReal)2*ui[mbs])/(ai[mbs]+mbs);
1335:   } else {
1336:     B->info.fill_ratio_needed = 0.0;
1337:   }
1338: #if defined(PETSC_USE_INFO)
1339:   if (ai[mbs] != 0.) {
1340:     PetscReal af = B->info.fill_ratio_needed;
1341:     PetscInfo(A,"Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
1342:     PetscInfo(A,"Run with -pc_factor_fill %g or use \n",(double)af);
1343:     PetscInfo(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
1344:   } else {
1345:     PetscInfo(A,"Empty matrix\n");
1346:   }
1347: #endif
1348:   if (perm_identity) {
1349:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1350:   } else {
1351:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
1352:   }
1353:   return 0;
1354: }

1356: PetscErrorCode MatSolve_SeqBAIJ_N_NaturalOrdering(Mat A,Vec bb,Vec xx)
1357: {
1358:   Mat_SeqBAIJ       *a=(Mat_SeqBAIJ*)A->data;
1359:   const PetscInt    *ai=a->i,*aj=a->j,*adiag=a->diag,*vi;
1360:   PetscInt          i,k,n=a->mbs;
1361:   PetscInt          nz,bs=A->rmap->bs,bs2=a->bs2;
1362:   const MatScalar   *aa=a->a,*v;
1363:   PetscScalar       *x,*s,*t,*ls;
1364:   const PetscScalar *b;

1366:   VecGetArrayRead(bb,&b);
1367:   VecGetArray(xx,&x);
1368:   t    = a->solve_work;

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

1373:   for (i=1; i<n; i++) {
1374:     v    = aa + bs2*ai[i];
1375:     vi   = aj + ai[i];
1376:     nz   = ai[i+1] - ai[i];
1377:     s    = t + bs*i;
1378:     PetscArraycpy(s,b+bs*i,bs); /* copy i_th block of b to t */
1379:     for (k=0;k<nz;k++) {
1380:       PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[k]);
1381:       v += bs2;
1382:     }
1383:   }

1385:   /* backward solve the upper triangular */
1386:   ls = a->solve_work + A->cmap->n;
1387:   for (i=n-1; i>=0; i--) {
1388:     v    = aa + bs2*(adiag[i+1]+1);
1389:     vi   = aj + adiag[i+1]+1;
1390:     nz   = adiag[i] - adiag[i+1]-1;
1391:     PetscArraycpy(ls,t+i*bs,bs);
1392:     for (k=0; k<nz; k++) {
1393:       PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[k]);
1394:       v += bs2;
1395:     }
1396:     PetscKernel_w_gets_A_times_v(bs,ls,aa+bs2*adiag[i],t+i*bs); /* *inv(diagonal[i]) */
1397:     PetscArraycpy(x+i*bs,t+i*bs,bs);
1398:   }

1400:   VecRestoreArrayRead(bb,&b);
1401:   VecRestoreArray(xx,&x);
1402:   PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);
1403:   return 0;
1404: }

1406: PetscErrorCode MatSolve_SeqBAIJ_N(Mat A,Vec bb,Vec xx)
1407: {
1408:   Mat_SeqBAIJ        *a   =(Mat_SeqBAIJ*)A->data;
1409:   IS                 iscol=a->col,isrow=a->row;
1410:   const PetscInt     *r,*c,*rout,*cout,*ai=a->i,*aj=a->j,*adiag=a->diag,*vi;
1411:   PetscInt           i,m,n=a->mbs;
1412:   PetscInt           nz,bs=A->rmap->bs,bs2=a->bs2;
1413:   const MatScalar    *aa=a->a,*v;
1414:   PetscScalar        *x,*s,*t,*ls;
1415:   const PetscScalar  *b;

1417:   VecGetArrayRead(bb,&b);
1418:   VecGetArray(xx,&x);
1419:   t    = a->solve_work;

1421:   ISGetIndices(isrow,&rout); r = rout;
1422:   ISGetIndices(iscol,&cout); c = cout;

1424:   /* forward solve the lower triangular */
1425:   PetscArraycpy(t,b+bs*r[0],bs);
1426:   for (i=1; i<n; i++) {
1427:     v    = aa + bs2*ai[i];
1428:     vi   = aj + ai[i];
1429:     nz   = ai[i+1] - ai[i];
1430:     s    = t + bs*i;
1431:     PetscArraycpy(s,b+bs*r[i],bs);
1432:     for (m=0; m<nz; m++) {
1433:       PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[m]);
1434:       v += bs2;
1435:     }
1436:   }

1438:   /* backward solve the upper triangular */
1439:   ls = a->solve_work + A->cmap->n;
1440:   for (i=n-1; i>=0; i--) {
1441:     v    = aa + bs2*(adiag[i+1]+1);
1442:     vi   = aj + adiag[i+1]+1;
1443:     nz   = adiag[i] - adiag[i+1] - 1;
1444:     PetscArraycpy(ls,t+i*bs,bs);
1445:     for (m=0; m<nz; m++) {
1446:       PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[m]);
1447:       v += bs2;
1448:     }
1449:     PetscKernel_w_gets_A_times_v(bs,ls,v,t+i*bs); /* *inv(diagonal[i]) */
1450:     PetscArraycpy(x + bs*c[i],t+i*bs,bs);
1451:   }
1452:   ISRestoreIndices(isrow,&rout);
1453:   ISRestoreIndices(iscol,&cout);
1454:   VecRestoreArrayRead(bb,&b);
1455:   VecRestoreArray(xx,&x);
1456:   PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);
1457:   return 0;
1458: }

1460: /*
1461:     For each block in an block array saves the largest absolute value in the block into another array
1462: */
1463: static PetscErrorCode MatBlockAbs_private(PetscInt nbs,PetscInt bs2,PetscScalar *blockarray,PetscReal *absarray)
1464: {
1465:   PetscInt       i,j;

1467:   PetscArrayzero(absarray,nbs+1);
1468:   for (i=0; i<nbs; i++) {
1469:     for (j=0; j<bs2; j++) {
1470:       if (absarray[i] < PetscAbsScalar(blockarray[i*nbs+j])) absarray[i] = PetscAbsScalar(blockarray[i*nbs+j]);
1471:     }
1472:   }
1473:   return 0;
1474: }

1476: /*
1477:      This needs to be renamed and called by the regular MatILUFactor_SeqBAIJ when drop tolerance is used
1478: */
1479: PetscErrorCode MatILUDTFactor_SeqBAIJ(Mat A,IS isrow,IS iscol,const MatFactorInfo *info,Mat *fact)
1480: {
1481:   Mat            B = *fact;
1482:   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data,*b;
1483:   IS             isicol;
1484:   const PetscInt *r,*ic;
1485:   PetscInt       i,mbs=a->mbs,bs=A->rmap->bs,bs2=a->bs2,*ai=a->i,*aj=a->j,*ajtmp,*adiag;
1486:   PetscInt       *bi,*bj,*bdiag;

1488:   PetscInt  row,nzi,nzi_bl,nzi_bu,*im,dtcount,nzi_al,nzi_au;
1489:   PetscInt  nlnk,*lnk;
1490:   PetscBT   lnkbt;
1491:   PetscBool row_identity,icol_identity;
1492:   MatScalar *aatmp,*pv,*batmp,*ba,*rtmp,*pc,*multiplier,*vtmp;
1493:   PetscInt  j,nz,*pj,*bjtmp,k,ncut,*jtmp;

1495:   PetscReal dt=info->dt;          /* shift=info->shiftamount; */
1496:   PetscInt  nnz_max;
1497:   PetscBool missing;
1498:   PetscReal *vtmp_abs;
1499:   MatScalar *v_work;
1500:   PetscInt  *v_pivots;
1501:   PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE;

1503:   /* ------- symbolic factorization, can be reused ---------*/
1504:   MatMissingDiagonal(A,&missing,&i);
1506:   adiag=a->diag;

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

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

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

1516:   /* allocate bj and ba; max num of nonzero entries is (ai[n]+2*n*dtcount+2) */
1517:   dtcount = (PetscInt)info->dtcount;
1518:   if (dtcount > mbs-1) dtcount = mbs-1;
1519:   nnz_max = ai[mbs]+2*mbs*dtcount +2;
1520:   /* printf("MatILUDTFactor_SeqBAIJ, bs %d, ai[mbs] %d, nnz_max  %d, dtcount %d\n",bs,ai[mbs],nnz_max,dtcount); */
1521:   PetscMalloc1(nnz_max,&bj);
1522:   nnz_max = nnz_max*bs2;
1523:   PetscMalloc1(nnz_max,&ba);

1525:   /* put together the new matrix */
1526:   MatSeqBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,NULL);
1527:   PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);

1529:   b               = (Mat_SeqBAIJ*)(B)->data;
1530:   b->free_a       = PETSC_TRUE;
1531:   b->free_ij      = PETSC_TRUE;
1532:   b->singlemalloc = PETSC_FALSE;

1534:   b->a    = ba;
1535:   b->j    = bj;
1536:   b->i    = bi;
1537:   b->diag = bdiag;
1538:   b->ilen = NULL;
1539:   b->imax = NULL;
1540:   b->row  = isrow;
1541:   b->col  = iscol;

1543:   PetscObjectReference((PetscObject)isrow);
1544:   PetscObjectReference((PetscObject)iscol);

1546:   b->icol  = isicol;
1547:   PetscMalloc1(bs*(mbs+1),&b->solve_work);
1548:   PetscLogObjectMemory((PetscObject)B,nnz_max*(sizeof(PetscInt)+sizeof(MatScalar)));
1549:   b->maxnz = nnz_max/bs2;

1551:   (B)->factortype            = MAT_FACTOR_ILUDT;
1552:   (B)->info.factor_mallocs   = 0;
1553:   (B)->info.fill_ratio_given = ((PetscReal)nnz_max)/((PetscReal)(ai[mbs]*bs2));
1554:   /* ------- end of symbolic factorization ---------*/
1555:   ISGetIndices(isrow,&r);
1556:   ISGetIndices(isicol,&ic);

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

1562:   /* im: used by PetscLLAddSortedLU(); jtmp: working array for column indices of active row */
1563:   PetscMalloc2(mbs,&im,mbs,&jtmp);
1564:   /* rtmp, vtmp: working arrays for sparse and contiguous row entries of active row */
1565:   PetscMalloc2(mbs*bs2,&rtmp,mbs*bs2,&vtmp);
1566:   PetscMalloc1(mbs+1,&vtmp_abs);
1567:   PetscMalloc3(bs,&v_work,bs2,&multiplier,bs,&v_pivots);

1569:   allowzeropivot = PetscNot(A->erroriffailure);
1570:   bi[0]       = 0;
1571:   bdiag[0]    = (nnz_max/bs2)-1; /* location of diagonal in factor B */
1572:   bi[2*mbs+1] = bdiag[0]+1; /* endof bj and ba array */
1573:   for (i=0; i<mbs; i++) {
1574:     /* copy initial fill into linked list */
1575:     nzi = ai[r[i]+1] - ai[r[i]];
1577:     nzi_al = adiag[r[i]] - ai[r[i]];
1578:     nzi_au = ai[r[i]+1] - adiag[r[i]] -1;

1580:     /* load in initial unfactored row */
1581:     ajtmp = aj + ai[r[i]];
1582:     PetscLLAddPerm(nzi,ajtmp,ic,mbs,&nlnk,lnk,lnkbt);
1583:     PetscArrayzero(rtmp,mbs*bs2);
1584:     aatmp = a->a + bs2*ai[r[i]];
1585:     for (j=0; j<nzi; j++) PetscArraycpy(rtmp+bs2*ic[ajtmp[j]],aatmp+bs2*j,bs2);

1587:     /* add pivot rows into linked list */
1588:     row = lnk[mbs];
1589:     while (row < i) {
1590:       nzi_bl = bi[row+1] - bi[row] + 1;
1591:       bjtmp  = bj + bdiag[row+1]+1; /* points to 1st column next to the diagonal in U */
1592:       PetscLLAddSortedLU(bjtmp,row,&nlnk,lnk,lnkbt,i,nzi_bl,im);
1593:       nzi   += nlnk;
1594:       row    = lnk[row];
1595:     }

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

1600:     /* numerical factorization */
1601:     bjtmp = jtmp;
1602:     row   = *bjtmp++; /* 1st pivot row */

1604:     while  (row < i) {
1605:       pc = rtmp + bs2*row;
1606:       pv = ba + bs2*bdiag[row]; /* inv(diag) of the pivot row */
1607:       PetscKernel_A_gets_A_times_B(bs,pc,pv,multiplier); /* pc= multiplier = pc*inv(diag[row]) */
1608:       MatBlockAbs_private(1,bs2,pc,vtmp_abs);
1609:       if (vtmp_abs[0] > dt) { /* apply tolerance dropping rule */
1610:         pj = bj + bdiag[row+1] + 1;         /* point to 1st entry of U(row,:) */
1611:         pv = ba + bs2*(bdiag[row+1] + 1);
1612:         nz = bdiag[row] - bdiag[row+1] - 1;         /* num of entries in U(row,:), excluding diagonal */
1613:         for (j=0; j<nz; j++) {
1614:           PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j);
1615:         }
1616:         /* PetscLogFlops(bslog*(nz+1.0)-bs); */
1617:       }
1618:       row = *bjtmp++;
1619:     }

1621:     /* copy sparse rtmp into contiguous vtmp; separate L and U part */
1622:     nzi_bl = 0; j = 0;
1623:     while (jtmp[j] < i) { /* L-part. Note: jtmp is sorted */
1624:       PetscArraycpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2);
1625:       nzi_bl++; j++;
1626:     }
1627:     nzi_bu = nzi - nzi_bl -1;

1629:     while (j < nzi) { /* U-part */
1630:       PetscArraycpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2);
1631:       j++;
1632:     }

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

1636:     bjtmp = bj + bi[i];
1637:     batmp = ba + bs2*bi[i];
1638:     /* apply level dropping rule to L part */
1639:     ncut = nzi_al + dtcount;
1640:     if (ncut < nzi_bl) {
1641:       PetscSortSplitReal(ncut,nzi_bl,vtmp_abs,jtmp);
1642:       PetscSortIntWithScalarArray(ncut,jtmp,vtmp);
1643:     } else {
1644:       ncut = nzi_bl;
1645:     }
1646:     for (j=0; j<ncut; j++) {
1647:       bjtmp[j] = jtmp[j];
1648:       PetscArraycpy(batmp+bs2*j,rtmp+bs2*bjtmp[j],bs2);
1649:     }
1650:     bi[i+1] = bi[i] + ncut;
1651:     nzi     = ncut + 1;

1653:     /* apply level dropping rule to U part */
1654:     ncut = nzi_au + dtcount;
1655:     if (ncut < nzi_bu) {
1656:       PetscSortSplitReal(ncut,nzi_bu,vtmp_abs+nzi_bl+1,jtmp+nzi_bl+1);
1657:       PetscSortIntWithScalarArray(ncut,jtmp+nzi_bl+1,vtmp+nzi_bl+1);
1658:     } else {
1659:       ncut = nzi_bu;
1660:     }
1661:     nzi += ncut;

1663:     /* mark bdiagonal */
1664:     bdiag[i+1]    = bdiag[i] - (ncut + 1);
1665:     bi[2*mbs - i] = bi[2*mbs - i +1] - (ncut + 1);

1667:     bjtmp  = bj + bdiag[i];
1668:     batmp  = ba + bs2*bdiag[i];
1669:     PetscArraycpy(batmp,rtmp+bs2*i,bs2);
1670:     *bjtmp = i;

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

1675:     for (k=0; k<ncut; k++) {
1676:       bjtmp[k] = jtmp[nzi_bl+1+k];
1677:       PetscArraycpy(batmp+bs2*k,rtmp+bs2*bjtmp[k],bs2);
1678:     }

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

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

1685:     PetscKernel_A_gets_inverse_A(bs,batmp,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
1686:     if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1687:   } /* for (i=0; i<mbs; i++) */
1688:   PetscFree3(v_work,multiplier,v_pivots);

1690:   /* printf("end of L %d, beginning of U %d\n",bi[mbs],bdiag[mbs]); */

1693:   ISRestoreIndices(isrow,&r);
1694:   ISRestoreIndices(isicol,&ic);

1696:   PetscLLDestroy(lnk,lnkbt);

1698:   PetscFree2(im,jtmp);
1699:   PetscFree2(rtmp,vtmp);

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

1704:   ISIdentity(isrow,&row_identity);
1705:   ISIdentity(isicol,&icol_identity);
1706:   if (row_identity && icol_identity) {
1707:     B->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
1708:   } else {
1709:     B->ops->solve = MatSolve_SeqBAIJ_N;
1710:   }

1712:   B->ops->solveadd          = NULL;
1713:   B->ops->solvetranspose    = NULL;
1714:   B->ops->solvetransposeadd = NULL;
1715:   B->ops->matsolve          = NULL;
1716:   B->assembled              = PETSC_TRUE;
1717:   B->preallocated           = PETSC_TRUE;
1718:   return 0;
1719: }