Actual source code: baijfact.c
1: /*
2: Factorization code for BAIJ format.
3: */
4: #include <../src/mat/impls/baij/seq/baij.h>
5: #include <petsc/private/kernels/blockinvert.h>
7: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2(Mat B, Mat A, const MatFactorInfo *info)
8: {
9: Mat C = B;
10: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
11: IS isrow = b->row, isicol = b->icol;
12: const PetscInt *r, *ic;
13: const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, bs2 = a->bs2;
14: const PetscInt *ajtmp, *bjtmp, *bdiag = b->diag;
15: MatScalar *rtmp, *pc, *mwork, *pv;
16: MatScalar *aa = a->a, *v;
17: PetscReal shift = info->shiftamount;
18: const PetscBool allowzeropivot = PetscNot(A->erroriffailure);
20: PetscFunctionBegin;
21: PetscCall(ISGetIndices(isrow, &r));
22: PetscCall(ISGetIndices(isicol, &ic));
24: /* generate work space needed by the factorization */
25: PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
26: PetscCall(PetscArrayzero(rtmp, bs2 * n));
28: for (PetscInt i = 0; i < n; i++) {
29: /* zero rtmp */
30: /* L part */
31: PetscInt *pj;
32: PetscInt nzL, nz = bi[i + 1] - bi[i];
33: bjtmp = bj + bi[i];
34: for (PetscInt j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
36: /* U part */
37: nz = bdiag[i] - bdiag[i + 1];
38: bjtmp = bj + bdiag[i + 1] + 1;
39: for (PetscInt j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
41: /* load in initial (unfactored row) */
42: nz = ai[r[i] + 1] - ai[r[i]];
43: ajtmp = aj + ai[r[i]];
44: v = aa + bs2 * ai[r[i]];
45: for (PetscInt j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ic[ajtmp[j]], v + bs2 * j, bs2));
47: /* elimination */
48: bjtmp = bj + bi[i];
49: nzL = bi[i + 1] - bi[i];
50: for (PetscInt k = 0; k < nzL; k++) {
51: PetscBool flg = PETSC_FALSE;
52: PetscInt row = bjtmp[k];
54: pc = rtmp + bs2 * row;
55: for (PetscInt j = 0; j < bs2; j++) {
56: if (pc[j] != (PetscScalar)0.0) {
57: flg = PETSC_TRUE;
58: break;
59: }
60: }
61: if (flg) {
62: pv = b->a + bs2 * bdiag[row];
63: /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
64: PetscCall(PetscKernel_A_gets_A_times_B_2(pc, pv, mwork));
66: pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
67: pv = b->a + bs2 * (bdiag[row + 1] + 1);
68: nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
69: for (PetscInt j = 0; j < nz; j++) {
70: /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
71: /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
72: v = rtmp + 4 * pj[j];
73: PetscCall(PetscKernel_A_gets_A_minus_B_times_C_2(v, pc, pv));
74: pv += 4;
75: }
76: PetscCall(PetscLogFlops(16.0 * nz + 12)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
77: }
78: }
80: /* finished row so stick it into b->a */
81: /* L part */
82: pv = b->a + bs2 * bi[i];
83: pj = b->j + bi[i];
84: nz = bi[i + 1] - bi[i];
85: for (PetscInt j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
87: /* Mark diagonal and invert diagonal for simpler triangular solves */
88: pv = b->a + bs2 * bdiag[i];
89: pj = b->j + bdiag[i];
90: PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
91: {
92: PetscBool zeropivotdetected;
94: PetscCall(PetscKernel_A_gets_inverse_A_2(pv, shift, allowzeropivot, &zeropivotdetected));
95: if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
96: }
98: /* U part */
99: pv = b->a + bs2 * (bdiag[i + 1] + 1);
100: pj = b->j + bdiag[i + 1] + 1;
101: nz = bdiag[i] - bdiag[i + 1] - 1;
102: for (PetscInt j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
103: }
105: PetscCall(PetscFree2(rtmp, mwork));
106: PetscCall(ISRestoreIndices(isicol, &ic));
107: PetscCall(ISRestoreIndices(isrow, &r));
109: C->ops->solve = MatSolve_SeqBAIJ_2;
110: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2;
111: C->assembled = PETSC_TRUE;
113: PetscCall(PetscLogFlops(1.333333333333 * 2 * 2 * 2 * n)); /* from inverting diagonal blocks */
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info)
118: {
119: Mat C = B;
120: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
121: PetscInt i, j, k, nz, nzL, row, *pj;
122: const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, bs2 = a->bs2;
123: const PetscInt *ajtmp, *bjtmp, *bdiag = b->diag;
124: MatScalar *rtmp, *pc, *mwork, *pv;
125: MatScalar *aa = a->a, *v;
126: PetscInt flg;
127: PetscReal shift = info->shiftamount;
128: PetscBool allowzeropivot, zeropivotdetected;
130: PetscFunctionBegin;
131: allowzeropivot = PetscNot(A->erroriffailure);
133: /* generate work space needed by the factorization */
134: PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
135: PetscCall(PetscArrayzero(rtmp, bs2 * n));
137: for (i = 0; i < n; i++) {
138: /* zero rtmp */
139: /* L part */
140: nz = bi[i + 1] - bi[i];
141: bjtmp = bj + bi[i];
142: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
144: /* U part */
145: nz = bdiag[i] - bdiag[i + 1];
146: bjtmp = bj + bdiag[i + 1] + 1;
147: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
149: /* load in initial (unfactored row) */
150: nz = ai[i + 1] - ai[i];
151: ajtmp = aj + ai[i];
152: v = aa + bs2 * ai[i];
153: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ajtmp[j], v + bs2 * j, bs2));
155: /* elimination */
156: bjtmp = bj + bi[i];
157: nzL = bi[i + 1] - bi[i];
158: for (k = 0; k < nzL; k++) {
159: row = bjtmp[k];
160: pc = rtmp + bs2 * row;
161: for (flg = 0, j = 0; j < bs2; j++) {
162: if (pc[j] != (PetscScalar)0.0) {
163: flg = 1;
164: break;
165: }
166: }
167: if (flg) {
168: pv = b->a + bs2 * bdiag[row];
169: /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
170: PetscCall(PetscKernel_A_gets_A_times_B_2(pc, pv, mwork));
172: pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
173: pv = b->a + bs2 * (bdiag[row + 1] + 1);
174: nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries in U(row,:) excluding diag */
175: for (j = 0; j < nz; j++) {
176: /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
177: /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
178: v = rtmp + 4 * pj[j];
179: PetscCall(PetscKernel_A_gets_A_minus_B_times_C_2(v, pc, pv));
180: pv += 4;
181: }
182: PetscCall(PetscLogFlops(16.0 * nz + 12)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
183: }
184: }
186: /* finished row so stick it into b->a */
187: /* L part */
188: pv = b->a + bs2 * bi[i];
189: pj = b->j + bi[i];
190: nz = bi[i + 1] - bi[i];
191: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
193: /* Mark diagonal and invert diagonal for simpler triangular solves */
194: pv = b->a + bs2 * bdiag[i];
195: pj = b->j + bdiag[i];
196: PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
197: PetscCall(PetscKernel_A_gets_inverse_A_2(pv, shift, allowzeropivot, &zeropivotdetected));
198: if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
200: /* U part */
201: /*
202: pv = b->a + bs2*bi[2*n-i];
203: pj = b->j + bi[2*n-i];
204: nz = bi[2*n-i+1] - bi[2*n-i] - 1;
205: */
206: pv = b->a + bs2 * (bdiag[i + 1] + 1);
207: pj = b->j + bdiag[i + 1] + 1;
208: nz = bdiag[i] - bdiag[i + 1] - 1;
209: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
210: }
211: PetscCall(PetscFree2(rtmp, mwork));
213: C->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering;
214: C->ops->forwardsolve = MatForwardSolve_SeqBAIJ_2_NaturalOrdering;
215: C->ops->backwardsolve = MatBackwardSolve_SeqBAIJ_2_NaturalOrdering;
216: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering;
217: C->assembled = PETSC_TRUE;
219: PetscCall(PetscLogFlops(1.333333333333 * 2 * 2 * 2 * n)); /* from inverting diagonal blocks */
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_inplace(Mat B, Mat A, const MatFactorInfo *info)
224: {
225: Mat C = B;
226: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
227: IS isrow = b->row, isicol = b->icol;
228: const PetscInt *r, *ic;
229: PetscInt i, j, n = a->mbs, *bi = b->i, *bj = b->j;
230: PetscInt *ajtmpold, *ajtmp, nz, row;
231: PetscInt *diag_offset = b->diag, idx, *ai = a->i, *aj = a->j, *pj;
232: MatScalar *pv, *v, *rtmp, m1, m2, m3, m4, *pc, *w, *x, x1, x2, x3, x4;
233: MatScalar p1, p2, p3, p4;
234: MatScalar *ba = b->a, *aa = a->a;
235: PetscReal shift = info->shiftamount;
236: PetscBool allowzeropivot, zeropivotdetected;
238: PetscFunctionBegin;
239: allowzeropivot = PetscNot(A->erroriffailure);
240: PetscCall(ISGetIndices(isrow, &r));
241: PetscCall(ISGetIndices(isicol, &ic));
242: PetscCall(PetscMalloc1(4 * (n + 1), &rtmp));
244: for (i = 0; i < n; i++) {
245: nz = bi[i + 1] - bi[i];
246: ajtmp = bj + bi[i];
247: for (j = 0; j < nz; j++) {
248: x = rtmp + 4 * ajtmp[j];
249: x[0] = x[1] = x[2] = x[3] = 0.0;
250: }
251: /* load in initial (unfactored row) */
252: idx = r[i];
253: nz = ai[idx + 1] - ai[idx];
254: ajtmpold = aj + ai[idx];
255: v = aa + 4 * ai[idx];
256: for (j = 0; j < nz; j++) {
257: x = rtmp + 4 * ic[ajtmpold[j]];
258: x[0] = v[0];
259: x[1] = v[1];
260: x[2] = v[2];
261: x[3] = v[3];
262: v += 4;
263: }
264: row = *ajtmp++;
265: while (row < i) {
266: pc = rtmp + 4 * row;
267: p1 = pc[0];
268: p2 = pc[1];
269: p3 = pc[2];
270: p4 = pc[3];
271: if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) {
272: pv = ba + 4 * diag_offset[row];
273: pj = bj + diag_offset[row] + 1;
274: x1 = pv[0];
275: x2 = pv[1];
276: x3 = pv[2];
277: x4 = pv[3];
278: pc[0] = m1 = p1 * x1 + p3 * x2;
279: pc[1] = m2 = p2 * x1 + p4 * x2;
280: pc[2] = m3 = p1 * x3 + p3 * x4;
281: pc[3] = m4 = p2 * x3 + p4 * x4;
282: nz = bi[row + 1] - diag_offset[row] - 1;
283: pv += 4;
284: for (j = 0; j < nz; j++) {
285: x1 = pv[0];
286: x2 = pv[1];
287: x3 = pv[2];
288: x4 = pv[3];
289: x = rtmp + 4 * pj[j];
290: x[0] -= m1 * x1 + m3 * x2;
291: x[1] -= m2 * x1 + m4 * x2;
292: x[2] -= m1 * x3 + m3 * x4;
293: x[3] -= m2 * x3 + m4 * x4;
294: pv += 4;
295: }
296: PetscCall(PetscLogFlops(16.0 * nz + 12.0));
297: }
298: row = *ajtmp++;
299: }
300: /* finished row so stick it into b->a */
301: pv = ba + 4 * bi[i];
302: pj = bj + bi[i];
303: nz = bi[i + 1] - bi[i];
304: for (j = 0; j < nz; j++) {
305: x = rtmp + 4 * pj[j];
306: pv[0] = x[0];
307: pv[1] = x[1];
308: pv[2] = x[2];
309: pv[3] = x[3];
310: pv += 4;
311: }
312: /* invert diagonal block */
313: w = ba + 4 * diag_offset[i];
314: PetscCall(PetscKernel_A_gets_inverse_A_2(w, shift, allowzeropivot, &zeropivotdetected));
315: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
316: }
318: PetscCall(PetscFree(rtmp));
319: PetscCall(ISRestoreIndices(isicol, &ic));
320: PetscCall(ISRestoreIndices(isrow, &r));
322: C->ops->solve = MatSolve_SeqBAIJ_2_inplace;
323: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_inplace;
324: C->assembled = PETSC_TRUE;
326: PetscCall(PetscLogFlops(1.333333333333 * 8 * b->mbs)); /* from inverting diagonal blocks */
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
329: /*
330: Version for when blocks are 2 by 2 Using natural ordering
331: */
332: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering_inplace(Mat C, Mat A, const MatFactorInfo *info)
333: {
334: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
335: PetscInt i, j, n = a->mbs, *bi = b->i, *bj = b->j;
336: PetscInt *ajtmpold, *ajtmp, nz, row;
337: PetscInt *diag_offset = b->diag, *ai = a->i, *aj = a->j, *pj;
338: MatScalar *pv, *v, *rtmp, *pc, *w, *x;
339: MatScalar p1, p2, p3, p4, m1, m2, m3, m4, x1, x2, x3, x4;
340: MatScalar *ba = b->a, *aa = a->a;
341: PetscReal shift = info->shiftamount;
342: PetscBool allowzeropivot, zeropivotdetected;
344: PetscFunctionBegin;
345: allowzeropivot = PetscNot(A->erroriffailure);
346: PetscCall(PetscMalloc1(4 * (n + 1), &rtmp));
347: for (i = 0; i < n; i++) {
348: nz = bi[i + 1] - bi[i];
349: ajtmp = bj + bi[i];
350: for (j = 0; j < nz; j++) {
351: x = rtmp + 4 * ajtmp[j];
352: x[0] = x[1] = x[2] = x[3] = 0.0;
353: }
354: /* load in initial (unfactored row) */
355: nz = ai[i + 1] - ai[i];
356: ajtmpold = aj + ai[i];
357: v = aa + 4 * ai[i];
358: for (j = 0; j < nz; j++) {
359: x = rtmp + 4 * ajtmpold[j];
360: x[0] = v[0];
361: x[1] = v[1];
362: x[2] = v[2];
363: x[3] = v[3];
364: v += 4;
365: }
366: row = *ajtmp++;
367: while (row < i) {
368: pc = rtmp + 4 * row;
369: p1 = pc[0];
370: p2 = pc[1];
371: p3 = pc[2];
372: p4 = pc[3];
373: if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) {
374: pv = ba + 4 * diag_offset[row];
375: pj = bj + diag_offset[row] + 1;
376: x1 = pv[0];
377: x2 = pv[1];
378: x3 = pv[2];
379: x4 = pv[3];
380: pc[0] = m1 = p1 * x1 + p3 * x2;
381: pc[1] = m2 = p2 * x1 + p4 * x2;
382: pc[2] = m3 = p1 * x3 + p3 * x4;
383: pc[3] = m4 = p2 * x3 + p4 * x4;
384: nz = bi[row + 1] - diag_offset[row] - 1;
385: pv += 4;
386: for (j = 0; j < nz; j++) {
387: x1 = pv[0];
388: x2 = pv[1];
389: x3 = pv[2];
390: x4 = pv[3];
391: x = rtmp + 4 * pj[j];
392: x[0] -= m1 * x1 + m3 * x2;
393: x[1] -= m2 * x1 + m4 * x2;
394: x[2] -= m1 * x3 + m3 * x4;
395: x[3] -= m2 * x3 + m4 * x4;
396: pv += 4;
397: }
398: PetscCall(PetscLogFlops(16.0 * nz + 12.0));
399: }
400: row = *ajtmp++;
401: }
402: /* finished row so stick it into b->a */
403: pv = ba + 4 * bi[i];
404: pj = bj + bi[i];
405: nz = bi[i + 1] - bi[i];
406: for (j = 0; j < nz; j++) {
407: x = rtmp + 4 * pj[j];
408: pv[0] = x[0];
409: pv[1] = x[1];
410: pv[2] = x[2];
411: pv[3] = x[3];
412: /*
413: printf(" col %d:",pj[j]);
414: PetscInt j1;
415: for (j1=0; j1<4; j1++) printf(" %g,",*(pv+j1));
416: printf("\n");
417: */
418: pv += 4;
419: }
420: /* invert diagonal block */
421: w = ba + 4 * diag_offset[i];
422: PetscCall(PetscKernel_A_gets_inverse_A_2(w, shift, allowzeropivot, &zeropivotdetected));
423: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
424: }
426: PetscCall(PetscFree(rtmp));
428: C->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering_inplace;
429: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering_inplace;
430: C->assembled = PETSC_TRUE;
432: PetscCall(PetscLogFlops(1.333333333333 * 8 * b->mbs)); /* from inverting diagonal blocks */
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: /*
437: Version for when blocks are 1 by 1.
438: */
439: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1(Mat B, Mat A, const MatFactorInfo *info)
440: {
441: Mat C = B;
442: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
443: IS isrow = b->row, isicol = b->icol;
444: const PetscInt *r, *ic, *ics;
445: const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *bdiag = b->diag;
446: PetscInt i, j, k, nz, nzL, row, *pj;
447: const PetscInt *ajtmp, *bjtmp;
448: MatScalar *rtmp, *pc, multiplier, *pv;
449: const MatScalar *aa = a->a, *v;
450: PetscBool row_identity, col_identity;
451: FactorShiftCtx sctx;
452: const PetscInt *ddiag;
453: PetscReal rs;
454: MatScalar d;
456: PetscFunctionBegin;
457: /* MatPivotSetUp(): initialize shift context sctx */
458: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
460: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
461: ddiag = a->diag;
462: sctx.shift_top = info->zeropivot;
463: for (i = 0; i < n; i++) {
464: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
465: d = (aa)[ddiag[i]];
466: rs = -PetscAbsScalar(d) - PetscRealPart(d);
467: v = aa + ai[i];
468: nz = ai[i + 1] - ai[i];
469: for (j = 0; j < nz; j++) rs += PetscAbsScalar(v[j]);
470: if (rs > sctx.shift_top) sctx.shift_top = rs;
471: }
472: sctx.shift_top *= 1.1;
473: sctx.nshift_max = 5;
474: sctx.shift_lo = 0.;
475: sctx.shift_hi = 1.;
476: }
478: PetscCall(ISGetIndices(isrow, &r));
479: PetscCall(ISGetIndices(isicol, &ic));
480: PetscCall(PetscMalloc1(n + 1, &rtmp));
481: ics = ic;
483: do {
484: sctx.newshift = PETSC_FALSE;
485: for (i = 0; i < n; i++) {
486: /* zero rtmp */
487: /* L part */
488: nz = bi[i + 1] - bi[i];
489: bjtmp = bj + bi[i];
490: for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;
492: /* U part */
493: nz = bdiag[i] - bdiag[i + 1];
494: bjtmp = bj + bdiag[i + 1] + 1;
495: for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;
497: /* load in initial (unfactored row) */
498: nz = ai[r[i] + 1] - ai[r[i]];
499: ajtmp = aj + ai[r[i]];
500: v = aa + ai[r[i]];
501: for (j = 0; j < nz; j++) rtmp[ics[ajtmp[j]]] = v[j];
503: /* ZeropivotApply() */
504: rtmp[i] += sctx.shift_amount; /* shift the diagonal of the matrix */
506: /* elimination */
507: bjtmp = bj + bi[i];
508: row = *bjtmp++;
509: nzL = bi[i + 1] - bi[i];
510: for (k = 0; k < nzL; k++) {
511: pc = rtmp + row;
512: if (*pc != (PetscScalar)0.0) {
513: pv = b->a + bdiag[row];
514: multiplier = *pc * (*pv);
515: *pc = multiplier;
517: pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
518: pv = b->a + bdiag[row + 1] + 1;
519: nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries in U(row,:) excluding diag */
520: for (j = 0; j < nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
521: PetscCall(PetscLogFlops(2.0 * nz));
522: }
523: row = *bjtmp++;
524: }
526: /* finished row so stick it into b->a */
527: rs = 0.0;
528: /* L part */
529: pv = b->a + bi[i];
530: pj = b->j + bi[i];
531: nz = bi[i + 1] - bi[i];
532: for (j = 0; j < nz; j++) {
533: pv[j] = rtmp[pj[j]];
534: rs += PetscAbsScalar(pv[j]);
535: }
537: /* U part */
538: pv = b->a + bdiag[i + 1] + 1;
539: pj = b->j + bdiag[i + 1] + 1;
540: nz = bdiag[i] - bdiag[i + 1] - 1;
541: for (j = 0; j < nz; j++) {
542: pv[j] = rtmp[pj[j]];
543: rs += PetscAbsScalar(pv[j]);
544: }
546: sctx.rs = rs;
547: sctx.pv = rtmp[i];
548: PetscCall(MatPivotCheck(B, A, info, &sctx, i));
549: if (sctx.newshift) break; /* break for-loop */
550: rtmp[i] = sctx.pv; /* sctx.pv might be updated in the case of MAT_SHIFT_INBLOCKS */
552: /* Mark diagonal and invert diagonal for simpler triangular solves */
553: pv = b->a + bdiag[i];
554: *pv = (PetscScalar)1.0 / rtmp[i];
556: } /* endof for (i=0; i<n; i++) { */
558: /* MatPivotRefine() */
559: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction > 0 && sctx.nshift < sctx.nshift_max) {
560: /*
561: * if no shift in this attempt & shifting & started shifting & can refine,
562: * then try lower shift
563: */
564: sctx.shift_hi = sctx.shift_fraction;
565: sctx.shift_fraction = (sctx.shift_hi + sctx.shift_lo) / 2.;
566: sctx.shift_amount = sctx.shift_fraction * sctx.shift_top;
567: sctx.newshift = PETSC_TRUE;
568: sctx.nshift++;
569: }
570: } while (sctx.newshift);
572: PetscCall(PetscFree(rtmp));
573: PetscCall(ISRestoreIndices(isicol, &ic));
574: PetscCall(ISRestoreIndices(isrow, &r));
576: PetscCall(ISIdentity(isrow, &row_identity));
577: PetscCall(ISIdentity(isicol, &col_identity));
578: if (row_identity && col_identity) {
579: C->ops->solve = MatSolve_SeqBAIJ_1_NaturalOrdering;
580: C->ops->forwardsolve = MatForwardSolve_SeqBAIJ_1_NaturalOrdering;
581: C->ops->backwardsolve = MatBackwardSolve_SeqBAIJ_1_NaturalOrdering;
582: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering;
583: } else {
584: C->ops->solve = MatSolve_SeqBAIJ_1;
585: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1;
586: }
587: C->assembled = PETSC_TRUE;
588: PetscCall(PetscLogFlops(C->cmap->n));
590: /* MatShiftView(A,info,&sctx) */
591: if (sctx.nshift) {
592: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
593: PetscCall(PetscInfo(A, "number of shift_pd tries %" PetscInt_FMT ", shift_amount %g, diagonal shifted up by %e fraction top_value %e\n", sctx.nshift, (double)sctx.shift_amount, (double)sctx.shift_fraction, (double)sctx.shift_top));
594: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
595: PetscCall(PetscInfo(A, "number of shift_nz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
596: } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
597: PetscCall(PetscInfo(A, "number of shift_inblocks applied %" PetscInt_FMT ", each shift_amount %g\n", sctx.nshift, (double)info->shiftamount));
598: }
599: }
600: PetscFunctionReturn(PETSC_SUCCESS);
601: }
603: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1_inplace(Mat C, Mat A, const MatFactorInfo *info)
604: {
605: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
606: IS isrow = b->row, isicol = b->icol;
607: const PetscInt *r, *ic;
608: PetscInt i, j, n = a->mbs, *bi = b->i, *bj = b->j;
609: PetscInt *ajtmpold, *ajtmp, nz, row, *ai = a->i, *aj = a->j;
610: PetscInt *diag_offset = b->diag, diag, *pj;
611: MatScalar *pv, *v, *rtmp, multiplier, *pc;
612: MatScalar *ba = b->a, *aa = a->a;
613: PetscBool row_identity, col_identity;
615: PetscFunctionBegin;
616: PetscCall(ISGetIndices(isrow, &r));
617: PetscCall(ISGetIndices(isicol, &ic));
618: PetscCall(PetscMalloc1(n + 1, &rtmp));
620: for (i = 0; i < n; i++) {
621: nz = bi[i + 1] - bi[i];
622: ajtmp = bj + bi[i];
623: for (j = 0; j < nz; j++) rtmp[ajtmp[j]] = 0.0;
625: /* load in initial (unfactored row) */
626: nz = ai[r[i] + 1] - ai[r[i]];
627: ajtmpold = aj + ai[r[i]];
628: v = aa + ai[r[i]];
629: for (j = 0; j < nz; j++) rtmp[ic[ajtmpold[j]]] = v[j];
631: row = *ajtmp++;
632: while (row < i) {
633: pc = rtmp + row;
634: if (*pc != 0.0) {
635: pv = ba + diag_offset[row];
636: pj = bj + diag_offset[row] + 1;
637: multiplier = *pc * *pv++;
638: *pc = multiplier;
639: nz = bi[row + 1] - diag_offset[row] - 1;
640: for (j = 0; j < nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
641: PetscCall(PetscLogFlops(1.0 + 2.0 * nz));
642: }
643: row = *ajtmp++;
644: }
645: /* finished row so stick it into b->a */
646: pv = ba + bi[i];
647: pj = bj + bi[i];
648: nz = bi[i + 1] - bi[i];
649: for (j = 0; j < nz; j++) pv[j] = rtmp[pj[j]];
650: diag = diag_offset[i] - bi[i];
651: /* check pivot entry for current row */
652: PetscCheck(pv[diag] != 0.0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot: row in original ordering %" PetscInt_FMT " in permuted ordering %" PetscInt_FMT, r[i], i);
653: pv[diag] = 1.0 / pv[diag];
654: }
656: PetscCall(PetscFree(rtmp));
657: PetscCall(ISRestoreIndices(isicol, &ic));
658: PetscCall(ISRestoreIndices(isrow, &r));
659: PetscCall(ISIdentity(isrow, &row_identity));
660: PetscCall(ISIdentity(isicol, &col_identity));
661: if (row_identity && col_identity) {
662: C->ops->solve = MatSolve_SeqBAIJ_1_NaturalOrdering_inplace;
663: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering_inplace;
664: } else {
665: C->ops->solve = MatSolve_SeqBAIJ_1_inplace;
666: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_inplace;
667: }
668: C->assembled = PETSC_TRUE;
669: PetscCall(PetscLogFlops(C->cmap->n));
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: static PetscErrorCode MatFactorGetSolverType_petsc(Mat A, MatSolverType *type)
674: {
675: PetscFunctionBegin;
676: *type = MATSOLVERPETSC;
677: PetscFunctionReturn(PETSC_SUCCESS);
678: }
680: PETSC_INTERN PetscErrorCode MatGetFactor_seqbaij_petsc(Mat A, MatFactorType ftype, Mat *B)
681: {
682: PetscInt n = A->rmap->n;
684: PetscFunctionBegin;
685: #if defined(PETSC_USE_COMPLEX)
686: PetscCheck(A->hermitian != PETSC_BOOL3_TRUE || !(ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC), PETSC_COMM_SELF, PETSC_ERR_SUP, "Hermitian Factor is not supported");
687: #endif
688: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
689: PetscCall(MatSetSizes(*B, n, n, n, n));
690: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU || ftype == MAT_FACTOR_ILUDT) {
691: PetscCall(MatSetType(*B, MATSEQBAIJ));
693: (*B)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqBAIJ;
694: (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqBAIJ;
695: PetscCall(PetscStrallocpy(MATORDERINGND, (char **)&(*B)->preferredordering[MAT_FACTOR_LU]));
696: PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ILU]));
697: PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ILUDT]));
698: } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
699: PetscCall(MatSetType(*B, MATSEQSBAIJ));
700: PetscCall(MatSeqSBAIJSetPreallocation(*B, A->rmap->bs, MAT_SKIP_ALLOCATION, NULL));
702: (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqBAIJ;
703: (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqBAIJ;
704: /* Future optimization would be direct symbolic and numerical factorization for BAIJ to support orderings and Cholesky, instead of first converting to SBAIJ */
705: PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]));
706: PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ICC]));
707: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");
708: (*B)->factortype = ftype;
709: (*B)->canuseordering = PETSC_TRUE;
711: PetscCall(PetscFree((*B)->solvertype));
712: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
713: PetscCall(PetscObjectComposeFunction((PetscObject)*B, "MatFactorGetSolverType_C", MatFactorGetSolverType_petsc));
714: PetscFunctionReturn(PETSC_SUCCESS);
715: }
717: PetscErrorCode MatLUFactor_SeqBAIJ(Mat A, IS row, IS col, const MatFactorInfo *info)
718: {
719: Mat C;
721: PetscFunctionBegin;
722: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &C));
723: PetscCall(MatLUFactorSymbolic(C, A, row, col, info));
724: PetscCall(MatLUFactorNumeric(C, A, info));
726: A->ops->solve = C->ops->solve;
727: A->ops->solvetranspose = C->ops->solvetranspose;
729: PetscCall(MatHeaderMerge(A, &C));
730: PetscFunctionReturn(PETSC_SUCCESS);
731: }
733: #include <../src/mat/impls/sbaij/seq/sbaij.h>
734: PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N(Mat C, Mat A, const MatFactorInfo *info)
735: {
736: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
737: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)C->data;
738: IS ip = b->row;
739: const PetscInt *rip;
740: PetscInt i, j, mbs = a->mbs, bs = A->rmap->bs, *bi = b->i, *bj = b->j, *bcol;
741: PetscInt *ai = a->i, *aj = a->j;
742: PetscInt k, jmin, jmax, *jl, *il, col, nexti, ili, nz;
743: MatScalar *rtmp, *ba = b->a, *bval, *aa = a->a, dk, uikdi;
744: PetscReal rs;
745: FactorShiftCtx sctx;
747: PetscFunctionBegin;
748: if (bs > 1) { /* convert A to a SBAIJ matrix and apply Cholesky factorization from it */
749: if (!a->sbaijMat) PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &a->sbaijMat));
750: PetscCall(a->sbaijMat->ops->choleskyfactornumeric(C, a->sbaijMat, info));
751: PetscCall(MatDestroy(&a->sbaijMat));
752: PetscFunctionReturn(PETSC_SUCCESS);
753: }
755: /* MatPivotSetUp(): initialize shift context sctx */
756: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
758: PetscCall(ISGetIndices(ip, &rip));
759: PetscCall(PetscMalloc3(mbs, &rtmp, mbs, &il, mbs, &jl));
761: sctx.shift_amount = 0.;
762: sctx.nshift = 0;
763: do {
764: sctx.newshift = PETSC_FALSE;
765: for (i = 0; i < mbs; i++) {
766: rtmp[i] = 0.0;
767: jl[i] = mbs;
768: il[0] = 0;
769: }
771: for (k = 0; k < mbs; k++) {
772: bval = ba + bi[k];
773: /* initialize k-th row by the perm[k]-th row of A */
774: jmin = ai[rip[k]];
775: jmax = ai[rip[k] + 1];
776: for (j = jmin; j < jmax; j++) {
777: col = rip[aj[j]];
778: if (col >= k) { /* only take upper triangular entry */
779: rtmp[col] = aa[j];
780: *bval++ = 0.0; /* for in-place factorization */
781: }
782: }
784: /* shift the diagonal of the matrix */
785: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
787: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
788: dk = rtmp[k];
789: i = jl[k]; /* first row to be added to k_th row */
791: while (i < k) {
792: nexti = jl[i]; /* next row to be added to k_th row */
794: /* compute multiplier, update diag(k) and U(i,k) */
795: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
796: uikdi = -ba[ili] * ba[bi[i]]; /* diagonal(k) */
797: dk += uikdi * ba[ili];
798: ba[ili] = uikdi; /* -U(i,k) */
800: /* add multiple of row i to k-th row */
801: jmin = ili + 1;
802: jmax = bi[i + 1];
803: if (jmin < jmax) {
804: for (j = jmin; j < jmax; j++) rtmp[bj[j]] += uikdi * ba[j];
805: /* update il and jl for row i */
806: il[i] = jmin;
807: j = bj[jmin];
808: jl[i] = jl[j];
809: jl[j] = i;
810: }
811: i = nexti;
812: }
814: /* shift the diagonals when zero pivot is detected */
815: /* compute rs=sum of abs(off-diagonal) */
816: rs = 0.0;
817: jmin = bi[k] + 1;
818: nz = bi[k + 1] - jmin;
819: if (nz) {
820: bcol = bj + jmin;
821: while (nz--) {
822: rs += PetscAbsScalar(rtmp[*bcol]);
823: bcol++;
824: }
825: }
827: sctx.rs = rs;
828: sctx.pv = dk;
829: PetscCall(MatPivotCheck(C, A, info, &sctx, k));
830: if (sctx.newshift) break;
831: dk = sctx.pv;
833: /* copy data into U(k,:) */
834: ba[bi[k]] = 1.0 / dk; /* U(k,k) */
835: jmin = bi[k] + 1;
836: jmax = bi[k + 1];
837: if (jmin < jmax) {
838: for (j = jmin; j < jmax; j++) {
839: col = bj[j];
840: ba[j] = rtmp[col];
841: rtmp[col] = 0.0;
842: }
843: /* add the k-th row into il and jl */
844: il[k] = jmin;
845: i = bj[jmin];
846: jl[k] = jl[i];
847: jl[i] = k;
848: }
849: }
850: } while (sctx.newshift);
851: PetscCall(PetscFree3(rtmp, il, jl));
853: PetscCall(ISRestoreIndices(ip, &rip));
855: C->assembled = PETSC_TRUE;
856: C->preallocated = PETSC_TRUE;
858: PetscCall(PetscLogFlops(C->rmap->N));
859: if (sctx.nshift) {
860: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
861: PetscCall(PetscInfo(A, "number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
862: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
863: PetscCall(PetscInfo(A, "number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
864: }
865: }
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering(Mat C, Mat A, const MatFactorInfo *info)
870: {
871: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
872: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)C->data;
873: PetscInt i, j, am = a->mbs;
874: PetscInt *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
875: PetscInt k, jmin, *jl, *il, nexti, ili, *acol, *bcol, nz;
876: MatScalar *rtmp, *ba = b->a, *aa = a->a, dk, uikdi, *aval, *bval;
877: PetscReal rs;
878: FactorShiftCtx sctx;
880: PetscFunctionBegin;
881: /* MatPivotSetUp(): initialize shift context sctx */
882: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
884: PetscCall(PetscMalloc3(am, &rtmp, am, &il, am, &jl));
886: do {
887: sctx.newshift = PETSC_FALSE;
888: for (i = 0; i < am; i++) {
889: rtmp[i] = 0.0;
890: jl[i] = am;
891: il[0] = 0;
892: }
894: for (k = 0; k < am; k++) {
895: /* initialize k-th row with elements nonzero in row perm(k) of A */
896: nz = ai[k + 1] - ai[k];
897: acol = aj + ai[k];
898: aval = aa + ai[k];
899: bval = ba + bi[k];
900: while (nz--) {
901: if (*acol < k) { /* skip lower triangular entries */
902: acol++;
903: aval++;
904: } else {
905: rtmp[*acol++] = *aval++;
906: *bval++ = 0.0; /* for in-place factorization */
907: }
908: }
910: /* shift the diagonal of the matrix */
911: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
913: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
914: dk = rtmp[k];
915: i = jl[k]; /* first row to be added to k_th row */
917: while (i < k) {
918: nexti = jl[i]; /* next row to be added to k_th row */
919: /* compute multiplier, update D(k) and U(i,k) */
920: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
921: uikdi = -ba[ili] * ba[bi[i]];
922: dk += uikdi * ba[ili];
923: ba[ili] = uikdi; /* -U(i,k) */
925: /* add multiple of row i to k-th row ... */
926: jmin = ili + 1;
927: nz = bi[i + 1] - jmin;
928: if (nz > 0) {
929: bcol = bj + jmin;
930: bval = ba + jmin;
931: while (nz--) rtmp[*bcol++] += uikdi * (*bval++);
932: /* update il and jl for i-th row */
933: il[i] = jmin;
934: j = bj[jmin];
935: jl[i] = jl[j];
936: jl[j] = i;
937: }
938: i = nexti;
939: }
941: /* shift the diagonals when zero pivot is detected */
942: /* compute rs=sum of abs(off-diagonal) */
943: rs = 0.0;
944: jmin = bi[k] + 1;
945: nz = bi[k + 1] - jmin;
946: if (nz) {
947: bcol = bj + jmin;
948: while (nz--) {
949: rs += PetscAbsScalar(rtmp[*bcol]);
950: bcol++;
951: }
952: }
954: sctx.rs = rs;
955: sctx.pv = dk;
956: PetscCall(MatPivotCheck(C, A, info, &sctx, k));
957: if (sctx.newshift) break; /* sctx.shift_amount is updated */
958: dk = sctx.pv;
960: /* copy data into U(k,:) */
961: ba[bi[k]] = 1.0 / dk;
962: jmin = bi[k] + 1;
963: nz = bi[k + 1] - jmin;
964: if (nz) {
965: bcol = bj + jmin;
966: bval = ba + jmin;
967: while (nz--) {
968: *bval++ = rtmp[*bcol];
969: rtmp[*bcol++] = 0.0;
970: }
971: /* add k-th row into il and jl */
972: il[k] = jmin;
973: i = bj[jmin];
974: jl[k] = jl[i];
975: jl[i] = k;
976: }
977: }
978: } while (sctx.newshift);
979: PetscCall(PetscFree3(rtmp, il, jl));
981: C->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
982: C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
983: C->assembled = PETSC_TRUE;
984: C->preallocated = PETSC_TRUE;
986: PetscCall(PetscLogFlops(C->rmap->N));
987: if (sctx.nshift) {
988: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
989: PetscCall(PetscInfo(A, "number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
990: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
991: PetscCall(PetscInfo(A, "number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
992: }
993: }
994: PetscFunctionReturn(PETSC_SUCCESS);
995: }
997: #include <petscbt.h>
998: #include <../src/mat/utils/freespace.h>
999: PetscErrorCode MatICCFactorSymbolic_SeqBAIJ(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
1000: {
1001: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1002: Mat_SeqSBAIJ *b;
1003: Mat B;
1004: PetscBool perm_identity, missing;
1005: PetscInt reallocs = 0, i, *ai = a->i, *aj = a->j, am = a->mbs, bs = A->rmap->bs, *ui;
1006: const PetscInt *rip;
1007: PetscInt jmin, jmax, nzk, k, j, *jl, prow, *il, nextprow;
1008: PetscInt nlnk, *lnk, *lnk_lvl = NULL, ncols, ncols_upper, *cols, *cols_lvl, *uj, **uj_ptr, **uj_lvl_ptr;
1009: PetscReal fill = info->fill, levels = info->levels;
1010: PetscFreeSpaceList free_space = NULL, current_space = NULL;
1011: PetscFreeSpaceList free_space_lvl = NULL, current_space_lvl = NULL;
1012: PetscBT lnkbt;
1014: PetscFunctionBegin;
1015: PetscCall(MatMissingDiagonal(A, &missing, &i));
1016: PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry %" PetscInt_FMT, i);
1018: if (bs > 1) {
1019: if (!a->sbaijMat) PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &a->sbaijMat));
1020: fact->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */
1022: PetscCall(MatICCFactorSymbolic(fact, a->sbaijMat, perm, info));
1023: PetscFunctionReturn(PETSC_SUCCESS);
1024: }
1026: PetscCall(ISIdentity(perm, &perm_identity));
1027: PetscCall(ISGetIndices(perm, &rip));
1029: /* special case that simply copies fill pattern */
1030: if (!levels && perm_identity) {
1031: PetscCall(PetscMalloc1(am + 1, &ui));
1032: for (i = 0; i < am; i++) ui[i] = ai[i + 1] - a->diag[i]; /* ui: rowlengths - changes when !perm_identity */
1033: B = fact;
1034: PetscCall(MatSeqSBAIJSetPreallocation(B, 1, 0, ui));
1036: b = (Mat_SeqSBAIJ *)B->data;
1037: uj = b->j;
1038: for (i = 0; i < am; i++) {
1039: aj = a->j + a->diag[i];
1040: for (j = 0; j < ui[i]; j++) *uj++ = *aj++;
1041: b->ilen[i] = ui[i];
1042: }
1043: PetscCall(PetscFree(ui));
1045: B->factortype = MAT_FACTOR_NONE;
1047: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1048: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1049: B->factortype = MAT_FACTOR_ICC;
1051: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1052: PetscFunctionReturn(PETSC_SUCCESS);
1053: }
1055: /* initialization */
1056: PetscCall(PetscMalloc1(am + 1, &ui));
1057: ui[0] = 0;
1058: PetscCall(PetscMalloc1(2 * am + 1, &cols_lvl));
1060: /* jl: linked list for storing indices of the pivot rows
1061: il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
1062: PetscCall(PetscMalloc4(am, &uj_ptr, am, &uj_lvl_ptr, am, &il, am, &jl));
1063: for (i = 0; i < am; i++) {
1064: jl[i] = am;
1065: il[i] = 0;
1066: }
1068: /* create and initialize a linked list for storing column indices of the active row k */
1069: nlnk = am + 1;
1070: PetscCall(PetscIncompleteLLCreate(am, am, nlnk, lnk, lnk_lvl, lnkbt));
1072: /* initial FreeSpace size is fill*(ai[am]+am)/2 */
1073: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am] / 2, am / 2)), &free_space));
1075: current_space = free_space;
1077: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am] / 2, am / 2)), &free_space_lvl));
1078: current_space_lvl = free_space_lvl;
1080: for (k = 0; k < am; k++) { /* for each active row k */
1081: /* initialize lnk by the column indices of row rip[k] of A */
1082: nzk = 0;
1083: ncols = ai[rip[k] + 1] - ai[rip[k]];
1084: ncols_upper = 0;
1085: cols = cols_lvl + am;
1086: for (j = 0; j < ncols; j++) {
1087: i = rip[*(aj + ai[rip[k]] + j)];
1088: if (i >= k) { /* only take upper triangular entry */
1089: cols[ncols_upper] = i;
1090: cols_lvl[ncols_upper] = -1; /* initialize level for nonzero entries */
1091: ncols_upper++;
1092: }
1093: }
1094: PetscCall(PetscIncompleteLLAdd(ncols_upper, cols, levels, cols_lvl, am, &nlnk, lnk, lnk_lvl, lnkbt));
1095: nzk += nlnk;
1097: /* update lnk by computing fill-in for each pivot row to be merged in */
1098: prow = jl[k]; /* 1st pivot row */
1100: while (prow < k) {
1101: nextprow = jl[prow];
1103: /* merge prow into k-th row */
1104: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
1105: jmax = ui[prow + 1];
1106: ncols = jmax - jmin;
1107: i = jmin - ui[prow];
1108: cols = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
1109: for (j = 0; j < ncols; j++) cols_lvl[j] = *(uj_lvl_ptr[prow] + i + j);
1110: PetscCall(PetscIncompleteLLAddSorted(ncols, cols, levels, cols_lvl, am, &nlnk, lnk, lnk_lvl, lnkbt));
1111: nzk += nlnk;
1113: /* update il and jl for prow */
1114: if (jmin < jmax) {
1115: il[prow] = jmin;
1117: j = *cols;
1118: jl[prow] = jl[j];
1119: jl[j] = prow;
1120: }
1121: prow = nextprow;
1122: }
1124: /* if free space is not available, make more free space */
1125: if (current_space->local_remaining < nzk) {
1126: i = am - k + 1; /* num of unfactored rows */
1127: i = PetscMin(PetscIntMultTruncate(i, nzk), PetscIntMultTruncate(i, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1128: PetscCall(PetscFreeSpaceGet(i, ¤t_space));
1129: PetscCall(PetscFreeSpaceGet(i, ¤t_space_lvl));
1130: reallocs++;
1131: }
1133: /* copy data into free_space and free_space_lvl, then initialize lnk */
1134: PetscCall(PetscIncompleteLLClean(am, am, nzk, lnk, lnk_lvl, current_space->array, current_space_lvl->array, lnkbt));
1136: /* add the k-th row into il and jl */
1137: if (nzk - 1 > 0) {
1138: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
1139: jl[k] = jl[i];
1140: jl[i] = k;
1141: il[k] = ui[k] + 1;
1142: }
1143: uj_ptr[k] = current_space->array;
1144: uj_lvl_ptr[k] = current_space_lvl->array;
1146: current_space->array += nzk;
1147: current_space->local_used += nzk;
1148: current_space->local_remaining -= nzk;
1150: current_space_lvl->array += nzk;
1151: current_space_lvl->local_used += nzk;
1152: current_space_lvl->local_remaining -= nzk;
1154: ui[k + 1] = ui[k] + nzk;
1155: }
1157: PetscCall(ISRestoreIndices(perm, &rip));
1158: PetscCall(PetscFree4(uj_ptr, uj_lvl_ptr, il, jl));
1159: PetscCall(PetscFree(cols_lvl));
1161: /* copy free_space into uj and free free_space; set uj in new datastructure; */
1162: PetscCall(PetscMalloc1(ui[am] + 1, &uj));
1163: PetscCall(PetscFreeSpaceContiguous(&free_space, uj));
1164: PetscCall(PetscIncompleteLLDestroy(lnk, lnkbt));
1165: PetscCall(PetscFreeSpaceDestroy(free_space_lvl));
1167: /* put together the new matrix in MATSEQSBAIJ format */
1168: B = fact;
1169: PetscCall(MatSeqSBAIJSetPreallocation(B, 1, MAT_SKIP_ALLOCATION, NULL));
1171: b = (Mat_SeqSBAIJ *)B->data;
1172: b->free_a = PETSC_TRUE;
1173: b->free_ij = PETSC_TRUE;
1175: PetscCall(PetscMalloc1(ui[am] + 1, &b->a));
1177: b->j = uj;
1178: b->i = ui;
1179: b->diag = NULL;
1180: b->ilen = NULL;
1181: b->imax = NULL;
1182: b->row = perm;
1183: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
1185: PetscCall(PetscObjectReference((PetscObject)perm));
1187: b->icol = perm;
1189: PetscCall(PetscObjectReference((PetscObject)perm));
1190: PetscCall(PetscMalloc1(am + 1, &b->solve_work));
1192: b->maxnz = b->nz = ui[am];
1194: B->info.factor_mallocs = reallocs;
1195: B->info.fill_ratio_given = fill;
1196: if (ai[am] != 0.) {
1197: /* nonzeros in lower triangular part of A (including diagonals)= (ai[am]+am)/2 */
1198: B->info.fill_ratio_needed = ((PetscReal)2 * ui[am]) / (ai[am] + am);
1199: } else {
1200: B->info.fill_ratio_needed = 0.0;
1201: }
1202: #if defined(PETSC_USE_INFO)
1203: if (ai[am] != 0) {
1204: PetscReal af = B->info.fill_ratio_needed;
1205: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
1206: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
1207: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
1208: } else {
1209: PetscCall(PetscInfo(A, "Empty matrix\n"));
1210: }
1211: #endif
1212: if (perm_identity) {
1213: B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1214: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1215: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1216: } else {
1217: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
1218: }
1219: PetscFunctionReturn(PETSC_SUCCESS);
1220: }
1222: PetscErrorCode MatCholeskyFactorSymbolic_SeqBAIJ(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
1223: {
1224: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1225: Mat_SeqSBAIJ *b;
1226: Mat B;
1227: PetscBool perm_identity, missing;
1228: PetscReal fill = info->fill;
1229: const PetscInt *rip;
1230: PetscInt i, mbs = a->mbs, bs = A->rmap->bs, *ai = a->i, *aj = a->j, reallocs = 0, prow;
1231: PetscInt *jl, jmin, jmax, nzk, *ui, k, j, *il, nextprow;
1232: PetscInt nlnk, *lnk, ncols, ncols_upper, *cols, *uj, **ui_ptr, *uj_ptr;
1233: PetscFreeSpaceList free_space = NULL, current_space = NULL;
1234: PetscBT lnkbt;
1236: PetscFunctionBegin;
1237: if (bs > 1) { /* convert to seqsbaij */
1238: if (!a->sbaijMat) PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &a->sbaijMat));
1239: fact->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */
1241: PetscCall(MatCholeskyFactorSymbolic(fact, a->sbaijMat, perm, info));
1242: PetscFunctionReturn(PETSC_SUCCESS);
1243: }
1245: PetscCall(MatMissingDiagonal(A, &missing, &i));
1246: PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry %" PetscInt_FMT, i);
1248: /* check whether perm is the identity mapping */
1249: PetscCall(ISIdentity(perm, &perm_identity));
1250: PetscCheck(perm_identity, PETSC_COMM_SELF, PETSC_ERR_SUP, "Matrix reordering is not supported");
1251: PetscCall(ISGetIndices(perm, &rip));
1253: /* initialization */
1254: PetscCall(PetscMalloc1(mbs + 1, &ui));
1255: ui[0] = 0;
1257: /* jl: linked list for storing indices of the pivot rows
1258: il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
1259: PetscCall(PetscMalloc4(mbs, &ui_ptr, mbs, &il, mbs, &jl, mbs, &cols));
1260: for (i = 0; i < mbs; i++) {
1261: jl[i] = mbs;
1262: il[i] = 0;
1263: }
1265: /* create and initialize a linked list for storing column indices of the active row k */
1266: nlnk = mbs + 1;
1267: PetscCall(PetscLLCreate(mbs, mbs, nlnk, lnk, lnkbt));
1269: /* initial FreeSpace size is fill* (ai[mbs]+mbs)/2 */
1270: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[mbs] / 2, mbs / 2)), &free_space));
1272: current_space = free_space;
1274: for (k = 0; k < mbs; k++) { /* for each active row k */
1275: /* initialize lnk by the column indices of row rip[k] of A */
1276: nzk = 0;
1277: ncols = ai[rip[k] + 1] - ai[rip[k]];
1278: ncols_upper = 0;
1279: for (j = 0; j < ncols; j++) {
1280: i = rip[*(aj + ai[rip[k]] + j)];
1281: if (i >= k) { /* only take upper triangular entry */
1282: cols[ncols_upper] = i;
1283: ncols_upper++;
1284: }
1285: }
1286: PetscCall(PetscLLAdd(ncols_upper, cols, mbs, &nlnk, lnk, lnkbt));
1287: nzk += nlnk;
1289: /* update lnk by computing fill-in for each pivot row to be merged in */
1290: prow = jl[k]; /* 1st pivot row */
1292: while (prow < k) {
1293: nextprow = jl[prow];
1294: /* merge prow into k-th row */
1295: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
1296: jmax = ui[prow + 1];
1297: ncols = jmax - jmin;
1298: uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
1299: PetscCall(PetscLLAddSorted(ncols, uj_ptr, mbs, &nlnk, lnk, lnkbt));
1300: nzk += nlnk;
1302: /* update il and jl for prow */
1303: if (jmin < jmax) {
1304: il[prow] = jmin;
1305: j = *uj_ptr;
1306: jl[prow] = jl[j];
1307: jl[j] = prow;
1308: }
1309: prow = nextprow;
1310: }
1312: /* if free space is not available, make more free space */
1313: if (current_space->local_remaining < nzk) {
1314: i = mbs - k + 1; /* num of unfactored rows */
1315: i = PetscMin(PetscIntMultTruncate(i, nzk), PetscIntMultTruncate(i, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1316: PetscCall(PetscFreeSpaceGet(i, ¤t_space));
1317: reallocs++;
1318: }
1320: /* copy data into free space, then initialize lnk */
1321: PetscCall(PetscLLClean(mbs, mbs, nzk, lnk, current_space->array, lnkbt));
1323: /* add the k-th row into il and jl */
1324: if (nzk - 1 > 0) {
1325: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
1326: jl[k] = jl[i];
1327: jl[i] = k;
1328: il[k] = ui[k] + 1;
1329: }
1330: ui_ptr[k] = current_space->array;
1331: current_space->array += nzk;
1332: current_space->local_used += nzk;
1333: current_space->local_remaining -= nzk;
1335: ui[k + 1] = ui[k] + nzk;
1336: }
1338: PetscCall(ISRestoreIndices(perm, &rip));
1339: PetscCall(PetscFree4(ui_ptr, il, jl, cols));
1341: /* copy free_space into uj and free free_space; set uj in new datastructure; */
1342: PetscCall(PetscMalloc1(ui[mbs] + 1, &uj));
1343: PetscCall(PetscFreeSpaceContiguous(&free_space, uj));
1344: PetscCall(PetscLLDestroy(lnk, lnkbt));
1346: /* put together the new matrix in MATSEQSBAIJ format */
1347: B = fact;
1348: PetscCall(MatSeqSBAIJSetPreallocation(B, bs, MAT_SKIP_ALLOCATION, NULL));
1350: b = (Mat_SeqSBAIJ *)B->data;
1351: b->free_a = PETSC_TRUE;
1352: b->free_ij = PETSC_TRUE;
1354: PetscCall(PetscMalloc1(ui[mbs] + 1, &b->a));
1356: b->j = uj;
1357: b->i = ui;
1358: b->diag = NULL;
1359: b->ilen = NULL;
1360: b->imax = NULL;
1361: b->row = perm;
1362: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
1364: PetscCall(PetscObjectReference((PetscObject)perm));
1365: b->icol = perm;
1366: PetscCall(PetscObjectReference((PetscObject)perm));
1367: PetscCall(PetscMalloc1(mbs + 1, &b->solve_work));
1368: b->maxnz = b->nz = ui[mbs];
1370: B->info.factor_mallocs = reallocs;
1371: B->info.fill_ratio_given = fill;
1372: if (ai[mbs] != 0.) {
1373: /* nonzeros in lower triangular part of A = (ai[mbs]+mbs)/2 */
1374: B->info.fill_ratio_needed = ((PetscReal)2 * ui[mbs]) / (ai[mbs] + mbs);
1375: } else {
1376: B->info.fill_ratio_needed = 0.0;
1377: }
1378: #if defined(PETSC_USE_INFO)
1379: if (ai[mbs] != 0.) {
1380: PetscReal af = B->info.fill_ratio_needed;
1381: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
1382: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
1383: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
1384: } else {
1385: PetscCall(PetscInfo(A, "Empty matrix\n"));
1386: }
1387: #endif
1388: if (perm_identity) {
1389: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1390: } else {
1391: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
1392: }
1393: PetscFunctionReturn(PETSC_SUCCESS);
1394: }
1396: PetscErrorCode MatSolve_SeqBAIJ_N_NaturalOrdering(Mat A, Vec bb, Vec xx)
1397: {
1398: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1399: const PetscInt *ai = a->i, *aj = a->j, *adiag = a->diag, *vi;
1400: PetscInt i, k, n = a->mbs;
1401: PetscInt nz, bs = A->rmap->bs, bs2 = a->bs2;
1402: const MatScalar *aa = a->a, *v;
1403: PetscScalar *x, *s, *t, *ls;
1404: const PetscScalar *b;
1406: PetscFunctionBegin;
1407: PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected bs %" PetscInt_FMT " > 0", bs);
1408: PetscCall(VecGetArrayRead(bb, &b));
1409: PetscCall(VecGetArray(xx, &x));
1410: t = a->solve_work;
1412: /* forward solve the lower triangular */
1413: PetscCall(PetscArraycpy(t, b, bs)); /* copy 1st block of b to t */
1415: for (i = 1; i < n; i++) {
1416: v = aa + bs2 * ai[i];
1417: vi = aj + ai[i];
1418: nz = ai[i + 1] - ai[i];
1419: s = t + bs * i;
1420: PetscCall(PetscArraycpy(s, b + bs * i, bs)); /* copy i_th block of b to t */
1421: for (k = 0; k < nz; k++) {
1422: PetscKernel_v_gets_v_minus_A_times_w(bs, s, v, t + bs * vi[k]);
1423: v += bs2;
1424: }
1425: }
1427: /* backward solve the upper triangular */
1428: ls = a->solve_work + A->cmap->n;
1429: for (i = n - 1; i >= 0; i--) {
1430: v = aa + bs2 * (adiag[i + 1] + 1);
1431: vi = aj + adiag[i + 1] + 1;
1432: nz = adiag[i] - adiag[i + 1] - 1;
1433: PetscCall(PetscArraycpy(ls, t + i * bs, bs));
1434: for (k = 0; k < nz; k++) {
1435: PetscKernel_v_gets_v_minus_A_times_w(bs, ls, v, t + bs * vi[k]);
1436: v += bs2;
1437: }
1438: PetscKernel_w_gets_A_times_v(bs, ls, aa + bs2 * adiag[i], t + i * bs); /* *inv(diagonal[i]) */
1439: PetscCall(PetscArraycpy(x + i * bs, t + i * bs, bs));
1440: }
1442: PetscCall(VecRestoreArrayRead(bb, &b));
1443: PetscCall(VecRestoreArray(xx, &x));
1444: PetscCall(PetscLogFlops(2.0 * (a->bs2) * (a->nz) - A->rmap->bs * A->cmap->n));
1445: PetscFunctionReturn(PETSC_SUCCESS);
1446: }
1448: PetscErrorCode MatSolve_SeqBAIJ_N(Mat A, Vec bb, Vec xx)
1449: {
1450: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1451: IS iscol = a->col, isrow = a->row;
1452: const PetscInt *r, *c, *rout, *cout, *ai = a->i, *aj = a->j, *adiag = a->diag, *vi;
1453: PetscInt i, m, n = a->mbs;
1454: PetscInt nz, bs = A->rmap->bs, bs2 = a->bs2;
1455: const MatScalar *aa = a->a, *v;
1456: PetscScalar *x, *s, *t, *ls;
1457: const PetscScalar *b;
1459: PetscFunctionBegin;
1460: PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected bs %" PetscInt_FMT " > 0", bs);
1461: PetscCall(VecGetArrayRead(bb, &b));
1462: PetscCall(VecGetArray(xx, &x));
1463: t = a->solve_work;
1465: PetscCall(ISGetIndices(isrow, &rout));
1466: r = rout;
1467: PetscCall(ISGetIndices(iscol, &cout));
1468: c = cout;
1470: /* forward solve the lower triangular */
1471: PetscCall(PetscArraycpy(t, b + bs * r[0], bs));
1472: for (i = 1; i < n; i++) {
1473: v = aa + bs2 * ai[i];
1474: vi = aj + ai[i];
1475: nz = ai[i + 1] - ai[i];
1476: s = t + bs * i;
1477: PetscCall(PetscArraycpy(s, b + bs * r[i], bs));
1478: for (m = 0; m < nz; m++) {
1479: PetscKernel_v_gets_v_minus_A_times_w(bs, s, v, t + bs * vi[m]);
1480: v += bs2;
1481: }
1482: }
1484: /* backward solve the upper triangular */
1485: ls = a->solve_work + A->cmap->n;
1486: for (i = n - 1; i >= 0; i--) {
1487: v = aa + bs2 * (adiag[i + 1] + 1);
1488: vi = aj + adiag[i + 1] + 1;
1489: nz = adiag[i] - adiag[i + 1] - 1;
1490: PetscCall(PetscArraycpy(ls, t + i * bs, bs));
1491: for (m = 0; m < nz; m++) {
1492: PetscKernel_v_gets_v_minus_A_times_w(bs, ls, v, t + bs * vi[m]);
1493: v += bs2;
1494: }
1495: PetscKernel_w_gets_A_times_v(bs, ls, v, t + i * bs); /* *inv(diagonal[i]) */
1496: PetscCall(PetscArraycpy(x + bs * c[i], t + i * bs, bs));
1497: }
1498: PetscCall(ISRestoreIndices(isrow, &rout));
1499: PetscCall(ISRestoreIndices(iscol, &cout));
1500: PetscCall(VecRestoreArrayRead(bb, &b));
1501: PetscCall(VecRestoreArray(xx, &x));
1502: PetscCall(PetscLogFlops(2.0 * (a->bs2) * (a->nz) - A->rmap->bs * A->cmap->n));
1503: PetscFunctionReturn(PETSC_SUCCESS);
1504: }
1506: /*
1507: For each block in an block array saves the largest absolute value in the block into another array
1508: */
1509: static PetscErrorCode MatBlockAbs_private(PetscInt nbs, PetscInt bs2, PetscScalar *blockarray, PetscReal *absarray)
1510: {
1511: PetscInt i, j;
1513: PetscFunctionBegin;
1514: PetscCall(PetscArrayzero(absarray, nbs + 1));
1515: for (i = 0; i < nbs; i++) {
1516: for (j = 0; j < bs2; j++) {
1517: if (absarray[i] < PetscAbsScalar(blockarray[i * nbs + j])) absarray[i] = PetscAbsScalar(blockarray[i * nbs + j]);
1518: }
1519: }
1520: PetscFunctionReturn(PETSC_SUCCESS);
1521: }
1523: /*
1524: This needs to be renamed and called by the regular MatILUFactor_SeqBAIJ when drop tolerance is used
1525: */
1526: PetscErrorCode MatILUDTFactor_SeqBAIJ(Mat A, IS isrow, IS iscol, const MatFactorInfo *info, Mat *fact)
1527: {
1528: Mat B = *fact;
1529: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b;
1530: IS isicol;
1531: const PetscInt *r, *ic;
1532: PetscInt i, mbs = a->mbs, bs = A->rmap->bs, bs2 = a->bs2, *ai = a->i, *aj = a->j, *ajtmp, *adiag;
1533: PetscInt *bi, *bj, *bdiag;
1535: PetscInt row, nzi, nzi_bl, nzi_bu, *im, dtcount, nzi_al, nzi_au;
1536: PetscInt nlnk, *lnk;
1537: PetscBT lnkbt;
1538: PetscBool row_identity, icol_identity;
1539: MatScalar *aatmp, *pv, *batmp, *ba, *rtmp, *pc, *multiplier, *vtmp;
1540: PetscInt j, nz, *pj, *bjtmp, k, ncut, *jtmp;
1542: PetscReal dt = info->dt; /* shift=info->shiftamount; */
1543: PetscInt nnz_max;
1544: PetscBool missing;
1545: PetscReal *vtmp_abs;
1546: MatScalar *v_work;
1547: PetscInt *v_pivots;
1548: PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE;
1550: PetscFunctionBegin;
1551: /* symbolic factorization, can be reused */
1552: PetscCall(MatMissingDiagonal(A, &missing, &i));
1553: PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry %" PetscInt_FMT, i);
1554: adiag = a->diag;
1556: PetscCall(ISInvertPermutation(iscol, PETSC_DECIDE, &isicol));
1558: /* bdiag is location of diagonal in factor */
1559: PetscCall(PetscMalloc1(mbs + 1, &bdiag));
1561: /* allocate row pointers bi */
1562: PetscCall(PetscMalloc1(2 * mbs + 2, &bi));
1564: /* allocate bj and ba; max num of nonzero entries is (ai[n]+2*n*dtcount+2) */
1565: dtcount = (PetscInt)info->dtcount;
1566: if (dtcount > mbs - 1) dtcount = mbs - 1;
1567: nnz_max = ai[mbs] + 2 * mbs * dtcount + 2;
1568: /* printf("MatILUDTFactor_SeqBAIJ, bs %d, ai[mbs] %d, nnz_max %d, dtcount %d\n",bs,ai[mbs],nnz_max,dtcount); */
1569: PetscCall(PetscMalloc1(nnz_max, &bj));
1570: nnz_max = nnz_max * bs2;
1571: PetscCall(PetscMalloc1(nnz_max, &ba));
1573: /* put together the new matrix */
1574: PetscCall(MatSeqBAIJSetPreallocation(B, bs, MAT_SKIP_ALLOCATION, NULL));
1576: b = (Mat_SeqBAIJ *)B->data;
1577: b->free_a = PETSC_TRUE;
1578: b->free_ij = PETSC_TRUE;
1580: b->a = ba;
1581: b->j = bj;
1582: b->i = bi;
1583: b->diag = bdiag;
1584: b->ilen = NULL;
1585: b->imax = NULL;
1586: b->row = isrow;
1587: b->col = iscol;
1589: PetscCall(PetscObjectReference((PetscObject)isrow));
1590: PetscCall(PetscObjectReference((PetscObject)iscol));
1592: b->icol = isicol;
1593: PetscCall(PetscMalloc1(bs * (mbs + 1), &b->solve_work));
1594: b->maxnz = nnz_max / bs2;
1596: B->factortype = MAT_FACTOR_ILUDT;
1597: B->info.factor_mallocs = 0;
1598: B->info.fill_ratio_given = ((PetscReal)nnz_max) / ((PetscReal)(ai[mbs] * bs2));
1599: /* end of symbolic factorization */
1600: PetscCall(ISGetIndices(isrow, &r));
1601: PetscCall(ISGetIndices(isicol, &ic));
1603: /* linked list for storing column indices of the active row */
1604: nlnk = mbs + 1;
1605: PetscCall(PetscLLCreate(mbs, mbs, nlnk, lnk, lnkbt));
1607: /* im: used by PetscLLAddSortedLU(); jtmp: working array for column indices of active row */
1608: PetscCall(PetscMalloc2(mbs, &im, mbs, &jtmp));
1609: /* rtmp, vtmp: working arrays for sparse and contiguous row entries of active row */
1610: PetscCall(PetscMalloc2(mbs * bs2, &rtmp, mbs * bs2, &vtmp));
1611: PetscCall(PetscMalloc1(mbs + 1, &vtmp_abs));
1612: PetscCall(PetscMalloc3(bs, &v_work, bs2, &multiplier, bs, &v_pivots));
1614: allowzeropivot = PetscNot(A->erroriffailure);
1615: bi[0] = 0;
1616: bdiag[0] = (nnz_max / bs2) - 1; /* location of diagonal in factor B */
1617: bi[2 * mbs + 1] = bdiag[0] + 1; /* endof bj and ba array */
1618: for (i = 0; i < mbs; i++) {
1619: /* copy initial fill into linked list */
1620: nzi = ai[r[i] + 1] - ai[r[i]];
1621: PetscCheck(nzi, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Empty row in matrix: row in original ordering %" PetscInt_FMT " in permuted ordering %" PetscInt_FMT, r[i], i);
1622: nzi_al = adiag[r[i]] - ai[r[i]];
1623: nzi_au = ai[r[i] + 1] - adiag[r[i]] - 1;
1625: /* load in initial unfactored row */
1626: ajtmp = aj + ai[r[i]];
1627: PetscCall(PetscLLAddPerm(nzi, ajtmp, ic, mbs, &nlnk, lnk, lnkbt));
1628: PetscCall(PetscArrayzero(rtmp, mbs * bs2));
1629: aatmp = a->a + bs2 * ai[r[i]];
1630: for (j = 0; j < nzi; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ic[ajtmp[j]], aatmp + bs2 * j, bs2));
1632: /* add pivot rows into linked list */
1633: row = lnk[mbs];
1634: while (row < i) {
1635: nzi_bl = bi[row + 1] - bi[row] + 1;
1636: bjtmp = bj + bdiag[row + 1] + 1; /* points to 1st column next to the diagonal in U */
1637: PetscCall(PetscLLAddSortedLU(bjtmp, row, &nlnk, lnk, lnkbt, i, nzi_bl, im));
1638: nzi += nlnk;
1639: row = lnk[row];
1640: }
1642: /* copy data from lnk into jtmp, then initialize lnk */
1643: PetscCall(PetscLLClean(mbs, mbs, nzi, lnk, jtmp, lnkbt));
1645: /* numerical factorization */
1646: bjtmp = jtmp;
1647: row = *bjtmp++; /* 1st pivot row */
1649: while (row < i) {
1650: pc = rtmp + bs2 * row;
1651: pv = ba + bs2 * bdiag[row]; /* inv(diag) of the pivot row */
1652: PetscKernel_A_gets_A_times_B(bs, pc, pv, multiplier); /* pc= multiplier = pc*inv(diag[row]) */
1653: PetscCall(MatBlockAbs_private(1, bs2, pc, vtmp_abs));
1654: if (vtmp_abs[0] > dt) { /* apply tolerance dropping rule */
1655: pj = bj + bdiag[row + 1] + 1; /* point to 1st entry of U(row,:) */
1656: pv = ba + bs2 * (bdiag[row + 1] + 1);
1657: nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries in U(row,:), excluding diagonal */
1658: for (j = 0; j < nz; j++) PetscKernel_A_gets_A_minus_B_times_C(bs, rtmp + bs2 * pj[j], pc, pv + bs2 * j);
1659: /* PetscCall(PetscLogFlops(bslog*(nz+1.0)-bs)); */
1660: }
1661: row = *bjtmp++;
1662: }
1664: /* copy sparse rtmp into contiguous vtmp; separate L and U part */
1665: nzi_bl = 0;
1666: j = 0;
1667: while (jtmp[j] < i) { /* L-part. Note: jtmp is sorted */
1668: PetscCall(PetscArraycpy(vtmp + bs2 * j, rtmp + bs2 * jtmp[j], bs2));
1669: nzi_bl++;
1670: j++;
1671: }
1672: nzi_bu = nzi - nzi_bl - 1;
1674: while (j < nzi) { /* U-part */
1675: PetscCall(PetscArraycpy(vtmp + bs2 * j, rtmp + bs2 * jtmp[j], bs2));
1676: j++;
1677: }
1679: PetscCall(MatBlockAbs_private(nzi, bs2, vtmp, vtmp_abs));
1681: bjtmp = bj + bi[i];
1682: batmp = ba + bs2 * bi[i];
1683: /* apply level dropping rule to L part */
1684: ncut = nzi_al + dtcount;
1685: if (ncut < nzi_bl) {
1686: PetscCall(PetscSortSplitReal(ncut, nzi_bl, vtmp_abs, jtmp));
1687: PetscCall(PetscSortIntWithScalarArray(ncut, jtmp, vtmp));
1688: } else {
1689: ncut = nzi_bl;
1690: }
1691: for (j = 0; j < ncut; j++) {
1692: bjtmp[j] = jtmp[j];
1693: PetscCall(PetscArraycpy(batmp + bs2 * j, rtmp + bs2 * bjtmp[j], bs2));
1694: }
1695: bi[i + 1] = bi[i] + ncut;
1696: nzi = ncut + 1;
1698: /* apply level dropping rule to U part */
1699: ncut = nzi_au + dtcount;
1700: if (ncut < nzi_bu) {
1701: PetscCall(PetscSortSplitReal(ncut, nzi_bu, vtmp_abs + nzi_bl + 1, jtmp + nzi_bl + 1));
1702: PetscCall(PetscSortIntWithScalarArray(ncut, jtmp + nzi_bl + 1, vtmp + nzi_bl + 1));
1703: } else {
1704: ncut = nzi_bu;
1705: }
1706: nzi += ncut;
1708: /* mark bdiagonal */
1709: bdiag[i + 1] = bdiag[i] - (ncut + 1);
1710: bi[2 * mbs - i] = bi[2 * mbs - i + 1] - (ncut + 1);
1712: bjtmp = bj + bdiag[i];
1713: batmp = ba + bs2 * bdiag[i];
1714: PetscCall(PetscArraycpy(batmp, rtmp + bs2 * i, bs2));
1715: *bjtmp = i;
1717: bjtmp = bj + bdiag[i + 1] + 1;
1718: batmp = ba + (bdiag[i + 1] + 1) * bs2;
1720: for (k = 0; k < ncut; k++) {
1721: bjtmp[k] = jtmp[nzi_bl + 1 + k];
1722: PetscCall(PetscArraycpy(batmp + bs2 * k, rtmp + bs2 * bjtmp[k], bs2));
1723: }
1725: im[i] = nzi; /* used by PetscLLAddSortedLU() */
1727: /* invert diagonal block for simpler triangular solves - add shift??? */
1728: batmp = ba + bs2 * bdiag[i];
1730: PetscCall(PetscKernel_A_gets_inverse_A(bs, batmp, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
1731: if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1732: } /* for (i=0; i<mbs; i++) */
1733: PetscCall(PetscFree3(v_work, multiplier, v_pivots));
1735: /* printf("end of L %d, beginning of U %d\n",bi[mbs],bdiag[mbs]); */
1736: PetscCheck(bi[mbs] < bdiag[mbs], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "end of L array %" PetscInt_FMT " cannot >= the beginning of U array %" PetscInt_FMT, bi[mbs], bdiag[mbs]);
1738: PetscCall(ISRestoreIndices(isrow, &r));
1739: PetscCall(ISRestoreIndices(isicol, &ic));
1741: PetscCall(PetscLLDestroy(lnk, lnkbt));
1743: PetscCall(PetscFree2(im, jtmp));
1744: PetscCall(PetscFree2(rtmp, vtmp));
1746: PetscCall(PetscLogFlops(bs2 * B->cmap->n));
1747: b->maxnz = b->nz = bi[mbs] + bdiag[0] - bdiag[mbs];
1749: PetscCall(ISIdentity(isrow, &row_identity));
1750: PetscCall(ISIdentity(isicol, &icol_identity));
1751: if (row_identity && icol_identity) {
1752: B->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
1753: } else {
1754: B->ops->solve = MatSolve_SeqBAIJ_N;
1755: }
1757: B->ops->solveadd = NULL;
1758: B->ops->solvetranspose = NULL;
1759: B->ops->solvetransposeadd = NULL;
1760: B->ops->matsolve = NULL;
1761: B->assembled = PETSC_TRUE;
1762: B->preallocated = PETSC_TRUE;
1763: PetscFunctionReturn(PETSC_SUCCESS);
1764: }