Actual source code: baijfact2.c
1: /*
2: Factorization code for BAIJ format.
3: */
5: #include <../src/mat/impls/baij/seq/baij.h>
6: #include <petsc/private/kernels/blockinvert.h>
7: #include <../src/mat/utils/freespace.h>
9: PETSC_INTERN PetscErrorCode MatDuplicateNoCreate_SeqBAIJ(Mat, Mat, MatDuplicateOption, PetscBool);
11: /*
12: This is not much faster than MatLUFactorNumeric_SeqBAIJ_N() but the solve is faster at least sometimes
13: */
14: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_15_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info)
15: {
16: Mat C = B;
17: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
18: PetscInt i, j, k, ipvt[15];
19: const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *ajtmp, *bjtmp, *bdiag = b->diag, *pj;
20: PetscInt nz, nzL, row;
21: MatScalar *rtmp, *pc, *mwork, *pv, *vv, work[225];
22: const MatScalar *v, *aa = a->a;
23: PetscInt bs2 = a->bs2, bs = A->rmap->bs, flg;
24: PetscInt sol_ver;
25: PetscBool allowzeropivot, zeropivotdetected;
27: PetscFunctionBegin;
28: allowzeropivot = PetscNot(A->erroriffailure);
29: PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)A)->prefix, "-sol_ver", &sol_ver, NULL));
31: /* generate work space needed by the factorization */
32: PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
33: PetscCall(PetscArrayzero(rtmp, bs2 * n));
35: for (i = 0; i < n; i++) {
36: /* zero rtmp */
37: /* L part */
38: nz = bi[i + 1] - bi[i];
39: bjtmp = bj + bi[i];
40: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
42: /* U part */
43: nz = bdiag[i] - bdiag[i + 1];
44: bjtmp = bj + bdiag[i + 1] + 1;
45: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
47: /* load in initial (unfactored row) */
48: nz = ai[i + 1] - ai[i];
49: ajtmp = aj + ai[i];
50: v = aa + bs2 * ai[i];
51: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ajtmp[j], v + bs2 * j, bs2));
53: /* elimination */
54: bjtmp = bj + bi[i];
55: nzL = bi[i + 1] - bi[i];
56: for (k = 0; k < nzL; k++) {
57: row = bjtmp[k];
58: pc = rtmp + bs2 * row;
59: for (flg = 0, j = 0; j < bs2; j++) {
60: if (pc[j] != 0.0) {
61: flg = 1;
62: break;
63: }
64: }
65: if (flg) {
66: pv = b->a + bs2 * bdiag[row];
67: PetscKernel_A_gets_A_times_B(bs, pc, pv, mwork);
68: /* PetscCall(PetscKernel_A_gets_A_times_B_15(pc,pv,mwork)); */
69: pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
70: pv = b->a + bs2 * (bdiag[row + 1] + 1);
71: nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
72: for (j = 0; j < nz; j++) {
73: vv = rtmp + bs2 * pj[j];
74: PetscKernel_A_gets_A_minus_B_times_C(bs, vv, pc, pv);
75: /* PetscCall(PetscKernel_A_gets_A_minus_B_times_C_15(vv,pc,pv)); */
76: pv += bs2;
77: }
78: PetscCall(PetscLogFlops(2.0 * bs2 * bs * (nz + 1) - bs2)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
79: }
80: }
82: /* finished row so stick it into b->a */
83: /* L part */
84: pv = b->a + bs2 * bi[i];
85: pj = b->j + bi[i];
86: nz = bi[i + 1] - bi[i];
87: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
89: /* Mark diagonal and invert diagonal for simpler triangular solves */
90: pv = b->a + bs2 * bdiag[i];
91: pj = b->j + bdiag[i];
92: PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
93: PetscCall(PetscKernel_A_gets_inverse_A_15(pv, ipvt, work, info->shiftamount, allowzeropivot, &zeropivotdetected));
94: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
96: /* U part */
97: pv = b->a + bs2 * (bdiag[i + 1] + 1);
98: pj = b->j + bdiag[i + 1] + 1;
99: nz = bdiag[i] - bdiag[i + 1] - 1;
100: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
101: }
103: PetscCall(PetscFree2(rtmp, mwork));
105: C->ops->solve = MatSolve_SeqBAIJ_15_NaturalOrdering_ver1;
106: C->ops->solvetranspose = MatSolve_SeqBAIJ_N_NaturalOrdering;
107: C->assembled = PETSC_TRUE;
109: PetscCall(PetscLogFlops(1.333333333333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_N(Mat B, Mat A, const MatFactorInfo *info)
114: {
115: Mat C = B;
116: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
117: IS isrow = b->row, isicol = b->icol;
118: const PetscInt *r, *ic;
119: PetscInt i, j, k, n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
120: PetscInt *ajtmp, *bjtmp, nz, nzL, row, *bdiag = b->diag, *pj;
121: MatScalar *rtmp, *pc, *mwork, *v, *pv, *aa = a->a;
122: PetscInt bs = A->rmap->bs, bs2 = a->bs2, *v_pivots, flg;
123: MatScalar *v_work;
124: PetscBool col_identity, row_identity, both_identity;
125: PetscBool allowzeropivot, zeropivotdetected;
127: PetscFunctionBegin;
128: PetscCall(ISGetIndices(isrow, &r));
129: PetscCall(ISGetIndices(isicol, &ic));
130: allowzeropivot = PetscNot(A->erroriffailure);
132: PetscCall(PetscCalloc1(bs2 * n, &rtmp));
134: /* generate work space needed by dense LU factorization */
135: PetscCall(PetscMalloc3(bs, &v_work, bs2, &mwork, bs, &v_pivots));
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[r[i] + 1] - ai[r[i]];
151: ajtmp = aj + ai[r[i]];
152: v = aa + bs2 * ai[r[i]];
153: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ic[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] != 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: pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
171: pv = b->a + bs2 * (bdiag[row + 1] + 1);
172: nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
173: for (j = 0; j < nz; j++) PetscKernel_A_gets_A_minus_B_times_C(bs, rtmp + bs2 * pj[j], pc, pv + bs2 * j);
174: PetscCall(PetscLogFlops(2.0 * bs2 * bs * (nz + 1) - bs2)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
175: }
176: }
178: /* finished row so stick it into b->a */
179: /* L part */
180: pv = b->a + bs2 * bi[i];
181: pj = b->j + bi[i];
182: nz = bi[i + 1] - bi[i];
183: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
185: /* Mark diagonal and invert diagonal for simpler triangular solves */
186: pv = b->a + bs2 * bdiag[i];
187: pj = b->j + bdiag[i];
188: PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
190: PetscCall(PetscKernel_A_gets_inverse_A(bs, pv, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
191: if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
193: /* U part */
194: pv = b->a + bs2 * (bdiag[i + 1] + 1);
195: pj = b->j + bdiag[i + 1] + 1;
196: nz = bdiag[i] - bdiag[i + 1] - 1;
197: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
198: }
200: PetscCall(PetscFree(rtmp));
201: PetscCall(PetscFree3(v_work, mwork, v_pivots));
202: PetscCall(ISRestoreIndices(isicol, &ic));
203: PetscCall(ISRestoreIndices(isrow, &r));
205: PetscCall(ISIdentity(isrow, &row_identity));
206: PetscCall(ISIdentity(isicol, &col_identity));
208: both_identity = (PetscBool)(row_identity && col_identity);
209: if (both_identity) {
210: switch (bs) {
211: case 9:
212: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
213: C->ops->solve = MatSolve_SeqBAIJ_9_NaturalOrdering;
214: #else
215: C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
216: #endif
217: break;
218: case 11:
219: C->ops->solve = MatSolve_SeqBAIJ_11_NaturalOrdering;
220: break;
221: case 12:
222: C->ops->solve = MatSolve_SeqBAIJ_12_NaturalOrdering;
223: break;
224: case 13:
225: C->ops->solve = MatSolve_SeqBAIJ_13_NaturalOrdering;
226: break;
227: case 14:
228: C->ops->solve = MatSolve_SeqBAIJ_14_NaturalOrdering;
229: break;
230: default:
231: C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
232: break;
233: }
234: } else {
235: C->ops->solve = MatSolve_SeqBAIJ_N;
236: }
237: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N;
239: C->assembled = PETSC_TRUE;
241: PetscCall(PetscLogFlops(1.333333333333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /*
246: ilu(0) with natural ordering under new data structure.
247: See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed description
248: because this code is almost identical to MatILUFactorSymbolic_SeqAIJ_ilu0_inplace().
249: */
251: static PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_ilu0(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
252: {
253: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b;
254: const PetscInt n = a->mbs, *ai = a->i, *aj, *adiag, bs2 = a->bs2;
255: PetscInt i, j, nz, *bi, *bj, *bdiag, bi_temp;
257: PetscFunctionBegin;
258: PetscCall(MatGetDiagonalMarkers_SeqBAIJ(A, &adiag, NULL));
259: PetscCall(MatDuplicateNoCreate_SeqBAIJ(fact, A, MAT_DO_NOT_COPY_VALUES, PETSC_FALSE));
260: b = (Mat_SeqBAIJ *)fact->data;
262: /* allocate matrix arrays for new data structure */
263: PetscCall(PetscShmgetAllocateArray(bs2 * ai[n], sizeof(PetscScalar), (void **)&b->a));
264: PetscCall(PetscShmgetAllocateArray(ai[n], sizeof(PetscInt), (void **)&b->j));
265: PetscCall(PetscShmgetAllocateArray(n + 1, sizeof(PetscInt), (void **)&b->i));
266: b->free_a = PETSC_TRUE;
267: b->free_ij = PETSC_TRUE;
268: fact->preallocated = PETSC_TRUE;
269: fact->assembled = PETSC_TRUE;
270: if (!b->diag) PetscCall(PetscMalloc1(n + 1, &b->diag));
271: bdiag = b->diag;
273: if (n > 0) PetscCall(PetscArrayzero(b->a, bs2 * ai[n]));
275: /* set bi and bj with new data structure */
276: bi = b->i;
277: bj = b->j;
279: /* L part */
280: bi[0] = 0;
281: for (i = 0; i < n; i++) {
282: nz = adiag[i] - ai[i];
283: bi[i + 1] = bi[i] + nz;
284: aj = a->j + ai[i];
285: for (j = 0; j < nz; j++) {
286: *bj = aj[j];
287: bj++;
288: }
289: }
291: /* U part */
292: bi_temp = bi[n];
293: bdiag[n] = bi[n] - 1;
294: for (i = n - 1; i >= 0; i--) {
295: nz = ai[i + 1] - adiag[i] - 1;
296: bi_temp = bi_temp + nz + 1;
297: aj = a->j + adiag[i] + 1;
298: for (j = 0; j < nz; j++) {
299: *bj = aj[j];
300: bj++;
301: }
302: /* diag[i] */
303: *bj = i;
304: bj++;
305: bdiag[i] = bi_temp - 1;
306: }
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
311: {
312: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b;
313: IS isicol;
314: const PetscInt *r, *ic;
315: PetscInt n = a->mbs, *ai = a->i, *aj = a->j;
316: PetscInt *bi, *cols, nnz, *cols_lvl;
317: PetscInt *bdiag, prow, fm, nzbd, reallocs = 0, dcount = 0;
318: PetscInt i, levels, diagonal_fill;
319: PetscBool col_identity, row_identity, both_identity;
320: PetscReal f;
321: PetscInt nlnk, *lnk, *lnk_lvl = NULL;
322: PetscBT lnkbt;
323: PetscInt nzi, *bj, **bj_ptr, **bjlvl_ptr;
324: PetscFreeSpaceList free_space = NULL, current_space = NULL;
325: PetscFreeSpaceList free_space_lvl = NULL, current_space_lvl = NULL;
326: PetscInt bs = A->rmap->bs, bs2 = a->bs2;
327: PetscBool diagDense;
329: PetscFunctionBegin;
330: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must be square matrix, rows %" PetscInt_FMT " columns %" PetscInt_FMT, A->rmap->n, A->cmap->n);
331: if (bs > 1) { /* check shifttype */
332: PetscCheck(info->shifttype != (PetscReal)MAT_SHIFT_NONZERO && info->shifttype != (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only MAT_SHIFT_NONE and MAT_SHIFT_INBLOCKS are supported for BAIJ matrix");
333: }
334: PetscCall(MatGetDiagonalMarkers_SeqBAIJ(A, NULL, &diagDense));
335: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry");
337: f = info->fill;
338: levels = (PetscInt)info->levels;
339: diagonal_fill = (PetscInt)info->diagonal_fill;
341: PetscCall(ISInvertPermutation(iscol, PETSC_DECIDE, &isicol));
343: PetscCall(ISIdentity(isrow, &row_identity));
344: PetscCall(ISIdentity(iscol, &col_identity));
346: both_identity = (PetscBool)(row_identity && col_identity);
348: if (!levels && both_identity) {
349: /* special case: ilu(0) with natural ordering */
350: PetscCall(MatILUFactorSymbolic_SeqBAIJ_ilu0(fact, A, isrow, iscol, info));
351: PetscCall(MatSeqBAIJSetNumericFactorization(fact, both_identity));
353: fact->factortype = MAT_FACTOR_ILU;
354: fact->info.factor_mallocs = 0;
355: fact->info.fill_ratio_given = info->fill;
356: fact->info.fill_ratio_needed = 1.0;
358: b = (Mat_SeqBAIJ *)fact->data;
359: b->row = isrow;
360: b->col = iscol;
361: b->icol = isicol;
362: PetscCall(PetscObjectReference((PetscObject)isrow));
363: PetscCall(PetscObjectReference((PetscObject)iscol));
364: b->pivotinblocks = info->pivotinblocks ? PETSC_TRUE : PETSC_FALSE;
366: PetscCall(PetscMalloc1((n + 1) * bs, &b->solve_work));
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: PetscCall(ISGetIndices(isrow, &r));
371: PetscCall(ISGetIndices(isicol, &ic));
373: /* get new row pointers */
374: PetscCall(PetscMalloc1(n + 1, &bi));
375: bi[0] = 0;
376: /* bdiag is location of diagonal in factor */
377: PetscCall(PetscMalloc1(n + 1, &bdiag));
378: bdiag[0] = 0;
380: PetscCall(PetscMalloc2(n, &bj_ptr, n, &bjlvl_ptr));
382: /* create a linked list for storing column indices of the active row */
383: nlnk = n + 1;
384: PetscCall(PetscIncompleteLLCreate(n, n, nlnk, lnk, lnk_lvl, lnkbt));
386: /* initial FreeSpace size is f*(ai[n]+1) */
387: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space));
388: current_space = free_space;
389: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space_lvl));
390: current_space_lvl = free_space_lvl;
392: for (i = 0; i < n; i++) {
393: nzi = 0;
394: /* copy current row into linked list */
395: nnz = ai[r[i] + 1] - ai[r[i]];
396: PetscCheck(nnz, 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);
397: cols = aj + ai[r[i]];
398: lnk[i] = -1; /* marker to indicate if diagonal exists */
399: PetscCall(PetscIncompleteLLInit(nnz, cols, n, ic, &nlnk, lnk, lnk_lvl, lnkbt));
400: nzi += nlnk;
402: /* make sure diagonal entry is included */
403: if (diagonal_fill && lnk[i] == -1) {
404: fm = n;
405: while (lnk[fm] < i) fm = lnk[fm];
406: lnk[i] = lnk[fm]; /* insert diagonal into linked list */
407: lnk[fm] = i;
408: lnk_lvl[i] = 0;
409: nzi++;
410: dcount++;
411: }
413: /* add pivot rows into the active row */
414: nzbd = 0;
415: prow = lnk[n];
416: while (prow < i) {
417: nnz = bdiag[prow];
418: cols = bj_ptr[prow] + nnz + 1;
419: cols_lvl = bjlvl_ptr[prow] + nnz + 1;
420: nnz = bi[prow + 1] - bi[prow] - nnz - 1;
422: PetscCall(PetscILULLAddSorted(nnz, cols, levels, cols_lvl, prow, &nlnk, lnk, lnk_lvl, lnkbt, prow));
423: nzi += nlnk;
424: prow = lnk[prow];
425: nzbd++;
426: }
427: bdiag[i] = nzbd;
428: bi[i + 1] = bi[i] + nzi;
430: /* if free space is not available, make more free space */
431: if (current_space->local_remaining < nzi) {
432: nnz = PetscIntMultTruncate(2, PetscIntMultTruncate(nzi, n - i)); /* estimated and max additional space needed */
433: PetscCall(PetscFreeSpaceGet(nnz, ¤t_space));
434: PetscCall(PetscFreeSpaceGet(nnz, ¤t_space_lvl));
435: reallocs++;
436: }
438: /* copy data into free_space and free_space_lvl, then initialize lnk */
439: PetscCall(PetscIncompleteLLClean(n, n, nzi, lnk, lnk_lvl, current_space->array, current_space_lvl->array, lnkbt));
441: bj_ptr[i] = current_space->array;
442: bjlvl_ptr[i] = current_space_lvl->array;
444: /* make sure the active row i has diagonal entry */
445: PetscCheck(*(bj_ptr[i] + bdiag[i]) == i, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Row %" PetscInt_FMT " has missing diagonal in factored matrix, try running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill", i);
447: current_space->array += nzi;
448: current_space->local_used += nzi;
449: current_space->local_remaining -= nzi;
451: current_space_lvl->array += nzi;
452: current_space_lvl->local_used += nzi;
453: current_space_lvl->local_remaining -= nzi;
454: }
456: PetscCall(ISRestoreIndices(isrow, &r));
457: PetscCall(ISRestoreIndices(isicol, &ic));
459: /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
460: PetscCall(PetscMalloc1(bi[n], &bj));
461: PetscCall(PetscFreeSpaceContiguous_LU(&free_space, bj, n, bi, bdiag));
463: PetscCall(PetscIncompleteLLDestroy(lnk, lnkbt));
464: PetscCall(PetscFreeSpaceDestroy(free_space_lvl));
465: PetscCall(PetscFree2(bj_ptr, bjlvl_ptr));
467: #if PetscDefined(USE_INFO)
468: {
469: PetscReal af = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);
470: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)f, (double)af));
471: PetscCall(PetscInfo(A, "Run with -[sub_]pc_factor_fill %g or use \n", (double)af));
472: PetscCall(PetscInfo(A, "PCFactorSetFill([sub]pc,%g);\n", (double)af));
473: PetscCall(PetscInfo(A, "for best performance.\n"));
474: if (diagonal_fill) PetscCall(PetscInfo(A, "Detected and replaced %" PetscInt_FMT " missing diagonals\n", dcount));
475: }
476: #endif
478: /* put together the new matrix */
479: PetscCall(MatSeqBAIJSetPreallocation(fact, bs, MAT_SKIP_ALLOCATION, NULL));
481: b = (Mat_SeqBAIJ *)fact->data;
482: b->free_ij = PETSC_TRUE;
483: PetscCall(PetscShmgetAllocateArray(bs2 * (bdiag[0] + 1), sizeof(PetscScalar), (void **)&b->a));
484: b->free_a = PETSC_TRUE;
486: b->j = bj;
487: b->i = bi;
488: b->diag = bdiag;
489: b->ilen = NULL;
490: b->imax = NULL;
491: b->row = isrow;
492: b->col = iscol;
493: PetscCall(PetscObjectReference((PetscObject)isrow));
494: PetscCall(PetscObjectReference((PetscObject)iscol));
495: b->icol = isicol;
497: PetscCall(PetscMalloc1(bs * n + bs, &b->solve_work));
498: /* In b structure: Free imax, ilen, old a, old j.
499: Allocate bdiag, solve_work, new a, new j */
500: b->maxnz = b->nz = bdiag[0] + 1;
502: fact->info.factor_mallocs = reallocs;
503: fact->info.fill_ratio_given = f;
504: fact->info.fill_ratio_needed = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);
506: PetscCall(MatSeqBAIJSetNumericFactorization(fact, both_identity));
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }