Actual source code: baijfact2.c


  2: /*
  3:     Factorization code for BAIJ format.
  4: */

  6: #include <../src/mat/impls/baij/seq/baij.h>
  7: #include <petsc/private/kernels/blockinvert.h>
  8: #include <petscbt.h>
  9: #include <../src/mat/utils/freespace.h>

 11: /* ----------------------------------------------------------------*/
 12: extern PetscErrorCode MatDuplicateNoCreate_SeqBAIJ(Mat, Mat, MatDuplicateOption, PetscBool);

 14: /*
 15:    This is not much faster than MatLUFactorNumeric_SeqBAIJ_N() but the solve is faster at least sometimes
 16: */
 17: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_15_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info)
 18: {
 19:   Mat              C = B;
 20:   Mat_SeqBAIJ     *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
 21:   PetscInt         i, j, k, ipvt[15];
 22:   const PetscInt   n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *ajtmp, *bjtmp, *bdiag = b->diag, *pj;
 23:   PetscInt         nz, nzL, row;
 24:   MatScalar       *rtmp, *pc, *mwork, *pv, *vv, work[225];
 25:   const MatScalar *v, *aa = a->a;
 26:   PetscInt         bs2 = a->bs2, bs = A->rmap->bs, flg;
 27:   PetscInt         sol_ver;
 28:   PetscBool        allowzeropivot, zeropivotdetected;

 30:   allowzeropivot = PetscNot(A->erroriffailure);
 31:   PetscOptionsGetInt(NULL, ((PetscObject)A)->prefix, "-sol_ver", &sol_ver, NULL);

 33:   /* generate work space needed by the factorization */
 34:   PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork);
 35:   PetscArrayzero(rtmp, bs2 * n);

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

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

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

 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] != 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);
 70:         /* PetscKernel_A_gets_A_times_B_15(pc,pv,mwork); */
 71:         pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
 72:         pv = b->a + bs2 * (bdiag[row + 1] + 1);
 73:         nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
 74:         for (j = 0; j < nz; j++) {
 75:           vv = rtmp + bs2 * pj[j];
 76:           PetscKernel_A_gets_A_minus_B_times_C(bs, vv, pc, pv);
 77:           /* PetscKernel_A_gets_A_minus_B_times_C_15(vv,pc,pv); */
 78:           pv += bs2;
 79:         }
 80:         PetscLogFlops(2.0 * bs2 * bs * (nz + 1) - bs2); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
 81:       }
 82:     }

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

 91:     /* Mark diagonal and invert diagonal for simpler triangular solves */
 92:     pv = b->a + bs2 * bdiag[i];
 93:     pj = b->j + bdiag[i];
 94:     PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2);
 95:     PetscKernel_A_gets_inverse_A_15(pv, ipvt, work, info->shiftamount, allowzeropivot, &zeropivotdetected);
 96:     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

 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 (j = 0; j < nz; j++) PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2);
103:   }

105:   PetscFree2(rtmp, mwork);

107:   C->ops->solve          = MatSolve_SeqBAIJ_15_NaturalOrdering_ver1;
108:   C->ops->solvetranspose = MatSolve_SeqBAIJ_N_NaturalOrdering;
109:   C->assembled           = PETSC_TRUE;

111:   PetscLogFlops(1.333333333333 * bs * bs2 * b->mbs); /* from inverting diagonal blocks */
112:   return 0;
113: }

