Actual source code: baijfact.c

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

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

 20:   PetscFunctionBegin;
 21:   PetscCall(ISGetIndices(isrow, &r));
 22:   PetscCall(ISGetIndices(isicol, &ic));

 24:   /* generate work space needed by the factorization */
 25:   PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
 26:   PetscCall(PetscArrayzero(rtmp, bs2 * n));

 28:   for (PetscInt i = 0; i < n; i++) {
 29:     /* zero rtmp */
 30:     /* L part */
 31:     PetscInt *pj;
 32:     PetscInt  nzL, nz = bi[i + 1] - bi[i];
 33:     bjtmp = bj + bi[i];
 34:     for (PetscInt j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));

 36:     /* U part */
 37:     nz    = bdiag[i] - bdiag[i + 1];
 38:     bjtmp = bj + bdiag[i + 1] + 1;
 39:     for (PetscInt j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));

 41:     /* load in initial (unfactored row) */
 42:     nz    = ai[r[i] + 1] - ai[r[i]];
 43:     ajtmp = aj + ai[r[i]];
 44:     v     = aa + bs2 * ai[r[i]];
 45:     for (PetscInt j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ic[ajtmp[j]], v + bs2 * j, bs2));

 47:     /* elimination */
 48:     bjtmp = bj + bi[i];
 49:     nzL   = bi[i + 1] - bi[i];
 50:     for (PetscInt k = 0; k < nzL; k++) {
 51:       PetscBool flg = PETSC_FALSE;
 52:       PetscInt  row = bjtmp[k];

 54:       pc = rtmp + bs2 * row;
 55:       for (PetscInt j = 0; j < bs2; j++) {
 56:         if (pc[j] != (PetscScalar)0.0) {
 57:           flg = PETSC_TRUE;
 58:           break;
 59:         }
 60:       }
 61:       if (flg) {
 62:         pv = b->a + bs2 * bdiag[row];
 63:         /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
 64:         PetscCall(PetscKernel_A_gets_A_times_B_2(pc, pv, mwork));

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

 80:     /* finished row so stick it into b->a */
 81:     /* L part */
 82:     pv = b->a + bs2 * bi[i];
 83:     pj = b->j + bi[i];
 84:     nz = bi[i + 1] - bi[i];
 85:     for (PetscInt j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));

 87:     /* Mark diagonal and invert diagonal for simpler triangular solves */
 88:     pv = b->a + bs2 * bdiag[i];
 89:     pj = b->j + bdiag[i];
 90:     PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
 91:     {
 92:       PetscBool zeropivotdetected;

 94:       PetscCall(PetscKernel_A_gets_inverse_A_2(pv, shift, allowzeropivot, &zeropivotdetected));
 95:       if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
 96:     }

 98:     /* U part */
 99:     pv = b->a + bs2 * (bdiag[i + 1] + 1);
100:     pj = b->j + bdiag[i + 1] + 1;
101:     nz = bdiag[i] - bdiag[i + 1] - 1;
102:     for (PetscInt j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
103:   }

105:   PetscCall(PetscFree2(rtmp, mwork));
106:   PetscCall(ISRestoreIndices(isicol, &ic));
107:   PetscCall(ISRestoreIndices(isrow, &r));

109:   C->ops->solve          = MatSolve_SeqBAIJ_2;
110:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2;
111:   C->assembled           = PETSC_TRUE;

113:   PetscCall(PetscLogFlops(1.333333333333 * 2 * 2 * 2 * n)); /* from inverting diagonal blocks */
114:   PetscFunctionReturn(PETSC_SUCCESS);
115: }

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

130:   PetscFunctionBegin;
131:   allowzeropivot = PetscNot(A->erroriffailure);

133:   /* generate work space needed by the factorization */
134:   PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
135:   PetscCall(PetscArrayzero(rtmp, bs2 * n));

137:   for (i = 0; i < n; i++) {
138:     /* zero rtmp */
139:     /* L part */
140:     nz    = bi[i + 1] - bi[i];
141:     bjtmp = bj + bi[i];
142:     for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));

144:     /* U part */
145:     nz    = bdiag[i] - bdiag[i + 1];
146:     bjtmp = bj + bdiag[i + 1] + 1;
147:     for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));

149:     /* load in initial (unfactored row) */
150:     nz    = ai[i + 1] - ai[i];
151:     ajtmp = aj + ai[i];
152:     v     = aa + bs2 * ai[i];
153:     for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ajtmp[j], v + bs2 * j, bs2));

155:     /* elimination */
156:     bjtmp = bj + bi[i];
157:     nzL   = bi[i + 1] - bi[i];
158:     for (k = 0; k < nzL; k++) {
159:       row = bjtmp[k];
160:       pc  = rtmp + bs2 * row;
161:       for (flg = 0, j = 0; j < bs2; j++) {
162:         if (pc[j] != (PetscScalar)0.0) {
163:           flg = 1;
164:           break;
165:         }
166:       }
167:       if (flg) {
168:         pv = b->a + bs2 * bdiag[row];
169:         /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
170:         PetscCall(PetscKernel_A_gets_A_times_B_2(pc, pv, mwork));

172:         pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
173:         pv = b->a + bs2 * (bdiag[row + 1] + 1);
174:         nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries in U(row,:) excluding diag */
175:         for (j = 0; j < nz; j++) {
176:           /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
177:           /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
178:           v = rtmp + 4 * pj[j];
179:           PetscCall(PetscKernel_A_gets_A_minus_B_times_C_2(v, pc, pv));
180:           pv += 4;
181:         }
182:         PetscCall(PetscLogFlops(16.0 * nz + 12)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
183:       }
184:     }

186:     /* finished row so stick it into b->a */
187:     /* L part */
188:     pv = b->a + bs2 * bi[i];
189:     pj = b->j + bi[i];
190:     nz = bi[i + 1] - bi[i];
191:     for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));

193:     /* Mark diagonal and invert diagonal for simpler triangular solves */
194:     pv = b->a + bs2 * bdiag[i];
195:     pj = b->j + bdiag[i];
196:     PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
197:     PetscCall(PetscKernel_A_gets_inverse_A_2(pv, shift, allowzeropivot, &zeropivotdetected));
198:     if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

200:     /* U part */
201:     /*
202:     pv = b->a + bs2*bi[2*n-i];
203:     pj = b->j + bi[2*n-i];
204:     nz = bi[2*n-i+1] - bi[2*n-i] - 1;
205:     */
206:     pv = b->a + bs2 * (bdiag[i + 1] + 1);
207:     pj = b->j + bdiag[i + 1] + 1;
208:     nz = bdiag[i] - bdiag[i + 1] - 1;
209:     for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
210:   }
211:   PetscCall(PetscFree2(rtmp, mwork));

213:   C->ops->solve          = MatSolve_SeqBAIJ_2_NaturalOrdering;
214:   C->ops->forwardsolve   = MatForwardSolve_SeqBAIJ_2_NaturalOrdering;
215:   C->ops->backwardsolve  = MatBackwardSolve_SeqBAIJ_2_NaturalOrdering;
216:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering;
217:   C->assembled           = PETSC_TRUE;

219:   PetscCall(PetscLogFlops(1.333333333333 * 2 * 2 * 2 * n)); /* from inverting diagonal blocks */
220:   PetscFunctionReturn(PETSC_SUCCESS);
221: }

