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 <petscbt.h>
8: #include <../src/mat/utils/freespace.h>
10: extern PetscErrorCode MatDuplicateNoCreate_SeqBAIJ(Mat, Mat, MatDuplicateOption, PetscBool);
12: /*
13: This is not much faster than MatLUFactorNumeric_SeqBAIJ_N() but the solve is faster at least sometimes
14: */
15: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_15_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info)
16: {
17: Mat C = B;
18: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
19: PetscInt i, j, k, ipvt[15];
20: const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *ajtmp, *bjtmp, *bdiag = b->diag, *pj;
21: PetscInt nz, nzL, row;
22: MatScalar *rtmp, *pc, *mwork, *pv, *vv, work[225];
23: const MatScalar *v, *aa = a->a;
24: PetscInt bs2 = a->bs2, bs = A->rmap->bs, flg;
25: PetscInt sol_ver;
26: PetscBool allowzeropivot, zeropivotdetected;
28: PetscFunctionBegin;
29: allowzeropivot = PetscNot(A->erroriffailure);
30: PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)A)->prefix, "-sol_ver", &sol_ver, NULL));
32: /* generate work space needed by the factorization */
33: PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
34: PetscCall(PetscArrayzero(rtmp, bs2 * n));
36: for (i = 0; i < n; i++) {
37: /* zero rtmp */
38: /* L part */
39: nz = bi[i + 1] - bi[i];
40: bjtmp = bj + bi[i];
41: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
43: /* U part */
44: nz = bdiag[i] - bdiag[i + 1];
45: bjtmp = bj + bdiag[i + 1] + 1;
46: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
48: /* load in initial (unfactored row) */
49: nz = ai[i + 1] - ai[i];
50: ajtmp = aj + ai[i];
51: v = aa + bs2 * ai[i];
52: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ajtmp[j], v + bs2 * j, bs2));
54: /* elimination */
55: bjtmp = bj + bi[i];
56: nzL = bi[i + 1] - bi[i];
57: for (k = 0; k < nzL; k++) {
58: row = bjtmp[k];
59: pc = rtmp + bs2 * row;
60: for (flg = 0, j = 0; j < bs2; j++) {
61: if (pc[j] != 0.0) {
62: flg = 1;
63: break;
64: }
65: }
66: if (flg) {
67: pv = b->a + bs2 * bdiag[row];
68: PetscKernel_A_gets_A_times_B(bs, pc, pv, mwork);
69: /* PetscCall(PetscKernel_A_gets_A_times_B_15(pc,pv,mwork)); */
70: pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
71: pv = b->a + bs2 * (bdiag[row + 1] + 1);
72: nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
73: for (j = 0; j < nz; j++) {
74: vv = rtmp + bs2 * pj[j];
75: PetscKernel_A_gets_A_minus_B_times_C(bs, vv, pc, pv);
76: /* PetscCall(PetscKernel_A_gets_A_minus_B_times_C_15(vv,pc,pv)); */
77: pv += bs2;
78: }
79: PetscCall(PetscLogFlops(2.0 * bs2 * bs * (nz + 1) - bs2)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
80: }
81: }
83: /* finished row so stick it into b->a */
84: /* L part */
85: pv = b->a + bs2 * bi[i];
86: pj = b->j + bi[i];
87: nz = bi[i + 1] - bi[i];
88: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
90: /* Mark diagonal and invert diagonal for simpler triangular solves */
91: pv = b->a + bs2 * bdiag[i];
92: pj = b->j + bdiag[i];
93: PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
94: PetscCall(PetscKernel_A_gets_inverse_A_15(pv, ipvt, work, info->shiftamount, allowzeropivot, &zeropivotdetected));
95: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
97: /* U part */
98: pv = b->a + bs2 * (bdiag[i + 1] + 1);
99: pj = b->j + bdiag[i + 1] + 1;
100: nz = bdiag[i] - bdiag[i + 1] - 1;
101: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
102: }
104: PetscCall(PetscFree2(rtmp, mwork));
106: C->ops->solve = MatSolve_SeqBAIJ_15_NaturalOrdering_ver1;
107: C->ops->solvetranspose = MatSolve_SeqBAIJ_N_NaturalOrdering;
108: C->assembled = PETSC_TRUE;
110: PetscCall(PetscLogFlops(1.333333333333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_N(Mat B, Mat A, const MatFactorInfo *info)
115: {
116: Mat C = B;
117: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
118: IS isrow = b->row, isicol = b->icol;
119: const PetscInt *r, *ic;
120: PetscInt i, j, k, n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
121: PetscInt *ajtmp, *bjtmp, nz, nzL, row, *bdiag = b->diag, *pj;
122: MatScalar *rtmp, *pc, *mwork, *v, *pv, *aa = a->a;
123: PetscInt bs = A->rmap->bs, bs2 = a->bs2, *v_pivots, flg;
124: MatScalar *v_work;
125: PetscBool col_identity, row_identity, both_identity;
126: PetscBool allowzeropivot, zeropivotdetected;
128: PetscFunctionBegin;
129: PetscCall(ISGetIndices(isrow, &r));
130: PetscCall(ISGetIndices(isicol, &ic));
131: allowzeropivot = PetscNot(A->erroriffailure);
133: PetscCall(PetscCalloc1(bs2 * n, &rtmp));
135: /* generate work space needed by dense LU factorization */
136: PetscCall(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++) PetscCall(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++) PetscCall(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++) PetscCall(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: PetscCall(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++) PetscCall(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: PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
191: PetscCall(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++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
199: }
201: PetscCall(PetscFree(rtmp));
202: PetscCall(PetscFree3(v_work, mwork, v_pivots));
203: PetscCall(ISRestoreIndices(isicol, &ic));
204: PetscCall(ISRestoreIndices(isrow, &r));
206: PetscCall(ISIdentity(isrow, &row_identity));
207: PetscCall(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: PetscCall(PetscLogFlops(1.333333333333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */
243: PetscFunctionReturn(PETSC_SUCCESS);
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: static 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: PetscFunctionBegin;
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] + 1, sizeof(PetscScalar), (void **)&b->a));
264: PetscCall(PetscShmgetAllocateArray(ai[n] + 1, 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, d;
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: PetscBool missing;
327: PetscInt bs = A->rmap->bs, bs2 = a->bs2;
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: }
335: PetscCall(MatMissingDiagonal(A, &missing, &d));
336: PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry %" PetscInt_FMT, d);
338: f = info->fill;
339: levels = (PetscInt)info->levels;
340: diagonal_fill = (PetscInt)info->diagonal_fill;
342: PetscCall(ISInvertPermutation(iscol, PETSC_DECIDE, &isicol));
344: PetscCall(ISIdentity(isrow, &row_identity));
345: PetscCall(ISIdentity(iscol, &col_identity));
347: both_identity = (PetscBool)(row_identity && col_identity);
349: if (!levels && both_identity) {
350: /* special case: ilu(0) with natural ordering */
351: PetscCall(MatILUFactorSymbolic_SeqBAIJ_ilu0(fact, A, isrow, iscol, info));
352: PetscCall(MatSeqBAIJSetNumericFactorization(fact, both_identity));
354: fact->factortype = MAT_FACTOR_ILU;
355: fact->info.factor_mallocs = 0;
356: fact->info.fill_ratio_given = info->fill;
357: fact->info.fill_ratio_needed = 1.0;
359: b = (Mat_SeqBAIJ *)fact->data;
360: b->row = isrow;
361: b->col = iscol;
362: b->icol = isicol;
363: PetscCall(PetscObjectReference((PetscObject)isrow));
364: PetscCall(PetscObjectReference((PetscObject)iscol));
365: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
367: PetscCall(PetscMalloc1((n + 1) * bs, &b->solve_work));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: PetscCall(ISGetIndices(isrow, &r));
372: PetscCall(ISGetIndices(isicol, &ic));
374: /* get new row pointers */
375: PetscCall(PetscMalloc1(n + 1, &bi));
376: bi[0] = 0;
377: /* bdiag is location of diagonal in factor */
378: PetscCall(PetscMalloc1(n + 1, &bdiag));
379: bdiag[0] = 0;
381: PetscCall(PetscMalloc2(n, &bj_ptr, n, &bjlvl_ptr));
383: /* create a linked list for storing column indices of the active row */
384: nlnk = n + 1;
385: PetscCall(PetscIncompleteLLCreate(n, n, nlnk, lnk, lnk_lvl, lnkbt));
387: /* initial FreeSpace size is f*(ai[n]+1) */
388: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space));
389: current_space = free_space;
390: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space_lvl));
391: current_space_lvl = free_space_lvl;
393: for (i = 0; i < n; i++) {
394: nzi = 0;
395: /* copy current row into linked list */
396: nnz = ai[r[i] + 1] - ai[r[i]];
397: 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);
398: cols = aj + ai[r[i]];
399: lnk[i] = -1; /* marker to indicate if diagonal exists */
400: PetscCall(PetscIncompleteLLInit(nnz, cols, n, ic, &nlnk, lnk, lnk_lvl, lnkbt));
401: nzi += nlnk;
403: /* make sure diagonal entry is included */
404: if (diagonal_fill && lnk[i] == -1) {
405: fm = n;
406: while (lnk[fm] < i) fm = lnk[fm];
407: lnk[i] = lnk[fm]; /* insert diagonal into linked list */
408: lnk[fm] = i;
409: lnk_lvl[i] = 0;
410: nzi++;
411: dcount++;
412: }
414: /* add pivot rows into the active row */
415: nzbd = 0;
416: prow = lnk[n];
417: while (prow < i) {
418: nnz = bdiag[prow];
419: cols = bj_ptr[prow] + nnz + 1;
420: cols_lvl = bjlvl_ptr[prow] + nnz + 1;
421: nnz = bi[prow + 1] - bi[prow] - nnz - 1;
423: PetscCall(PetscILULLAddSorted(nnz, cols, levels, cols_lvl, prow, &nlnk, lnk, lnk_lvl, lnkbt, prow));
424: nzi += nlnk;
425: prow = lnk[prow];
426: nzbd++;
427: }
428: bdiag[i] = nzbd;
429: bi[i + 1] = bi[i] + nzi;
431: /* if free space is not available, make more free space */
432: if (current_space->local_remaining < nzi) {
433: nnz = PetscIntMultTruncate(2, PetscIntMultTruncate(nzi, n - i)); /* estimated and max additional space needed */
434: PetscCall(PetscFreeSpaceGet(nnz, ¤t_space));
435: PetscCall(PetscFreeSpaceGet(nnz, ¤t_space_lvl));
436: reallocs++;
437: }
439: /* copy data into free_space and free_space_lvl, then initialize lnk */
440: PetscCall(PetscIncompleteLLClean(n, n, nzi, lnk, lnk_lvl, current_space->array, current_space_lvl->array, lnkbt));
442: bj_ptr[i] = current_space->array;
443: bjlvl_ptr[i] = current_space_lvl->array;
445: /* make sure the active row i has diagonal entry */
446: 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);
448: current_space->array += nzi;
449: current_space->local_used += nzi;
450: current_space->local_remaining -= nzi;
452: current_space_lvl->array += nzi;
453: current_space_lvl->local_used += nzi;
454: current_space_lvl->local_remaining -= nzi;
455: }
457: PetscCall(ISRestoreIndices(isrow, &r));
458: PetscCall(ISRestoreIndices(isicol, &ic));
460: /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
461: PetscCall(PetscMalloc1(bi[n] + 1, &bj));
462: PetscCall(PetscFreeSpaceContiguous_LU(&free_space, bj, n, bi, bdiag));
464: PetscCall(PetscIncompleteLLDestroy(lnk, lnkbt));
465: PetscCall(PetscFreeSpaceDestroy(free_space_lvl));
466: PetscCall(PetscFree2(bj_ptr, bjlvl_ptr));
468: #if defined(PETSC_USE_INFO)
469: {
470: PetscReal af = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);
471: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)f, (double)af));
472: PetscCall(PetscInfo(A, "Run with -[sub_]pc_factor_fill %g or use \n", (double)af));
473: PetscCall(PetscInfo(A, "PCFactorSetFill([sub]pc,%g);\n", (double)af));
474: PetscCall(PetscInfo(A, "for best performance.\n"));
475: if (diagonal_fill) PetscCall(PetscInfo(A, "Detected and replaced %" PetscInt_FMT " missing diagonals\n", dcount));
476: }
477: #endif
479: /* put together the new matrix */
480: PetscCall(MatSeqBAIJSetPreallocation(fact, bs, MAT_SKIP_ALLOCATION, NULL));
482: b = (Mat_SeqBAIJ *)fact->data;
483: b->free_ij = PETSC_TRUE;
484: PetscCall(PetscShmgetAllocateArray(bs2 * (bdiag[0] + 1), sizeof(PetscScalar), (void **)&b->a));
485: b->free_a = PETSC_TRUE;
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: PetscCall(PetscObjectReference((PetscObject)isrow));
496: PetscCall(PetscObjectReference((PetscObject)iscol));
497: b->icol = isicol;
499: PetscCall(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: PetscCall(MatSeqBAIJSetNumericFactorization(fact, both_identity));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: #if 0
513: // unused
514: /*
515: This code is virtually identical to MatILUFactorSymbolic_SeqAIJ
516: except that the data structure of Mat_SeqAIJ is slightly different.
517: Not a good example of code reuse.
518: */
519: static PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_inplace(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
520: {
521: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b;
522: IS isicol;
523: const PetscInt *r, *ic, *ai = a->i, *aj = a->j, *xi;
524: PetscInt prow, n = a->mbs, *ainew, *ajnew, jmax, *fill, nz, *im, *ajfill, *flev, *xitmp;
525: PetscInt *dloc, idx, row, m, fm, nzf, nzi, reallocate = 0, dcount = 0;
526: PetscInt incrlev, nnz, i, bs = A->rmap->bs, bs2 = a->bs2, levels, diagonal_fill, dd;
527: PetscBool col_identity, row_identity, both_identity, flg;
528: PetscReal f;
530: PetscFunctionBegin;
531: PetscCall(MatMissingDiagonal_SeqBAIJ(A, &flg, &dd));
532: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix A is missing diagonal entry in row %" PetscInt_FMT, dd);
534: f = info->fill;
535: levels = (PetscInt)info->levels;
536: diagonal_fill = (PetscInt)info->diagonal_fill;
538: PetscCall(ISInvertPermutation(iscol, PETSC_DECIDE, &isicol));
540: PetscCall(ISIdentity(isrow, &row_identity));
541: PetscCall(ISIdentity(iscol, &col_identity));
542: both_identity = (PetscBool)(row_identity && col_identity);
544: if (!levels && both_identity) { /* special case copy the nonzero structure */
545: PetscCall(MatDuplicateNoCreate_SeqBAIJ(fact, A, MAT_DO_NOT_COPY_VALUES, PETSC_TRUE));
546: PetscCall(MatSeqBAIJSetNumericFactorization_inplace(fact, both_identity));
548: fact->factortype = MAT_FACTOR_ILU;
549: b = (Mat_SeqBAIJ *)fact->data;
550: b->row = isrow;
551: b->col = iscol;
552: PetscCall(PetscObjectReference((PetscObject)isrow));
553: PetscCall(PetscObjectReference((PetscObject)iscol));
554: b->icol = isicol;
555: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
557: PetscCall(PetscMalloc1((n + 1) * bs, &b->solve_work));
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: /* general case perform the symbolic factorization */
562: PetscCall(ISGetIndices(isrow, &r));
563: PetscCall(ISGetIndices(isicol, &ic));
565: /* get new row pointers */
566: PetscCall(PetscMalloc1(n + 1, &ainew));
567: ainew[0] = 0;
568: /* don't know how many column pointers are needed so estimate */
569: jmax = (PetscInt)(f * ai[n] + 1);
570: PetscCall(PetscMalloc1(jmax, &ajnew));
571: /* ajfill is level of fill for each fill entry */
572: PetscCall(PetscMalloc1(jmax, &ajfill));
573: /* fill is a linked list of nonzeros in active row */
574: PetscCall(PetscMalloc1(n + 1, &fill));
575: /* im is level for each filled value */
576: PetscCall(PetscMalloc1(n + 1, &im));
577: /* dloc is location of diagonal in factor */
578: PetscCall(PetscMalloc1(n + 1, &dloc));
579: dloc[0] = 0;
580: for (prow = 0; prow < n; prow++) {
581: /* copy prow into linked list */
582: nzf = nz = ai[r[prow] + 1] - ai[r[prow]];
583: PetscCheck(nz, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Empty row in matrix: row in original ordering %" PetscInt_FMT " in permuted ordering %" PetscInt_FMT, r[prow], prow);
584: xi = aj + ai[r[prow]];
585: fill[n] = n;
586: fill[prow] = -1; /* marker for diagonal entry */
587: while (nz--) {
588: fm = n;
589: idx = ic[*xi++];
590: do {
591: m = fm;
592: fm = fill[m];
593: } while (fm < idx);
594: fill[m] = idx;
595: fill[idx] = fm;
596: im[idx] = 0;
597: }
599: /* make sure diagonal entry is included */
600: if (diagonal_fill && fill[prow] == -1) {
601: fm = n;
602: while (fill[fm] < prow) fm = fill[fm];
603: fill[prow] = fill[fm]; /* insert diagonal into linked list */
604: fill[fm] = prow;
605: im[prow] = 0;
606: nzf++;
607: dcount++;
608: }
610: nzi = 0;
611: row = fill[n];
612: while (row < prow) {
613: incrlev = im[row] + 1;
614: nz = dloc[row];
615: xi = ajnew + ainew[row] + nz + 1;
616: flev = ajfill + ainew[row] + nz + 1;
617: nnz = ainew[row + 1] - ainew[row] - nz - 1;
618: fm = row;
619: while (nnz-- > 0) {
620: idx = *xi++;
621: if (*flev + incrlev > levels) {
622: flev++;
623: continue;
624: }
625: do {
626: m = fm;
627: fm = fill[m];
628: } while (fm < idx);
629: if (fm != idx) {
630: im[idx] = *flev + incrlev;
631: fill[m] = idx;
632: fill[idx] = fm;
633: fm = idx;
634: nzf++;
635: } else if (im[idx] > *flev + incrlev) im[idx] = *flev + incrlev;
636: flev++;
637: }
638: row = fill[row];
639: nzi++;
640: }
641: /* copy new filled row into permanent storage */
642: ainew[prow + 1] = ainew[prow] + nzf;
643: if (ainew[prow + 1] > jmax) {
644: /* estimate how much additional space we will need */
645: /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
646: /* just double the memory each time */
647: PetscInt maxadd = jmax;
648: /* maxadd = (int)(((f*ai[n]+1)*(n-prow+5))/n); */
649: if (maxadd < nzf) maxadd = (n - prow) * (nzf + 1);
650: jmax += maxadd;
652: /* allocate a longer ajnew and ajfill */
653: PetscCall(PetscMalloc1(jmax, &xitmp));
654: PetscCall(PetscArraycpy(xitmp, ajnew, ainew[prow]));
655: PetscCall(PetscFree(ajnew));
656: ajnew = xitmp;
657: PetscCall(PetscMalloc1(jmax, &xitmp));
658: PetscCall(PetscArraycpy(xitmp, ajfill, ainew[prow]));
659: PetscCall(PetscFree(ajfill));
660: ajfill = xitmp;
661: reallocate++; /* count how many reallocations are needed */
662: }
663: xitmp = ajnew + ainew[prow];
664: flev = ajfill + ainew[prow];
665: dloc[prow] = nzi;
666: fm = fill[n];
667: while (nzf--) {
668: *xitmp++ = fm;
669: *flev++ = im[fm];
670: fm = fill[fm];
671: }
672: /* make sure row has diagonal entry */
673: PetscCheck(ajnew[ainew[prow] + dloc[prow]] == prow, 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", prow);
674: }
675: PetscCall(PetscFree(ajfill));
676: PetscCall(ISRestoreIndices(isrow, &r));
677: PetscCall(ISRestoreIndices(isicol, &ic));
678: PetscCall(PetscFree(fill));
679: PetscCall(PetscFree(im));
681: #if defined(PETSC_USE_INFO)
682: {
683: PetscReal af = ((PetscReal)ainew[n]) / ((PetscReal)ai[n]);
684: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocate, (double)f, (double)af));
685: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
686: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g);\n", (double)af));
687: PetscCall(PetscInfo(A, "for best performance.\n"));
688: if (diagonal_fill) PetscCall(PetscInfo(A, "Detected and replaced %" PetscInt_FMT " missing diagonals\n", dcount));
689: }
690: #endif
692: /* put together the new matrix */
693: PetscCall(MatSeqBAIJSetPreallocation(fact, bs, MAT_SKIP_ALLOCATION, NULL));
694: b = (Mat_SeqBAIJ *)fact->data;
695: b->free_ij = PETSC_TRUE;
696: PetscCall(PetscShmgetAllocateArray(bs2 * ainew[n],sizeof(PetscScalar),(void **)&b->a));
697: b->free_a = PETSC_TRUE;
698: b->j = ajnew;
699: b->i = ainew;
700: for (i = 0; i < n; i++) dloc[i] += ainew[i];
701: b->diag = dloc;
702: b->free_diag = PETSC_TRUE;
703: b->ilen = NULL;
704: b->imax = NULL;
705: b->row = isrow;
706: b->col = iscol;
707: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
709: PetscCall(PetscObjectReference((PetscObject)isrow));
710: PetscCall(PetscObjectReference((PetscObject)iscol));
711: b->icol = isicol;
712: PetscCall(PetscMalloc1(bs * n + bs, &b->solve_work));
713: /* In b structure: Free imax, ilen, old a, old j.
714: Allocate dloc, solve_work, new a, new j */
715: b->maxnz = b->nz = ainew[n];
717: fact->info.factor_mallocs = reallocate;
718: fact->info.fill_ratio_given = f;
719: fact->info.fill_ratio_needed = ((PetscReal)ainew[n]) / ((PetscReal)ai[prow]);
721: PetscCall(MatSeqBAIJSetNumericFactorization_inplace(fact, both_identity));
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
724: #endif
726: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE(Mat A)
727: {
728: /* Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; */
729: /* int i,*AJ=a->j,nz=a->nz; */
731: PetscFunctionBegin;
732: /* Undo Column scaling */
733: /* while (nz--) { */
734: /* AJ[i] = AJ[i]/4; */
735: /* } */
736: /* This should really invoke a push/pop logic, but we don't have that yet. */
737: A->ops->setunfactored = NULL;
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
741: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE_usj(Mat A)
742: {
743: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
744: PetscInt *AJ = a->j, nz = a->nz;
745: unsigned short *aj = (unsigned short *)AJ;
747: PetscFunctionBegin;
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: PetscFunctionReturn(PETSC_SUCCESS);
752: }