115: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_N(Mat B, Mat A, const MatFactorInfo *info)
116: {
117:   Mat             C = B;
118:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
119:   IS              isrow = b->row, isicol = b->icol;
120:   const PetscInt *r, *ic;
121:   PetscInt        i, j, k, n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
122:   PetscInt       *ajtmp, *bjtmp, nz, nzL, row, *bdiag = b->diag, *pj;
123:   MatScalar      *rtmp, *pc, *mwork, *v, *pv, *aa     = a->a;
124:   PetscInt        bs = A->rmap->bs, bs2 = a->bs2, *v_pivots, flg;
125:   MatScalar      *v_work;
126:   PetscBool       col_identity, row_identity, both_identity;
127:   PetscBool       allowzeropivot, zeropivotdetected;

129:   ISGetIndices(isrow, &r);
130:   ISGetIndices(isicol, &ic);
131:   allowzeropivot = PetscNot(A->erroriffailure);

133:   PetscCalloc1(bs2 * n, &rtmp);

135:   /* generate work space needed by dense LU factorization */
136:   PetscMalloc3(bs, &v_work, bs2, &mwork, bs, &v_pivots);

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

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

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

156:     /* elimination */
157:     bjtmp = bj + bi[i];
158:     nzL   = bi[i + 1] - bi[i];
159:     for (k = 0; k < nzL; k++) {
160:       row = bjtmp[k];
161:       pc  = rtmp + bs2 * row;
162:       for (flg = 0, j = 0; j < bs2; j++) {
163:         if (pc[j] != 0.0) {
164:           flg = 1;
165:           break;
166:         }
167:       }
168:       if (flg) {
169:         pv = b->a + bs2 * bdiag[row];
170:         PetscKernel_A_gets_A_times_B(bs, pc, pv, mwork); /* *pc = *pc * (*pv); */
171:         pj = b->j + bdiag[row + 1] + 1;                  /* beginning of U(row,:) */
172:         pv = b->a + bs2 * (bdiag[row + 1] + 1);
173:         nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
174:         for (j = 0; j < nz; j++) PetscKernel_A_gets_A_minus_B_times_C(bs, rtmp + bs2 * pj[j], pc, pv + bs2 * j);
175:         PetscLogFlops(2.0 * bs2 * bs * (nz + 1) - bs2); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
176:       }
177:     }

179:     /* finished row so stick it into b->a */
180:     /* L part */
181:     pv = b->a + bs2 * bi[i];
182:     pj = b->j + bi[i];
183:     nz = bi[i + 1] - bi[i];
184:     for (j = 0; j < nz; j++) PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2);

186:     /* Mark diagonal and invert diagonal for simpler triangular solves */
187:     pv = b->a + bs2 * bdiag[i];
188:     pj = b->j + bdiag[i];
189:     PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2);

191:     PetscKernel_A_gets_inverse_A(bs, pv, v_pivots, v_work, allowzeropivot, &zeropivotdetected);
192:     if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

194:     /* U part */
195:     pv = b->a + bs2 * (bdiag[i + 1] + 1);
196:     pj = b->j + bdiag[i + 1] + 1;
197:     nz = bdiag[i] - bdiag[i + 1] - 1;
198:     for (j = 0; j < nz; j++) PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2);
199:   }

201:   PetscFree(rtmp);
202:   PetscFree3(v_work, mwork, v_pivots);
203:   ISRestoreIndices(isicol, &ic);
204:   ISRestoreIndices(isrow, &r);

206:   ISIdentity(isrow, &row_identity);
207:   ISIdentity(isicol, &col_identity);

209:   both_identity = (PetscBool)(row_identity && col_identity);
210:   if (both_identity) {
211:     switch (bs) {
212:     case 9:
213: #if defined(PETSC_HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && defined(PETSC_USE_REAL_DOUBLE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_64BIT_INDICES)
214:       C->ops->solve = MatSolve_SeqBAIJ_9_NaturalOrdering;
215: #else
216:       C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
217: #endif
218:       break;
219:     case 11:
220:       C->ops->solve = MatSolve_SeqBAIJ_11_NaturalOrdering;
221:       break;
222:     case 12:
223:       C->ops->solve = MatSolve_SeqBAIJ_12_NaturalOrdering;
224:       break;
225:     case 13:
226:       C->ops->solve = MatSolve_SeqBAIJ_13_NaturalOrdering;
227:       break;
228:     case 14:
229:       C->ops->solve = MatSolve_SeqBAIJ_14_NaturalOrdering;
230:       break;
231:     default:
232:       C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
233:       break;
234:     }
235:   } else {
236:     C->ops->solve = MatSolve_SeqBAIJ_N;
237:   }
238:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N;

240:   C->assembled = PETSC_TRUE;

242:   PetscLogFlops(1.333333333333 * bs * bs2 * b->mbs); /* from inverting diagonal blocks */
243:   return 0;
244: }

246: /*
247:    ilu(0) with natural ordering under new data structure.
248:    See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed description
249:    because this code is almost identical to MatILUFactorSymbolic_SeqAIJ_ilu0_inplace().
250: */