223: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_inplace(Mat B, Mat A, const MatFactorInfo *info)
224: {
225:   Mat             C = B;
226:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
227:   IS              isrow = b->row, isicol = b->icol;
228:   const PetscInt *r, *ic;
229:   PetscInt        i, j, n = a->mbs, *bi = b->i, *bj = b->j;
230:   PetscInt       *ajtmpold, *ajtmp, nz, row;
231:   PetscInt       *diag_offset = b->diag, idx, *ai = a->i, *aj = a->j, *pj;
232:   MatScalar      *pv, *v, *rtmp, m1, m2, m3, m4, *pc, *w, *x, x1, x2, x3, x4;
233:   MatScalar       p1, p2, p3, p4;
234:   MatScalar      *ba = b->a, *aa = a->a;
235:   PetscReal       shift = info->shiftamount;
236:   PetscBool       allowzeropivot, zeropivotdetected;

238:   PetscFunctionBegin;
239:   allowzeropivot = PetscNot(A->erroriffailure);
240:   PetscCall(ISGetIndices(isrow, &r));
241:   PetscCall(ISGetIndices(isicol, &ic));
242:   PetscCall(PetscMalloc1(4 * (n + 1), &rtmp));

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

318:   PetscCall(PetscFree(rtmp));
319:   PetscCall(ISRestoreIndices(isicol, &ic));
320:   PetscCall(ISRestoreIndices(isrow, &r));

322:   C->ops->solve          = MatSolve_SeqBAIJ_2_inplace;
323:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_inplace;
324:   C->assembled           = PETSC_TRUE;

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

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

426:   PetscCall(PetscFree(rtmp));

428:   C->ops->solve          = MatSolve_SeqBAIJ_2_NaturalOrdering_inplace;
429:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering_inplace;
430:   C->assembled           = PETSC_TRUE;

432:   PetscCall(PetscLogFlops(1.333333333333 * 8 * b->mbs)); /* from inverting diagonal blocks */
433:   PetscFunctionReturn(PETSC_SUCCESS);
434: }

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

456:   PetscFunctionBegin;
457:   /* MatPivotSetUp(): initialize shift context sctx */
458:   PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));

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

478:   PetscCall(ISGetIndices(isrow, &r));
479:   PetscCall(ISGetIndices(isicol, &ic));
480:   PetscCall(PetscMalloc1(n + 1, &rtmp));
481:   ics = ic;

483:   do {
484:     sctx.newshift = PETSC_FALSE;
485:     for (i = 0; i < n; i++) {
486:       /* zero rtmp */
487:       /* L part */
488:       nz    = bi[i + 1] - bi[i];
489:       bjtmp = bj + bi[i];
490:       for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;

492:       /* U part */
493:       nz    = bdiag[i] - bdiag[i + 1];
494:       bjtmp = bj + bdiag[i + 1] + 1;
495:       for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;

497:       /* load in initial (unfactored row) */
498:       nz    = ai[r[i] + 1] - ai[r[i]];
499:       ajtmp = aj + ai[r[i]];
500:       v     = aa + ai[r[i]];
501:       for (j = 0; j < nz; j++) rtmp[ics[ajtmp[j]]] = v[j];

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

506:       /* elimination */
507:       bjtmp = bj + bi[i];
508:       row   = *bjtmp++;
509:       nzL   = bi[i + 1] - bi[i];
510:       for (k = 0; k < nzL; k++) {
511:         pc = rtmp + row;
512:         if (*pc != (PetscScalar)0.0) {
513:           pv         = b->a + bdiag[row];
514:           multiplier = *pc * (*pv);
515:           *pc        = multiplier;

517:           pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
518:           pv = b->a + bdiag[row + 1] + 1;
519:           nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries in U(row,:) excluding diag */
520:           for (j = 0; j < nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
521:           PetscCall(PetscLogFlops(2.0 * nz));
522:         }
523:         row = *bjtmp++;
524:       }

526:       /* finished row so stick it into b->a */
527:       rs = 0.0;
528:       /* L part */
529:       pv = b->a + bi[i];
530:       pj = b->j + bi[i];
531:       nz = bi[i + 1] - bi[i];
532:       for (j = 0; j < nz; j++) {
533:         pv[j] = rtmp[pj[j]];
534:         rs += PetscAbsScalar(pv[j]);
535:       }

537:       /* U part */
538:       pv = b->a + bdiag[i + 1] + 1;
539:       pj = b->j + bdiag[i + 1] + 1;
540:       nz = bdiag[i] - bdiag[i + 1] - 1;
541:       for (j = 0; j < nz; j++) {
542:         pv[j] = rtmp[pj[j]];
543:         rs += PetscAbsScalar(pv[j]);
544:       }

546:       sctx.rs = rs;
547:       sctx.pv = rtmp[i];
548:       PetscCall(MatPivotCheck(B, A, info, &sctx, i));
549:       if (sctx.newshift) break; /* break for-loop */
550:       rtmp[i] = sctx.pv;        /* sctx.pv might be updated in the case of MAT_SHIFT_INBLOCKS */

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

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

558:     /* MatPivotRefine() */
559:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction > 0 && sctx.nshift < sctx.nshift_max) {
560:       /*
561:        * if no shift in this attempt & shifting & started shifting & can refine,
562:        * then try lower shift
563:        */
564:       sctx.shift_hi       = sctx.shift_fraction;
565:       sctx.shift_fraction = (sctx.shift_hi + sctx.shift_lo) / 2.;
566:       sctx.shift_amount   = sctx.shift_fraction * sctx.shift_top;
567:       sctx.newshift       = PETSC_TRUE;
568:       sctx.nshift++;
569:     }
570:   } while (sctx.newshift);

572:   PetscCall(PetscFree(rtmp));
573:   PetscCall(ISRestoreIndices(isicol, &ic));
574:   PetscCall(ISRestoreIndices(isrow, &r));

576:   PetscCall(ISIdentity(isrow, &row_identity));
577:   PetscCall(ISIdentity(isicol, &col_identity));
578:   if (row_identity && col_identity) {
579:     C->ops->solve          = MatSolve_SeqBAIJ_1_NaturalOrdering;
580:     C->ops->forwardsolve   = MatForwardSolve_SeqBAIJ_1_NaturalOrdering;
581:     C->ops->backwardsolve  = MatBackwardSolve_SeqBAIJ_1_NaturalOrdering;
582:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering;
583:   } else {
584:     C->ops->solve          = MatSolve_SeqBAIJ_1;
585:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1;
586:   }
587:   C->assembled = PETSC_TRUE;
588:   PetscCall(PetscLogFlops(C->cmap->n));

590:   /* MatShiftView(A,info,&sctx) */
591:   if (sctx.nshift) {
592:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
593:       PetscCall(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));
594:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
595:       PetscCall(PetscInfo(A, "number of shift_nz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
596:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
597:       PetscCall(PetscInfo(A, "number of shift_inblocks applied %" PetscInt_FMT ", each shift_amount %g\n", sctx.nshift, (double)info->shiftamount));
598:     }
599:   }
600:   PetscFunctionReturn(PETSC_SUCCESS);
601: }

