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, ¤t_space);
433: PetscFreeSpaceGet(nnz, ¤t_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: }