252: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_ilu0(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
253: {
254:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b;
255:   PetscInt     n = a->mbs, *ai = a->i, *aj, *adiag = a->diag, bs2 = a->bs2;
256:   PetscInt     i, j, nz, *bi, *bj, *bdiag, bi_temp;

258:   MatDuplicateNoCreate_SeqBAIJ(fact, A, MAT_DO_NOT_COPY_VALUES, PETSC_FALSE);
259:   b = (Mat_SeqBAIJ *)(fact)->data;

261:   /* allocate matrix arrays for new data structure */
262:   PetscMalloc3(bs2 * ai[n] + 1, &b->a, ai[n] + 1, &b->j, n + 1, &b->i);

264:   b->singlemalloc    = PETSC_TRUE;
265:   b->free_a          = PETSC_TRUE;
266:   b->free_ij         = PETSC_TRUE;
267:   fact->preallocated = PETSC_TRUE;
268:   fact->assembled    = PETSC_TRUE;
269:   if (!b->diag) { PetscMalloc1(n + 1, &b->diag); }
270:   bdiag = b->diag;

272:   if (n > 0) PetscArrayzero(b->a, bs2 * ai[n]);

274:   /* set bi and bj with new data structure */
275:   bi = b->i;
276:   bj = b->j;

278:   /* L part */
279:   bi[0] = 0;
280:   for (i = 0; i < n; i++) {
281:     nz        = adiag[i] - ai[i];
282:     bi[i + 1] = bi[i] + nz;
283:     aj        = a->j + ai[i];
284:     for (j = 0; j < nz; j++) {
285:       *bj = aj[j];
286:       bj++;
287:     }
288:   }

290:   /* U part */
291:   bi_temp  = bi[n];
292:   bdiag[n] = bi[n] - 1;
293:   for (i = n - 1; i >= 0; i--) {
294:     nz      = ai[i + 1] - adiag[i] - 1;
295:     bi_temp = bi_temp + nz + 1;
296:     aj      = a->j + adiag[i] + 1;
297:     for (j = 0; j < nz; j++) {
298:       *bj = aj[j];
299:       bj++;
300:     }
301:     /* diag[i] */
302:     *bj = i;
303:     bj++;
304:     bdiag[i] = bi_temp - 1;
305:   }
306:   return 0;
307: }

309: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
310: {
311:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ *)A->data, *b;
312:   IS                 isicol;
313:   const PetscInt    *r, *ic;
314:   PetscInt           n = a->mbs, *ai = a->i, *aj = a->j, d;
315:   PetscInt          *bi, *cols, nnz, *cols_lvl;
316:   PetscInt          *bdiag, prow, fm, nzbd, reallocs = 0, dcount = 0;
317:   PetscInt           i, levels, diagonal_fill;
318:   PetscBool          col_identity, row_identity, both_identity;
319:   PetscReal          f;
320:   PetscInt           nlnk, *lnk, *lnk_lvl = NULL;
321:   PetscBT            lnkbt;
322:   PetscInt           nzi, *bj, **bj_ptr, **bjlvl_ptr;
323:   PetscFreeSpaceList free_space = NULL, current_space = NULL;
324:   PetscFreeSpaceList free_space_lvl = NULL, current_space_lvl = NULL;
325:   PetscBool          missing;
326:   PetscInt           bs = A->rmap->bs, bs2 = a->bs2;

329:   if (bs > 1) { /* check shifttype */
331:   }

333:   MatMissingDiagonal(A, &missing, &d);

336:   f             = info->fill;
337:   levels        = (PetscInt)info->levels;
338:   diagonal_fill = (PetscInt)info->diagonal_fill;

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

342:   ISIdentity(isrow, &row_identity);
343:   ISIdentity(iscol, &col_identity);

345:   both_identity = (PetscBool)(row_identity && col_identity);

347:   if (!levels && both_identity) {
348:     /* special case: ilu(0) with natural ordering */
349:     MatILUFactorSymbolic_SeqBAIJ_ilu0(fact, A, isrow, iscol, info);
350:     MatSeqBAIJSetNumericFactorization(fact, both_identity);

352:     fact->factortype               = MAT_FACTOR_ILU;
353:     (fact)->info.factor_mallocs    = 0;
354:     (fact)->info.fill_ratio_given  = info->fill;
355:     (fact)->info.fill_ratio_needed = 1.0;

357:     b       = (Mat_SeqBAIJ *)(fact)->data;
358:     b->row  = isrow;
359:     b->col  = iscol;
360:     b->icol = isicol;
361:     PetscObjectReference((PetscObject)isrow);
362:     PetscObjectReference((PetscObject)iscol);
363:     b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;

365:     PetscMalloc1((n + 1) * bs, &b->solve_work);
366:     return 0;
367:   }

369:   ISGetIndices(isrow, &r);
370:   ISGetIndices(isicol, &ic);

372:   /* get new row pointers */
373:   PetscMalloc1(n + 1, &bi);
374:   bi[0] = 0;
375:   /* bdiag is location of diagonal in factor */
376:   PetscMalloc1(n + 1, &bdiag);
377:   bdiag[0] = 0;

379:   PetscMalloc2(n, &bj_ptr, n, &bjlvl_ptr);

381:   /* create a linked list for storing column indices of the active row */
382:   nlnk = n + 1;
383:   PetscIncompleteLLCreate(n, n, nlnk, lnk, lnk_lvl, lnkbt);

385:   /* initial FreeSpace size is f*(ai[n]+1) */
386:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space);
387:   current_space = free_space;
388:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space_lvl);
389:   current_space_lvl = free_space_lvl;