603: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1_inplace(Mat C, Mat A, const MatFactorInfo *info)
604: {
605:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
606:   IS              isrow = b->row, isicol = b->icol;
607:   const PetscInt *r, *ic;
608:   PetscInt        i, j, n = a->mbs, *bi = b->i, *bj = b->j;
609:   PetscInt       *ajtmpold, *ajtmp, nz, row, *ai = a->i, *aj = a->j;
610:   PetscInt       *diag_offset = b->diag, diag, *pj;
611:   MatScalar      *pv, *v, *rtmp, multiplier, *pc;
612:   MatScalar      *ba = b->a, *aa = a->a;
613:   PetscBool       row_identity, col_identity;

615:   PetscFunctionBegin;
616:   PetscCall(ISGetIndices(isrow, &r));
617:   PetscCall(ISGetIndices(isicol, &ic));
618:   PetscCall(PetscMalloc1(n + 1, &rtmp));

620:   for (i = 0; i < n; i++) {
621:     nz    = bi[i + 1] - bi[i];
622:     ajtmp = bj + bi[i];
623:     for (j = 0; j < nz; j++) rtmp[ajtmp[j]] = 0.0;

625:     /* load in initial (unfactored row) */
626:     nz       = ai[r[i] + 1] - ai[r[i]];
627:     ajtmpold = aj + ai[r[i]];
628:     v        = aa + ai[r[i]];
629:     for (j = 0; j < nz; j++) rtmp[ic[ajtmpold[j]]] = v[j];

631:     row = *ajtmp++;
632:     while (row < i) {
633:       pc = rtmp + row;
634:       if (*pc != 0.0) {
635:         pv         = ba + diag_offset[row];
636:         pj         = bj + diag_offset[row] + 1;
637:         multiplier = *pc * *pv++;
638:         *pc        = multiplier;
639:         nz         = bi[row + 1] - diag_offset[row] - 1;
640:         for (j = 0; j < nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
641:         PetscCall(PetscLogFlops(1.0 + 2.0 * nz));
642:       }
643:       row = *ajtmp++;
644:     }
645:     /* finished row so stick it into b->a */
646:     pv = ba + bi[i];
647:     pj = bj + bi[i];
648:     nz = bi[i + 1] - bi[i];
649:     for (j = 0; j < nz; j++) pv[j] = rtmp[pj[j]];
650:     diag = diag_offset[i] - bi[i];
651:     /* check pivot entry for current row */
652:     PetscCheck(pv[diag] != 0.0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot: row in original ordering %" PetscInt_FMT " in permuted ordering %" PetscInt_FMT, r[i], i);
653:     pv[diag] = 1.0 / pv[diag];
654:   }

656:   PetscCall(PetscFree(rtmp));
657:   PetscCall(ISRestoreIndices(isicol, &ic));
658:   PetscCall(ISRestoreIndices(isrow, &r));
659:   PetscCall(ISIdentity(isrow, &row_identity));
660:   PetscCall(ISIdentity(isicol, &col_identity));
661:   if (row_identity && col_identity) {
662:     C->ops->solve          = MatSolve_SeqBAIJ_1_NaturalOrdering_inplace;
663:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering_inplace;
664:   } else {
665:     C->ops->solve          = MatSolve_SeqBAIJ_1_inplace;
666:     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_inplace;
667:   }
668:   C->assembled = PETSC_TRUE;
669:   PetscCall(PetscLogFlops(C->cmap->n));
670:   PetscFunctionReturn(PETSC_SUCCESS);
671: }

673: static PetscErrorCode MatFactorGetSolverType_petsc(Mat A, MatSolverType *type)
674: {
675:   PetscFunctionBegin;
676:   *type = MATSOLVERPETSC;
677:   PetscFunctionReturn(PETSC_SUCCESS);
678: }

680: PETSC_INTERN PetscErrorCode MatGetFactor_seqbaij_petsc(Mat A, MatFactorType ftype, Mat *B)
681: {
682:   PetscInt n = A->rmap->n;

684:   PetscFunctionBegin;
685: #if defined(PETSC_USE_COMPLEX)
686:   PetscCheck(A->hermitian != PETSC_BOOL3_TRUE || !(ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC), PETSC_COMM_SELF, PETSC_ERR_SUP, "Hermitian Factor is not supported");
687: #endif
688:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
689:   PetscCall(MatSetSizes(*B, n, n, n, n));
690:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU || ftype == MAT_FACTOR_ILUDT) {
691:     PetscCall(MatSetType(*B, MATSEQBAIJ));

693:     (*B)->ops->lufactorsymbolic  = MatLUFactorSymbolic_SeqBAIJ;
694:     (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqBAIJ;
695:     PetscCall(PetscStrallocpy(MATORDERINGND, (char **)&(*B)->preferredordering[MAT_FACTOR_LU]));
696:     PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ILU]));
697:     PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ILUDT]));
698:   } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
699:     PetscCall(MatSetType(*B, MATSEQSBAIJ));
700:     PetscCall(MatSeqSBAIJSetPreallocation(*B, A->rmap->bs, MAT_SKIP_ALLOCATION, NULL));

702:     (*B)->ops->iccfactorsymbolic      = MatICCFactorSymbolic_SeqBAIJ;
703:     (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqBAIJ;
704:     /*  Future optimization would be direct symbolic and numerical factorization for BAIJ to support orderings and Cholesky, instead of first converting to SBAIJ */
705:     PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]));
706:     PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ICC]));
707:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");
708:   (*B)->factortype     = ftype;
709:   (*B)->canuseordering = PETSC_TRUE;

711:   PetscCall(PetscFree((*B)->solvertype));
712:   PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
713:   PetscCall(PetscObjectComposeFunction((PetscObject)*B, "MatFactorGetSolverType_C", MatFactorGetSolverType_petsc));
714:   PetscFunctionReturn(PETSC_SUCCESS);
715: }

717: PetscErrorCode MatLUFactor_SeqBAIJ(Mat A, IS row, IS col, const MatFactorInfo *info)
718: {
719:   Mat C;

721:   PetscFunctionBegin;
722:   PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &C));
723:   PetscCall(MatLUFactorSymbolic(C, A, row, col, info));
724:   PetscCall(MatLUFactorNumeric(C, A, info));

726:   A->ops->solve          = C->ops->solve;
727:   A->ops->solvetranspose = C->ops->solvetranspose;

729:   PetscCall(MatHeaderMerge(A, &C));
730:   PetscFunctionReturn(PETSC_SUCCESS);
731: }

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

747:   PetscFunctionBegin;
748:   if (bs > 1) { /* convert A to a SBAIJ matrix and apply Cholesky factorization from it */
749:     if (!a->sbaijMat) PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &a->sbaijMat));
750:     PetscCall(a->sbaijMat->ops->choleskyfactornumeric(C, a->sbaijMat, info));
751:     PetscCall(MatDestroy(&a->sbaijMat));
752:     PetscFunctionReturn(PETSC_SUCCESS);
753:   }

755:   /* MatPivotSetUp(): initialize shift context sctx */
756:   PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));

758:   PetscCall(ISGetIndices(ip, &rip));
759:   PetscCall(PetscMalloc3(mbs, &rtmp, mbs, &il, mbs, &jl));

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

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

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

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

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

794:         /* compute multiplier, update diag(k) and U(i,k) */
795:         ili   = il[i];                /* index of first nonzero element in U(i,k:bms-1) */
796:         uikdi = -ba[ili] * ba[bi[i]]; /* diagonal(k) */
797:         dk += uikdi * ba[ili];
798:         ba[ili] = uikdi; /* -U(i,k) */

800:         /* add multiple of row i to k-th row */
801:         jmin = ili + 1;
802:         jmax = bi[i + 1];
803:         if (jmin < jmax) {
804:           for (j = jmin; j < jmax; j++) rtmp[bj[j]] += uikdi * ba[j];
805:           /* update il and jl for row i */
806:           il[i] = jmin;
807:           j     = bj[jmin];
808:           jl[i] = jl[j];
809:           jl[j] = i;
810:         }
811:         i = nexti;
812:       }

814:       /* shift the diagonals when zero pivot is detected */
815:       /* compute rs=sum of abs(off-diagonal) */
816:       rs   = 0.0;
817:       jmin = bi[k] + 1;
818:       nz   = bi[k + 1] - jmin;
819:       if (nz) {
820:         bcol = bj + jmin;
821:         while (nz--) {
822:           rs += PetscAbsScalar(rtmp[*bcol]);
823:           bcol++;
824:         }
825:       }

827:       sctx.rs = rs;
828:       sctx.pv = dk;
829:       PetscCall(MatPivotCheck(C, A, info, &sctx, k));
830:       if (sctx.newshift) break;
831:       dk = sctx.pv;

833:       /* copy data into U(k,:) */
834:       ba[bi[k]] = 1.0 / dk; /* U(k,k) */
835:       jmin      = bi[k] + 1;
836:       jmax      = bi[k + 1];
837:       if (jmin < jmax) {
838:         for (j = jmin; j < jmax; j++) {
839:           col       = bj[j];
840:           ba[j]     = rtmp[col];
841:           rtmp[col] = 0.0;
842:         }
843:         /* add the k-th row into il and jl */
844:         il[k] = jmin;
845:         i     = bj[jmin];
846:         jl[k] = jl[i];
847:         jl[i] = k;
848:       }
849:     }
850:   } while (sctx.newshift);
851:   PetscCall(PetscFree3(rtmp, il, jl));

853:   PetscCall(ISRestoreIndices(ip, &rip));

855:   C->assembled    = PETSC_TRUE;
856:   C->preallocated = PETSC_TRUE;