391:   for (i = 0; i < n; i++) {
392:     nzi = 0;
393:     /* copy current row into linked list */
394:     nnz = ai[r[i] + 1] - ai[r[i]];
396:     cols   = aj + ai[r[i]];
397:     lnk[i] = -1; /* marker to indicate if diagonal exists */
398:     PetscIncompleteLLInit(nnz, cols, n, ic, &nlnk, lnk, lnk_lvl, lnkbt);
399:     nzi += nlnk;

401:     /* make sure diagonal entry is included */
402:     if (diagonal_fill && lnk[i] == -1) {
403:       fm = n;
404:       while (lnk[fm] < i) fm = lnk[fm];
405:       lnk[i]     = lnk[fm]; /* insert diagonal into linked list */
406:       lnk[fm]    = i;
407:       lnk_lvl[i] = 0;
408:       nzi++;
409:       dcount++;
410:     }

412:     /* add pivot rows into the active row */
413:     nzbd = 0;
414:     prow = lnk[n];
415:     while (prow < i) {
416:       nnz      = bdiag[prow];
417:       cols     = bj_ptr[prow] + nnz + 1;
418:       cols_lvl = bjlvl_ptr[prow] + nnz + 1;
419:       nnz      = bi[prow + 1] - bi[prow] - nnz - 1;

421:       PetscILULLAddSorted(nnz, cols, levels, cols_lvl, prow, &nlnk, lnk, lnk_lvl, lnkbt, prow);
422:       nzi += nlnk;
423:       prow = lnk[prow];
424:       nzbd++;
425:     }
426:     bdiag[i]  = nzbd;
427:     bi[i + 1] = bi[i] + nzi;

429:     /* if free space is not available, make more free space */
430:     if (current_space->local_remaining < nzi) {
431:       nnz = PetscIntMultTruncate(2, PetscIntMultTruncate(nzi, (n - i))); /* estimated and max additional space needed */
432:       PetscFreeSpaceGet(nnz, &current_space);
433:       PetscFreeSpaceGet(nnz, &current_space_lvl);
434:       reallocs++;
435:     }

437:     /* copy data into free_space and free_space_lvl, then initialize lnk */
438:     PetscIncompleteLLClean(n, n, nzi, lnk, lnk_lvl, current_space->array, current_space_lvl->array, lnkbt);

440:     bj_ptr[i]    = current_space->array;
441:     bjlvl_ptr[i] = current_space_lvl->array;

443:     /* make sure the active row i has diagonal entry */

446:     current_space->array += nzi;
447:     current_space->local_used += nzi;
448:     current_space->local_remaining -= nzi;

450:     current_space_lvl->array += nzi;
451:     current_space_lvl->local_used += nzi;
452:     current_space_lvl->local_remaining -= nzi;
453:   }

455:   ISRestoreIndices(isrow, &r);
456:   ISRestoreIndices(isicol, &ic);

458:   /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
459:   PetscMalloc1(bi[n] + 1, &bj);
460:   PetscFreeSpaceContiguous_LU(&free_space, bj, n, bi, bdiag);

462:   PetscIncompleteLLDestroy(lnk, lnkbt);
463:   PetscFreeSpaceDestroy(free_space_lvl);
464:   PetscFree2(bj_ptr, bjlvl_ptr);

466: #if defined(PETSC_USE_INFO)
467:   {
468:     PetscReal af = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);
469:     PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)f, (double)af);
470:     PetscInfo(A, "Run with -[sub_]pc_factor_fill %g or use \n", (double)af);
471:     PetscInfo(A, "PCFactorSetFill([sub]pc,%g);\n", (double)af);
472:     PetscInfo(A, "for best performance.\n");
473:     if (diagonal_fill) PetscInfo(A, "Detected and replaced %" PetscInt_FMT " missing diagonals\n", dcount);
474:   }
475: #endif

477:   /* put together the new matrix */
478:   MatSeqBAIJSetPreallocation(fact, bs, MAT_SKIP_ALLOCATION, NULL);

480:   b               = (Mat_SeqBAIJ *)(fact)->data;
481:   b->free_a       = PETSC_TRUE;
482:   b->free_ij      = PETSC_TRUE;
483:   b->singlemalloc = PETSC_FALSE;

485:   PetscMalloc1(bs2 * (bdiag[0] + 1), &b->a);

487:   b->j         = bj;
488:   b->i         = bi;
489:   b->diag      = bdiag;
490:   b->free_diag = PETSC_TRUE;
491:   b->ilen      = NULL;
492:   b->imax      = NULL;
493:   b->row       = isrow;
494:   b->col       = iscol;
495:   PetscObjectReference((PetscObject)isrow);
496:   PetscObjectReference((PetscObject)iscol);
497:   b->icol = isicol;

499:   PetscMalloc1(bs * n + bs, &b->solve_work);
500:   /* In b structure:  Free imax, ilen, old a, old j.
501:      Allocate bdiag, solve_work, new a, new j */
502:   b->maxnz = b->nz = bdiag[0] + 1;

504:   fact->info.factor_mallocs    = reallocs;
505:   fact->info.fill_ratio_given  = f;
506:   fact->info.fill_ratio_needed = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);

508:   MatSeqBAIJSetNumericFactorization(fact, both_identity);
509:   return 0;
510: }