858:   PetscCall(PetscLogFlops(C->rmap->N));
859:   if (sctx.nshift) {
860:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
861:       PetscCall(PetscInfo(A, "number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
862:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
863:       PetscCall(PetscInfo(A, "number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
864:     }
865:   }
866:   PetscFunctionReturn(PETSC_SUCCESS);
867: }

869: PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering(Mat C, Mat A, const MatFactorInfo *info)
870: {
871:   Mat_SeqBAIJ   *a = (Mat_SeqBAIJ *)A->data;
872:   Mat_SeqSBAIJ  *b = (Mat_SeqSBAIJ *)C->data;
873:   PetscInt       i, j, am = a->mbs;
874:   PetscInt      *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
875:   PetscInt       k, jmin, *jl, *il, nexti, ili, *acol, *bcol, nz;
876:   MatScalar     *rtmp, *ba = b->a, *aa = a->a, dk, uikdi, *aval, *bval;
877:   PetscReal      rs;
878:   FactorShiftCtx sctx;

880:   PetscFunctionBegin;
881:   /* MatPivotSetUp(): initialize shift context sctx */
882:   PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));

884:   PetscCall(PetscMalloc3(am, &rtmp, am, &il, am, &jl));

886:   do {
887:     sctx.newshift = PETSC_FALSE;
888:     for (i = 0; i < am; i++) {
889:       rtmp[i] = 0.0;
890:       jl[i]   = am;
891:       il[0]   = 0;
892:     }

894:     for (k = 0; k < am; k++) {
895:       /* initialize k-th row with elements nonzero in row perm(k) of A */
896:       nz   = ai[k + 1] - ai[k];
897:       acol = aj + ai[k];
898:       aval = aa + ai[k];
899:       bval = ba + bi[k];
900:       while (nz--) {
901:         if (*acol < k) { /* skip lower triangular entries */
902:           acol++;
903:           aval++;
904:         } else {
905:           rtmp[*acol++] = *aval++;
906:           *bval++       = 0.0; /* for in-place factorization */
907:         }
908:       }

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

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

917:       while (i < k) {
918:         nexti = jl[i]; /* next row to be added to k_th row */
919:         /* compute multiplier, update D(k) and U(i,k) */
920:         ili   = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
921:         uikdi = -ba[ili] * ba[bi[i]];
922:         dk += uikdi * ba[ili];
923:         ba[ili] = uikdi; /* -U(i,k) */

925:         /* add multiple of row i to k-th row ... */
926:         jmin = ili + 1;
927:         nz   = bi[i + 1] - jmin;
928:         if (nz > 0) {
929:           bcol = bj + jmin;
930:           bval = ba + jmin;
931:           while (nz--) rtmp[*bcol++] += uikdi * (*bval++);
932:           /* update il and jl for i-th row */
933:           il[i] = jmin;
934:           j     = bj[jmin];
935:           jl[i] = jl[j];
936:           jl[j] = i;
937:         }
938:         i = nexti;
939:       }

941:       /* shift the diagonals when zero pivot is detected */
942:       /* compute rs=sum of abs(off-diagonal) */
943:       rs   = 0.0;
944:       jmin = bi[k] + 1;
945:       nz   = bi[k + 1] - jmin;
946:       if (nz) {
947:         bcol = bj + jmin;
948:         while (nz--) {
949:           rs += PetscAbsScalar(rtmp[*bcol]);
950:           bcol++;
951:         }
952:       }

954:       sctx.rs = rs;
955:       sctx.pv = dk;
956:       PetscCall(MatPivotCheck(C, A, info, &sctx, k));
957:       if (sctx.newshift) break; /* sctx.shift_amount is updated */
958:       dk = sctx.pv;

960:       /* copy data into U(k,:) */
961:       ba[bi[k]] = 1.0 / dk;
962:       jmin      = bi[k] + 1;
963:       nz        = bi[k + 1] - jmin;
964:       if (nz) {
965:         bcol = bj + jmin;
966:         bval = ba + jmin;
967:         while (nz--) {
968:           *bval++       = rtmp[*bcol];
969:           rtmp[*bcol++] = 0.0;
970:         }
971:         /* add k-th row into il and jl */
972:         il[k] = jmin;
973:         i     = bj[jmin];
974:         jl[k] = jl[i];
975:         jl[i] = k;
976:       }
977:     }
978:   } while (sctx.newshift);
979:   PetscCall(PetscFree3(rtmp, il, jl));

981:   C->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
982:   C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
983:   C->assembled           = PETSC_TRUE;
984:   C->preallocated        = PETSC_TRUE;

986:   PetscCall(PetscLogFlops(C->rmap->N));
987:   if (sctx.nshift) {
988:     if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
989:       PetscCall(PetscInfo(A, "number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
990:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
991:       PetscCall(PetscInfo(A, "number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
992:     }
993:   }
994:   PetscFunctionReturn(PETSC_SUCCESS);
995: }

997: #include <petscbt.h>
998: #include <../src/mat/utils/freespace.h>
999: PetscErrorCode MatICCFactorSymbolic_SeqBAIJ(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
1000: {
1001:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ *)A->data;
1002:   Mat_SeqSBAIJ      *b;
1003:   Mat                B;
1004:   PetscBool          perm_identity, missing;
1005:   PetscInt           reallocs = 0, i, *ai = a->i, *aj = a->j, am = a->mbs, bs = A->rmap->bs, *ui;
1006:   const PetscInt    *rip;
1007:   PetscInt           jmin, jmax, nzk, k, j, *jl, prow, *il, nextprow;
1008:   PetscInt           nlnk, *lnk, *lnk_lvl = NULL, ncols, ncols_upper, *cols, *cols_lvl, *uj, **uj_ptr, **uj_lvl_ptr;
1009:   PetscReal          fill = info->fill, levels = info->levels;
1010:   PetscFreeSpaceList free_space = NULL, current_space = NULL;
1011:   PetscFreeSpaceList free_space_lvl = NULL, current_space_lvl = NULL;
1012:   PetscBT            lnkbt;

1014:   PetscFunctionBegin;
1015:   PetscCall(MatMissingDiagonal(A, &missing, &i));
1016:   PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry %" PetscInt_FMT, i);

1018:   if (bs > 1) {
1019:     if (!a->sbaijMat) PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &a->sbaijMat));
1020:     fact->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */

1022:     PetscCall(MatICCFactorSymbolic(fact, a->sbaijMat, perm, info));
1023:     PetscFunctionReturn(PETSC_SUCCESS);
1024:   }

1026:   PetscCall(ISIdentity(perm, &perm_identity));
1027:   PetscCall(ISGetIndices(perm, &rip));

1029:   /* special case that simply copies fill pattern */
1030:   if (!levels && perm_identity) {
1031:     PetscCall(PetscMalloc1(am + 1, &ui));
1032:     for (i = 0; i < am; i++) ui[i] = ai[i + 1] - a->diag[i]; /* ui: rowlengths - changes when !perm_identity */
1033:     B = fact;
1034:     PetscCall(MatSeqSBAIJSetPreallocation(B, 1, 0, ui));

1036:     b  = (Mat_SeqSBAIJ *)B->data;
1037:     uj = b->j;
1038:     for (i = 0; i < am; i++) {
1039:       aj = a->j + a->diag[i];
1040:       for (j = 0; j < ui[i]; j++) *uj++ = *aj++;
1041:       b->ilen[i] = ui[i];
1042:     }
1043:     PetscCall(PetscFree(ui));

1045:     B->factortype = MAT_FACTOR_NONE;

1047:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1048:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1049:     B->factortype = MAT_FACTOR_ICC;

1051:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1052:     PetscFunctionReturn(PETSC_SUCCESS);
1053:   }

1055:   /* initialization */
1056:   PetscCall(PetscMalloc1(am + 1, &ui));
1057:   ui[0] = 0;
1058:   PetscCall(PetscMalloc1(2 * am + 1, &cols_lvl));

1060:   /* jl: linked list for storing indices of the pivot rows
1061:      il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
1062:   PetscCall(PetscMalloc4(am, &uj_ptr, am, &uj_lvl_ptr, am, &il, am, &jl));
1063:   for (i = 0; i < am; i++) {
1064:     jl[i] = am;
1065:     il[i] = 0;
1066:   }

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

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

1075:   current_space = free_space;

1077:   PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am] / 2, am / 2)), &free_space_lvl));
1078:   current_space_lvl = free_space_lvl;

1080:   for (k = 0; k < am; k++) { /* for each active row k */
1081:     /* initialize lnk by the column indices of row rip[k] of A */
1082:     nzk         = 0;
1083:     ncols       = ai[rip[k] + 1] - ai[rip[k]];
1084:     ncols_upper = 0;
1085:     cols        = cols_lvl + am;
1086:     for (j = 0; j < ncols; j++) {
1087:       i = rip[*(aj + ai[rip[k]] + j)];
1088:       if (i >= k) { /* only take upper triangular entry */
1089:         cols[ncols_upper]     = i;
1090:         cols_lvl[ncols_upper] = -1; /* initialize level for nonzero entries */
1091:         ncols_upper++;
1092:       }
1093:     }
1094:     PetscCall(PetscIncompleteLLAdd(ncols_upper, cols, levels, cols_lvl, am, &nlnk, lnk, lnk_lvl, lnkbt));
1095:     nzk += nlnk;

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

1100:     while (prow < k) {
1101:       nextprow = jl[prow];

1103:       /* merge prow into k-th row */
1104:       jmin  = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
1105:       jmax  = ui[prow + 1];
1106:       ncols = jmax - jmin;
1107:       i     = jmin - ui[prow];
1108:       cols  = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
1109:       for (j = 0; j < ncols; j++) cols_lvl[j] = *(uj_lvl_ptr[prow] + i + j);
1110:       PetscCall(PetscIncompleteLLAddSorted(ncols, cols, levels, cols_lvl, am, &nlnk, lnk, lnk_lvl, lnkbt));
1111:       nzk += nlnk;

1113:       /* update il and jl for prow */
1114:       if (jmin < jmax) {
1115:         il[prow] = jmin;

1117:         j        = *cols;
1118:         jl[prow] = jl[j];
1119:         jl[j]    = prow;
1120:       }
1121:       prow = nextprow;
1122:     }

1124:     /* if free space is not available, make more free space */
1125:     if (current_space->local_remaining < nzk) {
1126:       i = am - k + 1;                                                             /* num of unfactored rows */
1127:       i = PetscMin(PetscIntMultTruncate(i, nzk), PetscIntMultTruncate(i, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1128:       PetscCall(PetscFreeSpaceGet(i, &current_space));
1129:       PetscCall(PetscFreeSpaceGet(i, &current_space_lvl));
1130:       reallocs++;
1131:     }

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

1136:     /* add the k-th row into il and jl */
1137:     if (nzk - 1 > 0) {
1138:       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
1139:       jl[k] = jl[i];
1140:       jl[i] = k;
1141:       il[k] = ui[k] + 1;
1142:     }
1143:     uj_ptr[k]     = current_space->array;
1144:     uj_lvl_ptr[k] = current_space_lvl->array;

1146:     current_space->array += nzk;
1147:     current_space->local_used += nzk;
1148:     current_space->local_remaining -= nzk;

1150:     current_space_lvl->array += nzk;
1151:     current_space_lvl->local_used += nzk;
1152:     current_space_lvl->local_remaining -= nzk;

1154:     ui[k + 1] = ui[k] + nzk;
1155:   }

1157:   PetscCall(ISRestoreIndices(perm, &rip));
1158:   PetscCall(PetscFree4(uj_ptr, uj_lvl_ptr, il, jl));
1159:   PetscCall(PetscFree(cols_lvl));

1161:   /* copy free_space into uj and free free_space; set uj in new datastructure; */
1162:   PetscCall(PetscMalloc1(ui[am] + 1, &uj));
1163:   PetscCall(PetscFreeSpaceContiguous(&free_space, uj));
1164:   PetscCall(PetscIncompleteLLDestroy(lnk, lnkbt));
1165:   PetscCall(PetscFreeSpaceDestroy(free_space_lvl));

1167:   /* put together the new matrix in MATSEQSBAIJ format */
1168:   B = fact;
1169:   PetscCall(MatSeqSBAIJSetPreallocation(B, 1, MAT_SKIP_ALLOCATION, NULL));

1171:   b          = (Mat_SeqSBAIJ *)B->data;
1172:   b->free_a  = PETSC_TRUE;
1173:   b->free_ij = PETSC_TRUE;

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

1177:   b->j             = uj;
1178:   b->i             = ui;
1179:   b->diag          = NULL;
1180:   b->ilen          = NULL;
1181:   b->imax          = NULL;
1182:   b->row           = perm;
1183:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */

1185:   PetscCall(PetscObjectReference((PetscObject)perm));

1187:   b->icol = perm;

1189:   PetscCall(PetscObjectReference((PetscObject)perm));
1190:   PetscCall(PetscMalloc1(am + 1, &b->solve_work));

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

1194:   B->info.factor_mallocs   = reallocs;
1195:   B->info.fill_ratio_given = fill;
1196:   if (ai[am] != 0.) {
1197:     /* nonzeros in lower triangular part of A (including diagonals)= (ai[am]+am)/2 */
1198:     B->info.fill_ratio_needed = ((PetscReal)2 * ui[am]) / (ai[am] + am);
1199:   } else {
1200:     B->info.fill_ratio_needed = 0.0;
1201:   }
1202: #if defined(PETSC_USE_INFO)
1203:   if (ai[am] != 0) {
1204:     PetscReal af = B->info.fill_ratio_needed;
1205:     PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
1206:     PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
1207:     PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
1208:   } else {
1209:     PetscCall(PetscInfo(A, "Empty matrix\n"));
1210:   }
1211: #endif
1212:   if (perm_identity) {
1213:     B->ops->solve                 = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1214:     B->ops->solvetranspose        = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1215:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1216:   } else {
1217:     fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
1218:   }
1219:   PetscFunctionReturn(PETSC_SUCCESS);
1220: }

1222: PetscErrorCode MatCholeskyFactorSymbolic_SeqBAIJ(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
1223: {
1224:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ *)A->data;
1225:   Mat_SeqSBAIJ      *b;
1226:   Mat                B;
1227:   PetscBool          perm_identity, missing;
1228:   PetscReal          fill = info->fill;
1229:   const PetscInt    *rip;
1230:   PetscInt           i, mbs = a->mbs, bs = A->rmap->bs, *ai = a->i, *aj = a->j, reallocs = 0, prow;
1231:   PetscInt          *jl, jmin, jmax, nzk, *ui, k, j, *il, nextprow;
1232:   PetscInt           nlnk, *lnk, ncols, ncols_upper, *cols, *uj, **ui_ptr, *uj_ptr;
1233:   PetscFreeSpaceList free_space = NULL, current_space = NULL;
1234:   PetscBT            lnkbt;

1236:   PetscFunctionBegin;
1237:   if (bs > 1) { /* convert to seqsbaij */
1238:     if (!a->sbaijMat) PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &a->sbaijMat));
1239:     fact->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */

1241:     PetscCall(MatCholeskyFactorSymbolic(fact, a->sbaijMat, perm, info));
1242:     PetscFunctionReturn(PETSC_SUCCESS);
1243:   }

1245:   PetscCall(MatMissingDiagonal(A, &missing, &i));
1246:   PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry %" PetscInt_FMT, i);

1248:   /* check whether perm is the identity mapping */
1249:   PetscCall(ISIdentity(perm, &perm_identity));
1250:   PetscCheck(perm_identity, PETSC_COMM_SELF, PETSC_ERR_SUP, "Matrix reordering is not supported");
1251:   PetscCall(ISGetIndices(perm, &rip));

1253:   /* initialization */
1254:   PetscCall(PetscMalloc1(mbs + 1, &ui));
1255:   ui[0] = 0;

1257:   /* jl: linked list for storing indices of the pivot rows
1258:      il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
1259:   PetscCall(PetscMalloc4(mbs, &ui_ptr, mbs, &il, mbs, &jl, mbs, &cols));
1260:   for (i = 0; i < mbs; i++) {
1261:     jl[i] = mbs;
1262:     il[i] = 0;
1263:   }

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

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

1272:   current_space = free_space;

1274:   for (k = 0; k < mbs; k++) { /* for each active row k */
1275:     /* initialize lnk by the column indices of row rip[k] of A */
1276:     nzk         = 0;
1277:     ncols       = ai[rip[k] + 1] - ai[rip[k]];
1278:     ncols_upper = 0;
1279:     for (j = 0; j < ncols; j++) {
1280:       i = rip[*(aj + ai[rip[k]] + j)];
1281:       if (i >= k) { /* only take upper triangular entry */
1282:         cols[ncols_upper] = i;
1283:         ncols_upper++;
1284:       }
1285:     }
1286:     PetscCall(PetscLLAdd(ncols_upper, cols, mbs, &nlnk, lnk, lnkbt));
1287:     nzk += nlnk;

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

1292:     while (prow < k) {
1293:       nextprow = jl[prow];
1294:       /* merge prow into k-th row */
1295:       jmin   = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
1296:       jmax   = ui[prow + 1];
1297:       ncols  = jmax - jmin;
1298:       uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
1299:       PetscCall(PetscLLAddSorted(ncols, uj_ptr, mbs, &nlnk, lnk, lnkbt));
1300:       nzk += nlnk;

1302:       /* update il and jl for prow */
1303:       if (jmin < jmax) {
1304:         il[prow] = jmin;
1305:         j        = *uj_ptr;
1306:         jl[prow] = jl[j];
1307:         jl[j]    = prow;
1308:       }
1309:       prow = nextprow;
1310:     }

1312:     /* if free space is not available, make more free space */
1313:     if (current_space->local_remaining < nzk) {
1314:       i = mbs - k + 1;                                                            /* num of unfactored rows */
1315:       i = PetscMin(PetscIntMultTruncate(i, nzk), PetscIntMultTruncate(i, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1316:       PetscCall(PetscFreeSpaceGet(i, &current_space));
1317:       reallocs++;
1318:     }

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

1323:     /* add the k-th row into il and jl */
1324:     if (nzk - 1 > 0) {
1325:       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
1326:       jl[k] = jl[i];
1327:       jl[i] = k;
1328:       il[k] = ui[k] + 1;
1329:     }
1330:     ui_ptr[k] = current_space->array;
1331:     current_space->array += nzk;
1332:     current_space->local_used += nzk;
1333:     current_space->local_remaining -= nzk;

1335:     ui[k + 1] = ui[k] + nzk;
1336:   }

1338:   PetscCall(ISRestoreIndices(perm, &rip));
1339:   PetscCall(PetscFree4(ui_ptr, il, jl, cols));

1341:   /* copy free_space into uj and free free_space; set uj in new datastructure; */
1342:   PetscCall(PetscMalloc1(ui[mbs] + 1, &uj));
1343:   PetscCall(PetscFreeSpaceContiguous(&free_space, uj));
1344:   PetscCall(PetscLLDestroy(lnk, lnkbt));

1346:   /* put together the new matrix in MATSEQSBAIJ format */
1347:   B = fact;
1348:   PetscCall(MatSeqSBAIJSetPreallocation(B, bs, MAT_SKIP_ALLOCATION, NULL));

1350:   b          = (Mat_SeqSBAIJ *)B->data;
1351:   b->free_a  = PETSC_TRUE;
1352:   b->free_ij = PETSC_TRUE;

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

1356:   b->j             = uj;
1357:   b->i             = ui;
1358:   b->diag          = NULL;
1359:   b->ilen          = NULL;
1360:   b->imax          = NULL;
1361:   b->row           = perm;
1362:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */

1364:   PetscCall(PetscObjectReference((PetscObject)perm));
1365:   b->icol = perm;
1366:   PetscCall(PetscObjectReference((PetscObject)perm));
1367:   PetscCall(PetscMalloc1(mbs + 1, &b->solve_work));
1368:   b->maxnz = b->nz = ui[mbs];

1370:   B->info.factor_mallocs   = reallocs;
1371:   B->info.fill_ratio_given = fill;
1372:   if (ai[mbs] != 0.) {
1373:     /* nonzeros in lower triangular part of A = (ai[mbs]+mbs)/2 */
1374:     B->info.fill_ratio_needed = ((PetscReal)2 * ui[mbs]) / (ai[mbs] + mbs);
1375:   } else {
1376:     B->info.fill_ratio_needed = 0.0;
1377:   }
1378: #if defined(PETSC_USE_INFO)
1379:   if (ai[mbs] != 0.) {
1380:     PetscReal af = B->info.fill_ratio_needed;
1381:     PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
1382:     PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
1383:     PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
1384:   } else {
1385:     PetscCall(PetscInfo(A, "Empty matrix\n"));
1386:   }
1387: #endif
1388:   if (perm_identity) {
1389:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1390:   } else {
1391:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
1392:   }
1393:   PetscFunctionReturn(PETSC_SUCCESS);
1394: }

1396: PetscErrorCode MatSolve_SeqBAIJ_N_NaturalOrdering(Mat A, Vec bb, Vec xx)
1397: {
1398:   Mat_SeqBAIJ       *a  = (Mat_SeqBAIJ *)A->data;
1399:   const PetscInt    *ai = a->i, *aj = a->j, *adiag = a->diag, *vi;
1400:   PetscInt           i, k, n                       = a->mbs;
1401:   PetscInt           nz, bs = A->rmap->bs, bs2 = a->bs2;
1402:   const MatScalar   *aa = a->a, *v;
1403:   PetscScalar       *x, *s, *t, *ls;
1404:   const PetscScalar *b;

1406:   PetscFunctionBegin;
1407:   PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected bs %" PetscInt_FMT " > 0", bs);
1408:   PetscCall(VecGetArrayRead(bb, &b));
1409:   PetscCall(VecGetArray(xx, &x));
1410:   t = a->solve_work;

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

1415:   for (i = 1; i < n; i++) {
1416:     v  = aa + bs2 * ai[i];
1417:     vi = aj + ai[i];
1418:     nz = ai[i + 1] - ai[i];
1419:     s  = t + bs * i;
1420:     PetscCall(PetscArraycpy(s, b + bs * i, bs)); /* copy i_th block of b to t */
1421:     for (k = 0; k < nz; k++) {
1422:       PetscKernel_v_gets_v_minus_A_times_w(bs, s, v, t + bs * vi[k]);
1423:       v += bs2;
1424:     }
1425:   }

1427:   /* backward solve the upper triangular */
1428:   ls = a->solve_work + A->cmap->n;
1429:   for (i = n - 1; i >= 0; i--) {
1430:     v  = aa + bs2 * (adiag[i + 1] + 1);
1431:     vi = aj + adiag[i + 1] + 1;
1432:     nz = adiag[i] - adiag[i + 1] - 1;
1433:     PetscCall(PetscArraycpy(ls, t + i * bs, bs));
1434:     for (k = 0; k < nz; k++) {
1435:       PetscKernel_v_gets_v_minus_A_times_w(bs, ls, v, t + bs * vi[k]);
1436:       v += bs2;
1437:     }
1438:     PetscKernel_w_gets_A_times_v(bs, ls, aa + bs2 * adiag[i], t + i * bs); /* *inv(diagonal[i]) */
1439:     PetscCall(PetscArraycpy(x + i * bs, t + i * bs, bs));
1440:   }

1442:   PetscCall(VecRestoreArrayRead(bb, &b));
1443:   PetscCall(VecRestoreArray(xx, &x));
1444:   PetscCall(PetscLogFlops(2.0 * (a->bs2) * (a->nz) - A->rmap->bs * A->cmap->n));
1445:   PetscFunctionReturn(PETSC_SUCCESS);
1446: }

1448: PetscErrorCode MatSolve_SeqBAIJ_N(Mat A, Vec bb, Vec xx)
1449: {
1450:   Mat_SeqBAIJ       *a     = (Mat_SeqBAIJ *)A->data;
1451:   IS                 iscol = a->col, isrow = a->row;
1452:   const PetscInt    *r, *c, *rout, *cout, *ai = a->i, *aj = a->j, *adiag = a->diag, *vi;
1453:   PetscInt           i, m, n = a->mbs;
1454:   PetscInt           nz, bs = A->rmap->bs, bs2 = a->bs2;
1455:   const MatScalar   *aa = a->a, *v;
1456:   PetscScalar       *x, *s, *t, *ls;
1457:   const PetscScalar *b;

1459:   PetscFunctionBegin;
1460:   PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected bs %" PetscInt_FMT " > 0", bs);
1461:   PetscCall(VecGetArrayRead(bb, &b));
1462:   PetscCall(VecGetArray(xx, &x));
1463:   t = a->solve_work;

1465:   PetscCall(ISGetIndices(isrow, &rout));
1466:   r = rout;
1467:   PetscCall(ISGetIndices(iscol, &cout));
1468:   c = cout;

1470:   /* forward solve the lower triangular */
1471:   PetscCall(PetscArraycpy(t, b + bs * r[0], bs));
1472:   for (i = 1; i < n; i++) {
1473:     v  = aa + bs2 * ai[i];
1474:     vi = aj + ai[i];
1475:     nz = ai[i + 1] - ai[i];
1476:     s  = t + bs * i;
1477:     PetscCall(PetscArraycpy(s, b + bs * r[i], bs));
1478:     for (m = 0; m < nz; m++) {
1479:       PetscKernel_v_gets_v_minus_A_times_w(bs, s, v, t + bs * vi[m]);
1480:       v += bs2;
1481:     }
1482:   }

1484:   /* backward solve the upper triangular */
1485:   ls = a->solve_work + A->cmap->n;
1486:   for (i = n - 1; i >= 0; i--) {
1487:     v  = aa + bs2 * (adiag[i + 1] + 1);
1488:     vi = aj + adiag[i + 1] + 1;
1489:     nz = adiag[i] - adiag[i + 1] - 1;
1490:     PetscCall(PetscArraycpy(ls, t + i * bs, bs));
1491:     for (m = 0; m < nz; m++) {
1492:       PetscKernel_v_gets_v_minus_A_times_w(bs, ls, v, t + bs * vi[m]);
1493:       v += bs2;
1494:     }
1495:     PetscKernel_w_gets_A_times_v(bs, ls, v, t + i * bs); /* *inv(diagonal[i]) */
1496:     PetscCall(PetscArraycpy(x + bs * c[i], t + i * bs, bs));
1497:   }
1498:   PetscCall(ISRestoreIndices(isrow, &rout));
1499:   PetscCall(ISRestoreIndices(iscol, &cout));
1500:   PetscCall(VecRestoreArrayRead(bb, &b));
1501:   PetscCall(VecRestoreArray(xx, &x));
1502:   PetscCall(PetscLogFlops(2.0 * (a->bs2) * (a->nz) - A->rmap->bs * A->cmap->n));
1503:   PetscFunctionReturn(PETSC_SUCCESS);
1504: }

1506: /*
1507:     For each block in an block array saves the largest absolute value in the block into another array
1508: */
1509: static PetscErrorCode MatBlockAbs_private(PetscInt nbs, PetscInt bs2, PetscScalar *blockarray, PetscReal *absarray)
1510: {
1511:   PetscInt i, j;

1513:   PetscFunctionBegin;
1514:   PetscCall(PetscArrayzero(absarray, nbs + 1));
1515:   for (i = 0; i < nbs; i++) {
1516:     for (j = 0; j < bs2; j++) {
1517:       if (absarray[i] < PetscAbsScalar(blockarray[i * nbs + j])) absarray[i] = PetscAbsScalar(blockarray[i * nbs + j]);
1518:     }
1519:   }
1520:   PetscFunctionReturn(PETSC_SUCCESS);
1521: }

1523: /*
1524:      This needs to be renamed and called by the regular MatILUFactor_SeqBAIJ when drop tolerance is used
1525: */
1526: PetscErrorCode MatILUDTFactor_SeqBAIJ(Mat A, IS isrow, IS iscol, const MatFactorInfo *info, Mat *fact)
1527: {
1528:   Mat             B = *fact;
1529:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ *)A->data, *b;
1530:   IS              isicol;
1531:   const PetscInt *r, *ic;
1532:   PetscInt        i, mbs = a->mbs, bs = A->rmap->bs, bs2 = a->bs2, *ai = a->i, *aj = a->j, *ajtmp, *adiag;
1533:   PetscInt       *bi, *bj, *bdiag;

1535:   PetscInt   row, nzi, nzi_bl, nzi_bu, *im, dtcount, nzi_al, nzi_au;
1536:   PetscInt   nlnk, *lnk;
1537:   PetscBT    lnkbt;
1538:   PetscBool  row_identity, icol_identity;
1539:   MatScalar *aatmp, *pv, *batmp, *ba, *rtmp, *pc, *multiplier, *vtmp;
1540:   PetscInt   j, nz, *pj, *bjtmp, k, ncut, *jtmp;

1542:   PetscReal  dt = info->dt; /* shift=info->shiftamount; */
1543:   PetscInt   nnz_max;
1544:   PetscBool  missing;
1545:   PetscReal *vtmp_abs;
1546:   MatScalar *v_work;
1547:   PetscInt  *v_pivots;
1548:   PetscBool  allowzeropivot, zeropivotdetected = PETSC_FALSE;

1550:   PetscFunctionBegin;
1551:   /* symbolic factorization, can be reused */
1552:   PetscCall(MatMissingDiagonal(A, &missing, &i));
1553:   PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry %" PetscInt_FMT, i);
1554:   adiag = a->diag;

1556:   PetscCall(ISInvertPermutation(iscol, PETSC_DECIDE, &isicol));

1558:   /* bdiag is location of diagonal in factor */
1559:   PetscCall(PetscMalloc1(mbs + 1, &bdiag));

1561:   /* allocate row pointers bi */
1562:   PetscCall(PetscMalloc1(2 * mbs + 2, &bi));

1564:   /* allocate bj and ba; max num of nonzero entries is (ai[n]+2*n*dtcount+2) */
1565:   dtcount = (PetscInt)info->dtcount;
1566:   if (dtcount > mbs - 1) dtcount = mbs - 1;
1567:   nnz_max = ai[mbs] + 2 * mbs * dtcount + 2;
1568:   /* printf("MatILUDTFactor_SeqBAIJ, bs %d, ai[mbs] %d, nnz_max  %d, dtcount %d\n",bs,ai[mbs],nnz_max,dtcount); */
1569:   PetscCall(PetscMalloc1(nnz_max, &bj));
1570:   nnz_max = nnz_max * bs2;
1571:   PetscCall(PetscMalloc1(nnz_max, &ba));

1573:   /* put together the new matrix */
1574:   PetscCall(MatSeqBAIJSetPreallocation(B, bs, MAT_SKIP_ALLOCATION, NULL));

1576:   b          = (Mat_SeqBAIJ *)B->data;
1577:   b->free_a  = PETSC_TRUE;
1578:   b->free_ij = PETSC_TRUE;

1580:   b->a    = ba;
1581:   b->j    = bj;
1582:   b->i    = bi;
1583:   b->diag = bdiag;
1584:   b->ilen = NULL;
1585:   b->imax = NULL;
1586:   b->row  = isrow;
1587:   b->col  = iscol;

1589:   PetscCall(PetscObjectReference((PetscObject)isrow));
1590:   PetscCall(PetscObjectReference((PetscObject)iscol));

1592:   b->icol = isicol;
1593:   PetscCall(PetscMalloc1(bs * (mbs + 1), &b->solve_work));
1594:   b->maxnz = nnz_max / bs2;

1596:   B->factortype            = MAT_FACTOR_ILUDT;
1597:   B->info.factor_mallocs   = 0;
1598:   B->info.fill_ratio_given = ((PetscReal)nnz_max) / ((PetscReal)(ai[mbs] * bs2));
1599:   /* end of symbolic factorization */
1600:   PetscCall(ISGetIndices(isrow, &r));
1601:   PetscCall(ISGetIndices(isicol, &ic));

1603:   /* linked list for storing column indices of the active row */
1604:   nlnk = mbs + 1;
1605:   PetscCall(PetscLLCreate(mbs, mbs, nlnk, lnk, lnkbt));

1607:   /* im: used by PetscLLAddSortedLU(); jtmp: working array for column indices of active row */
1608:   PetscCall(PetscMalloc2(mbs, &im, mbs, &jtmp));
1609:   /* rtmp, vtmp: working arrays for sparse and contiguous row entries of active row */
1610:   PetscCall(PetscMalloc2(mbs * bs2, &rtmp, mbs * bs2, &vtmp));
1611:   PetscCall(PetscMalloc1(mbs + 1, &vtmp_abs));
1612:   PetscCall(PetscMalloc3(bs, &v_work, bs2, &multiplier, bs, &v_pivots));

1614:   allowzeropivot  = PetscNot(A->erroriffailure);
1615:   bi[0]           = 0;
1616:   bdiag[0]        = (nnz_max / bs2) - 1; /* location of diagonal in factor B */
1617:   bi[2 * mbs + 1] = bdiag[0] + 1;        /* endof bj and ba array */
1618:   for (i = 0; i < mbs; i++) {
1619:     /* copy initial fill into linked list */
1620:     nzi = ai[r[i] + 1] - ai[r[i]];
1621:     PetscCheck(nzi, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Empty row in matrix: row in original ordering %" PetscInt_FMT " in permuted ordering %" PetscInt_FMT, r[i], i);
1622:     nzi_al = adiag[r[i]] - ai[r[i]];
1623:     nzi_au = ai[r[i] + 1] - adiag[r[i]] - 1;

1625:     /* load in initial unfactored row */
1626:     ajtmp = aj + ai[r[i]];
1627:     PetscCall(PetscLLAddPerm(nzi, ajtmp, ic, mbs, &nlnk, lnk, lnkbt));
1628:     PetscCall(PetscArrayzero(rtmp, mbs * bs2));
1629:     aatmp = a->a + bs2 * ai[r[i]];
1630:     for (j = 0; j < nzi; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ic[ajtmp[j]], aatmp + bs2 * j, bs2));

1632:     /* add pivot rows into linked list */
1633:     row = lnk[mbs];
1634:     while (row < i) {
1635:       nzi_bl = bi[row + 1] - bi[row] + 1;
1636:       bjtmp  = bj + bdiag[row + 1] + 1; /* points to 1st column next to the diagonal in U */
1637:       PetscCall(PetscLLAddSortedLU(bjtmp, row, &nlnk, lnk, lnkbt, i, nzi_bl, im));
1638:       nzi += nlnk;
1639:       row = lnk[row];
1640:     }

1642:     /* copy data from lnk into jtmp, then initialize lnk */
1643:     PetscCall(PetscLLClean(mbs, mbs, nzi, lnk, jtmp, lnkbt));

1645:     /* numerical factorization */
1646:     bjtmp = jtmp;
1647:     row   = *bjtmp++; /* 1st pivot row */

1649:     while (row < i) {
1650:       pc = rtmp + bs2 * row;
1651:       pv = ba + bs2 * bdiag[row];                           /* inv(diag) of the pivot row */
1652:       PetscKernel_A_gets_A_times_B(bs, pc, pv, multiplier); /* pc= multiplier = pc*inv(diag[row]) */
1653:       PetscCall(MatBlockAbs_private(1, bs2, pc, vtmp_abs));
1654:       if (vtmp_abs[0] > dt) {         /* apply tolerance dropping rule */
1655:         pj = bj + bdiag[row + 1] + 1; /* point to 1st entry of U(row,:) */
1656:         pv = ba + bs2 * (bdiag[row + 1] + 1);
1657:         nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries in U(row,:), excluding diagonal */
1658:         for (j = 0; j < nz; j++) PetscKernel_A_gets_A_minus_B_times_C(bs, rtmp + bs2 * pj[j], pc, pv + bs2 * j);
1659:         /* PetscCall(PetscLogFlops(bslog*(nz+1.0)-bs)); */
1660:       }
1661:       row = *bjtmp++;
1662:     }

1664:     /* copy sparse rtmp into contiguous vtmp; separate L and U part */
1665:     nzi_bl = 0;
1666:     j      = 0;
1667:     while (jtmp[j] < i) { /* L-part. Note: jtmp is sorted */
1668:       PetscCall(PetscArraycpy(vtmp + bs2 * j, rtmp + bs2 * jtmp[j], bs2));
1669:       nzi_bl++;
1670:       j++;
1671:     }
1672:     nzi_bu = nzi - nzi_bl - 1;

1674:     while (j < nzi) { /* U-part */
1675:       PetscCall(PetscArraycpy(vtmp + bs2 * j, rtmp + bs2 * jtmp[j], bs2));
1676:       j++;
1677:     }

1679:     PetscCall(MatBlockAbs_private(nzi, bs2, vtmp, vtmp_abs));

1681:     bjtmp = bj + bi[i];
1682:     batmp = ba + bs2 * bi[i];
1683:     /* apply level dropping rule to L part */
1684:     ncut = nzi_al + dtcount;
1685:     if (ncut < nzi_bl) {
1686:       PetscCall(PetscSortSplitReal(ncut, nzi_bl, vtmp_abs, jtmp));
1687:       PetscCall(PetscSortIntWithScalarArray(ncut, jtmp, vtmp));
1688:     } else {
1689:       ncut = nzi_bl;
1690:     }
1691:     for (j = 0; j < ncut; j++) {
1692:       bjtmp[j] = jtmp[j];
1693:       PetscCall(PetscArraycpy(batmp + bs2 * j, rtmp + bs2 * bjtmp[j], bs2));
1694:     }
1695:     bi[i + 1] = bi[i] + ncut;
1696:     nzi       = ncut + 1;

1698:     /* apply level dropping rule to U part */
1699:     ncut = nzi_au + dtcount;
1700:     if (ncut < nzi_bu) {
1701:       PetscCall(PetscSortSplitReal(ncut, nzi_bu, vtmp_abs + nzi_bl + 1, jtmp + nzi_bl + 1));
1702:       PetscCall(PetscSortIntWithScalarArray(ncut, jtmp + nzi_bl + 1, vtmp + nzi_bl + 1));
1703:     } else {
1704:       ncut = nzi_bu;
1705:     }
1706:     nzi += ncut;

1708:     /* mark bdiagonal */
1709:     bdiag[i + 1]    = bdiag[i] - (ncut + 1);
1710:     bi[2 * mbs - i] = bi[2 * mbs - i + 1] - (ncut + 1);

1712:     bjtmp = bj + bdiag[i];
1713:     batmp = ba + bs2 * bdiag[i];
1714:     PetscCall(PetscArraycpy(batmp, rtmp + bs2 * i, bs2));
1715:     *bjtmp = i;

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

1720:     for (k = 0; k < ncut; k++) {
1721:       bjtmp[k] = jtmp[nzi_bl + 1 + k];
1722:       PetscCall(PetscArraycpy(batmp + bs2 * k, rtmp + bs2 * bjtmp[k], bs2));
1723:     }

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

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

1730:     PetscCall(PetscKernel_A_gets_inverse_A(bs, batmp, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
1731:     if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1732:   } /* for (i=0; i<mbs; i++) */
1733:   PetscCall(PetscFree3(v_work, multiplier, v_pivots));

1735:   /* printf("end of L %d, beginning of U %d\n",bi[mbs],bdiag[mbs]); */
1736:   PetscCheck(bi[mbs] < bdiag[mbs], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "end of L array %" PetscInt_FMT " cannot >= the beginning of U array %" PetscInt_FMT, bi[mbs], bdiag[mbs]);

1738:   PetscCall(ISRestoreIndices(isrow, &r));
1739:   PetscCall(ISRestoreIndices(isicol, &ic));

1741:   PetscCall(PetscLLDestroy(lnk, lnkbt));

1743:   PetscCall(PetscFree2(im, jtmp));
1744:   PetscCall(PetscFree2(rtmp, vtmp));

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

1749:   PetscCall(ISIdentity(isrow, &row_identity));
1750:   PetscCall(ISIdentity(isicol, &icol_identity));
1751:   if (row_identity && icol_identity) {
1752:     B->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
1753:   } else {
1754:     B->ops->solve = MatSolve_SeqBAIJ_N;
1755:   }

1757:   B->ops->solveadd          = NULL;
1758:   B->ops->solvetranspose    = NULL;
1759:   B->ops->solvetransposeadd = NULL;
1760:   B->ops->matsolve          = NULL;
1761:   B->assembled              = PETSC_TRUE;
1762:   B->preallocated           = PETSC_TRUE;
1763:   PetscFunctionReturn(PETSC_SUCCESS);
1764: }