512: /*
513:      This code is virtually identical to MatILUFactorSymbolic_SeqAIJ
514:    except that the data structure of Mat_SeqAIJ is slightly different.
515:    Not a good example of code reuse.
516: */
517: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_inplace(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
518: {
519:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ *)A->data, *b;
520:   IS              isicol;
521:   const PetscInt *r, *ic, *ai = a->i, *aj = a->j, *xi;
522:   PetscInt        prow, n = a->mbs, *ainew, *ajnew, jmax, *fill, nz, *im, *ajfill, *flev, *xitmp;
523:   PetscInt       *dloc, idx, row, m, fm, nzf, nzi, reallocate = 0, dcount = 0;
524:   PetscInt        incrlev, nnz, i, bs = A->rmap->bs, bs2 = a->bs2, levels, diagonal_fill, dd;
525:   PetscBool       col_identity, row_identity, both_identity, flg;
526:   PetscReal       f;

528:   MatMissingDiagonal_SeqBAIJ(A, &flg, &dd);

531:   f             = info->fill;
532:   levels        = (PetscInt)info->levels;
533:   diagonal_fill = (PetscInt)info->diagonal_fill;

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

537:   ISIdentity(isrow, &row_identity);
538:   ISIdentity(iscol, &col_identity);
539:   both_identity = (PetscBool)(row_identity && col_identity);

541:   if (!levels && both_identity) { /* special case copy the nonzero structure */
542:     MatDuplicateNoCreate_SeqBAIJ(fact, A, MAT_DO_NOT_COPY_VALUES, PETSC_TRUE);
543:     MatSeqBAIJSetNumericFactorization_inplace(fact, both_identity);

545:     fact->factortype = MAT_FACTOR_ILU;
546:     b                = (Mat_SeqBAIJ *)fact->data;
547:     b->row           = isrow;
548:     b->col           = iscol;
549:     PetscObjectReference((PetscObject)isrow);
550:     PetscObjectReference((PetscObject)iscol);
551:     b->icol          = isicol;
552:     b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;

554:     PetscMalloc1((n + 1) * bs, &b->solve_work);
555:     return 0;
556:   }

558:   /* general case perform the symbolic factorization */
559:   ISGetIndices(isrow, &r);
560:   ISGetIndices(isicol, &ic);

562:   /* get new row pointers */
563:   PetscMalloc1(n + 1, &ainew);
564:   ainew[0] = 0;
565:   /* don't know how many column pointers are needed so estimate */
566:   jmax = (PetscInt)(f * ai[n] + 1);
567:   PetscMalloc1(jmax, &ajnew);
568:   /* ajfill is level of fill for each fill entry */
569:   PetscMalloc1(jmax, &ajfill);
570:   /* fill is a linked list of nonzeros in active row */
571:   PetscMalloc1(n + 1, &fill);
572:   /* im is level for each filled value */
573:   PetscMalloc1(n + 1, &im);
574:   /* dloc is location of diagonal in factor */
575:   PetscMalloc1(n + 1, &dloc);
576:   dloc[0] = 0;
577:   for (prow = 0; prow < n; prow++) {
578:     /* copy prow into linked list */
579:     nzf = nz = ai[r[prow] + 1] - ai[r[prow]];
581:     xi         = aj + ai[r[prow]];
582:     fill[n]    = n;
583:     fill[prow] = -1; /* marker for diagonal entry */
584:     while (nz--) {
585:       fm  = n;
586:       idx = ic[*xi++];
587:       do {
588:         m  = fm;
589:         fm = fill[m];
590:       } while (fm < idx);
591:       fill[m]   = idx;
592:       fill[idx] = fm;
593:       im[idx]   = 0;
594:     }

596:     /* make sure diagonal entry is included */
597:     if (diagonal_fill && fill[prow] == -1) {
598:       fm = n;
599:       while (fill[fm] < prow) fm = fill[fm];
600:       fill[prow] = fill[fm]; /* insert diagonal into linked list */
601:       fill[fm]   = prow;
602:       im[prow]   = 0;
603:       nzf++;
604:       dcount++;
605:     }

607:     nzi = 0;
608:     row = fill[n];
609:     while (row < prow) {
610:       incrlev = im[row] + 1;
611:       nz      = dloc[row];
612:       xi      = ajnew + ainew[row] + nz + 1;
613:       flev    = ajfill + ainew[row] + nz + 1;
614:       nnz     = ainew[row + 1] - ainew[row] - nz - 1;
615:       fm      = row;
616:       while (nnz-- > 0) {
617:         idx = *xi++;
618:         if (*flev + incrlev > levels) {
619:           flev++;
620:           continue;
621:         }
622:         do {
623:           m  = fm;
624:           fm = fill[m];
625:         } while (fm < idx);
626:         if (fm != idx) {
627:           im[idx]   = *flev + incrlev;
628:           fill[m]   = idx;
629:           fill[idx] = fm;
630:           fm        = idx;
631:           nzf++;
632:         } else if (im[idx] > *flev + incrlev) im[idx] = *flev + incrlev;
633:         flev++;
634:       }
635:       row = fill[row];
636:       nzi++;
637:     }
638:     /* copy new filled row into permanent storage */
639:     ainew[prow + 1] = ainew[prow] + nzf;
640:     if (ainew[prow + 1] > jmax) {
641:       /* estimate how much additional space we will need */
642:       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
643:       /* just double the memory each time */
644:       PetscInt maxadd = jmax;
645:       /* maxadd = (int)(((f*ai[n]+1)*(n-prow+5))/n); */
646:       if (maxadd < nzf) maxadd = (n - prow) * (nzf + 1);
647:       jmax += maxadd;

649:       /* allocate a longer ajnew and ajfill */
650:       PetscMalloc1(jmax, &xitmp);
651:       PetscArraycpy(xitmp, ajnew, ainew[prow]);
652:       PetscFree(ajnew);
653:       ajnew = xitmp;
654:       PetscMalloc1(jmax, &xitmp);
655:       PetscArraycpy(xitmp, ajfill, ainew[prow]);
656:       PetscFree(ajfill);
657:       ajfill = xitmp;
658:       reallocate++; /* count how many reallocations are needed */
659:     }
660:     xitmp      = ajnew + ainew[prow];
661:     flev       = ajfill + ainew[prow];
662:     dloc[prow] = nzi;
663:     fm         = fill[n];
664:     while (nzf--) {
665:       *xitmp++ = fm;
666:       *flev++  = im[fm];
667:       fm       = fill[fm];
668:     }
669:     /* make sure row has diagonal entry */
671:                                                         try running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill",
672:                prow);
673:   }
674:   PetscFree(ajfill);
675:   ISRestoreIndices(isrow, &r);
676:   ISRestoreIndices(isicol, &ic);
677:   PetscFree(fill);
678:   PetscFree(im);

680: #if defined(PETSC_USE_INFO)
681:   {
682:     PetscReal af = ((PetscReal)ainew[n]) / ((PetscReal)ai[n]);
683:     PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocate, (double)f, (double)af);
684:     PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af);
685:     PetscInfo(A, "PCFactorSetFill(pc,%g);\n", (double)af);
686:     PetscInfo(A, "for best performance.\n");
687:     if (diagonal_fill) PetscInfo(A, "Detected and replaced %" PetscInt_FMT " missing diagonals\n", dcount);
688:   }
689: #endif

691:   /* put together the new matrix */
692:   MatSeqBAIJSetPreallocation(fact, bs, MAT_SKIP_ALLOCATION, NULL);
693:   b = (Mat_SeqBAIJ *)fact->data;

695:   b->free_a       = PETSC_TRUE;
696:   b->free_ij      = PETSC_TRUE;
697:   b->singlemalloc = PETSC_FALSE;

699:   PetscMalloc1(bs2 * ainew[n], &b->a);

701:   b->j = ajnew;
702:   b->i = ainew;
703:   for (i = 0; i < n; i++) dloc[i] += ainew[i];
704:   b->diag          = dloc;
705:   b->free_diag     = PETSC_TRUE;
706:   b->ilen          = NULL;
707:   b->imax          = NULL;
708:   b->row           = isrow;
709:   b->col           = iscol;
710:   b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;

712:   PetscObjectReference((PetscObject)isrow);
713:   PetscObjectReference((PetscObject)iscol);
714:   b->icol = isicol;
715:   PetscMalloc1(bs * n + bs, &b->solve_work);
716:   /* In b structure:  Free imax, ilen, old a, old j.
717:      Allocate dloc, solve_work, new a, new j */
718:   b->maxnz = b->nz = ainew[n];

720:   fact->info.factor_mallocs    = reallocate;
721:   fact->info.fill_ratio_given  = f;
722:   fact->info.fill_ratio_needed = ((PetscReal)ainew[n]) / ((PetscReal)ai[prow]);

724:   MatSeqBAIJSetNumericFactorization_inplace(fact, both_identity);
725:   return 0;
726: }

728: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE(Mat A)
729: {
730:   /* Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; */
731:   /* int i,*AJ=a->j,nz=a->nz; */

733:   /* Undo Column scaling */
734:   /*    while (nz--) { */
735:   /*      AJ[i] = AJ[i]/4; */
736:   /*    } */
737:   /* This should really invoke a push/pop logic, but we don't have that yet. */
738:   A->ops->setunfactored = NULL;
739:   return 0;
740: }

742: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE_usj(Mat A)
743: {
744:   Mat_SeqBAIJ    *a  = (Mat_SeqBAIJ *)A->data;
745:   PetscInt       *AJ = a->j, nz = a->nz;
746:   unsigned short *aj = (unsigned short *)AJ;

748:   /* Is this really necessary? */
749:   while (nz--) { AJ[nz] = (int)((unsigned int)aj[nz]); /* First extend, then convert to signed. */ }
750:   A->ops->setunfactored = NULL;
751:   return 0;
752: }