Actual source code: aijfact.c
1: #include <../src/mat/impls/aij/seq/aij.h>
2: #include <../src/mat/impls/sbaij/seq/sbaij.h>
3: #include <../src/mat/utils/freespace.h>
5: static PetscErrorCode MatFactorGetSolverType_petsc(Mat A, MatSolverType *type)
6: {
7: PetscFunctionBegin;
8: *type = MATSOLVERPETSC;
9: PetscFunctionReturn(PETSC_SUCCESS);
10: }
12: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_petsc(Mat A, MatFactorType ftype, Mat *B)
13: {
14: PetscInt n = A->rmap->n;
16: PetscFunctionBegin;
17: if (PetscDefined(USE_COMPLEX) && (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) && A->hermitian == PETSC_BOOL3_TRUE && A->symmetric != PETSC_BOOL3_TRUE) {
18: PetscCall(PetscInfo(A, "Hermitian MAT_FACTOR_CHOLESKY or MAT_FACTOR_ICC are not supported. Use MAT_FACTOR_LU instead.\n"));
19: *B = NULL;
20: PetscFunctionReturn(PETSC_SUCCESS);
21: }
23: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
24: PetscCall(MatSetSizes(*B, n, n, n, n));
25: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU || ftype == MAT_FACTOR_ILUDT) {
26: PetscCall(MatSetType(*B, MATSEQAIJ));
28: (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqAIJ;
29: (*B)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ;
31: PetscCall(MatSetBlockSizesFromMats(*B, A, A));
32: PetscCall(PetscStrallocpy(MATORDERINGND, (char **)&(*B)->preferredordering[MAT_FACTOR_LU]));
33: PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ILU]));
34: PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ILUDT]));
35: } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
36: PetscCall(MatSetType(*B, MATSEQSBAIJ));
37: PetscCall(MatSeqSBAIJSetPreallocation(*B, 1, MAT_SKIP_ALLOCATION, NULL));
39: (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqAIJ;
40: (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqAIJ;
41: PetscCall(PetscStrallocpy(MATORDERINGND, (char **)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]));
42: PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ICC]));
43: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");
44: (*B)->factortype = ftype;
46: PetscCall(PetscFree((*B)->solvertype));
47: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
48: (*B)->canuseordering = PETSC_TRUE;
49: PetscCall(PetscObjectComposeFunction((PetscObject)*B, "MatFactorGetSolverType_C", MatFactorGetSolverType_petsc));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: PetscErrorCode MatLUFactorSymbolic_SeqAIJ(Mat B, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
54: {
55: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b;
56: IS isicol;
57: const PetscInt *r, *ic, *ai = a->i, *aj = a->j, *ajtmp;
58: PetscInt i, n = A->rmap->n;
59: PetscInt *bi, *bj;
60: PetscInt *bdiag, row, nnz, nzi, reallocs = 0, nzbd, *im;
61: PetscReal f;
62: PetscInt nlnk, *lnk, k, **bi_ptr;
63: PetscFreeSpaceList free_space = NULL, current_space = NULL;
64: PetscBT lnkbt;
65: PetscBool diagDense;
67: PetscFunctionBegin;
68: PetscCheck(A->rmap->N == A->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "matrix must be square");
69: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, NULL, &diagDense));
70: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entries");
72: PetscCall(ISInvertPermutation(iscol, PETSC_DECIDE, &isicol));
73: PetscCall(ISGetIndices(isrow, &r));
74: PetscCall(ISGetIndices(isicol, &ic));
76: /* get new row and diagonal pointers, must be allocated separately because they will be given to the Mat_SeqAIJ and freed separately */
77: PetscCall(PetscShmgetAllocateArray(n + 1, sizeof(PetscInt), (void **)&bi));
78: PetscCall(PetscMalloc1(n + 1, &bdiag));
79: bi[0] = bdiag[0] = 0;
81: /* linked list for storing column indices of the active row */
82: nlnk = n + 1;
83: PetscCall(PetscLLCreate(n, n, nlnk, lnk, lnkbt));
85: PetscCall(PetscMalloc2(n + 1, &bi_ptr, n + 1, &im));
87: /* initial FreeSpace size is f*(ai[n]+1) */
88: f = info->fill;
89: if (n == 1) f = 1; /* prevent failure in corner case of 1x1 matrix with fill < 0.5 */
90: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space));
91: current_space = free_space;
93: for (i = 0; i < n; i++) {
94: /* copy previous fill into linked list */
95: nzi = 0;
96: nnz = ai[r[i] + 1] - ai[r[i]];
97: ajtmp = aj + ai[r[i]];
98: PetscCall(PetscLLAddPerm(nnz, ajtmp, ic, n, &nlnk, lnk, lnkbt));
99: nzi += nlnk;
101: /* add pivot rows into linked list */
102: row = lnk[n];
103: while (row < i) {
104: nzbd = bdiag[row] + 1; /* num of entries in the row with column index <= row */
105: ajtmp = bi_ptr[row] + nzbd; /* points to the entry next to the diagonal */
106: PetscCall(PetscLLAddSortedLU(ajtmp, row, &nlnk, lnk, lnkbt, i, nzbd, im));
107: nzi += nlnk;
108: row = lnk[row];
109: }
110: bi[i + 1] = bi[i] + nzi;
111: im[i] = nzi;
113: /* mark bdiag */
114: nzbd = 0;
115: nnz = nzi;
116: k = lnk[n];
117: while (nnz-- && k < i) {
118: nzbd++;
119: k = lnk[k];
120: }
121: bdiag[i] = nzbd; /* note: bdiag[i] = nnzL as input for PetscFreeSpaceContiguous_LU() */
123: /* if free space is not available, make more free space */
124: if (current_space->local_remaining < nzi) {
125: /* estimated additional space needed */
126: nnz = PetscIntMultTruncate(2, PetscIntMultTruncate(n - 1, nzi));
127: PetscCall(PetscFreeSpaceGet(nnz, ¤t_space));
128: reallocs++;
129: }
131: /* copy data into free space, then initialize lnk */
132: PetscCall(PetscLLClean(n, n, nzi, lnk, current_space->array, lnkbt));
134: bi_ptr[i] = current_space->array;
135: current_space->array += nzi;
136: current_space->local_used += nzi;
137: current_space->local_remaining -= nzi;
138: }
140: PetscCall(ISRestoreIndices(isrow, &r));
141: PetscCall(ISRestoreIndices(isicol, &ic));
143: /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
144: PetscCall(PetscShmgetAllocateArray(bi[n], sizeof(PetscInt), (void **)&bj));
145: PetscCall(PetscFreeSpaceContiguous_LU(&free_space, bj, n, bi, bdiag));
146: PetscCall(PetscLLDestroy(lnk, lnkbt));
147: PetscCall(PetscFree2(bi_ptr, im));
149: /* put together the new matrix */
150: PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(B, MAT_SKIP_ALLOCATION, NULL));
151: b = (Mat_SeqAIJ *)B->data;
152: b->free_ij = PETSC_TRUE;
153: PetscCall(PetscShmgetAllocateArray(bdiag[0] + 1, sizeof(PetscScalar), (void **)&b->a));
154: b->free_a = PETSC_TRUE;
155: b->j = bj;
156: b->i = bi;
157: b->diag = bdiag;
158: b->ilen = NULL;
159: b->imax = NULL;
160: b->row = isrow;
161: b->col = iscol;
162: PetscCall(PetscObjectReference((PetscObject)isrow));
163: PetscCall(PetscObjectReference((PetscObject)iscol));
164: b->icol = isicol;
165: PetscCall(PetscMalloc1(n, &b->solve_work));
167: /* In b structure: Free imax, ilen, old a, old j. Allocate solve_work, new a, new j */
168: b->maxnz = b->nz = bdiag[0] + 1;
170: B->factortype = MAT_FACTOR_LU;
171: B->info.factor_mallocs = reallocs;
172: B->info.fill_ratio_given = f;
174: if (ai[n]) {
175: B->info.fill_ratio_needed = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);
176: } else {
177: B->info.fill_ratio_needed = 0.0;
178: }
179: #if PetscDefined(USE_INFO)
180: if (ai[n] != 0) {
181: PetscReal af = B->info.fill_ratio_needed;
182: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)f, (double)af));
183: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
184: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g);\n", (double)af));
185: PetscCall(PetscInfo(A, "for best performance.\n"));
186: } else PetscCall(PetscInfo(A, "Empty matrix\n"));
187: #endif
188: B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ;
189: if (a->inode.size_csr) B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
190: PetscCall(MatSeqAIJCheckInode_FactorLU(B));
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: /*
195: Trouble in factorization, should we dump the original matrix?
196: */
197: PetscErrorCode MatFactorDumpMatrix(Mat A)
198: {
199: PetscBool flg = PETSC_FALSE;
201: PetscFunctionBegin;
202: PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, NULL, "-mat_factor_dump_on_error", &flg, NULL));
203: if (flg) {
204: PetscViewer viewer;
205: char filename[PETSC_MAX_PATH_LEN];
207: PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN, "matrix_factor_error.%d", PetscGlobalRank));
208: PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)A), filename, FILE_MODE_WRITE, &viewer));
209: PetscCall(MatView(A, viewer));
210: PetscCall(PetscViewerDestroy(&viewer));
211: }
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: PetscErrorCode MatLUFactorNumeric_SeqAIJ(Mat B, Mat A, const MatFactorInfo *info)
216: {
217: Mat C = B;
218: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)C->data;
219: IS isrow = b->row, isicol = b->icol;
220: const PetscInt *r, *ic, *ics;
221: const PetscInt n = A->rmap->n, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *bdiag = b->diag;
222: PetscInt i, j, k, nz, nzL, row, *pj;
223: const PetscInt *ajtmp, *bjtmp;
224: MatScalar *rtmp, *pc, multiplier, *pv;
225: const MatScalar *aa, *v;
226: MatScalar *ba;
227: PetscBool row_identity, col_identity;
228: FactorShiftCtx sctx;
229: const PetscInt *ddiag;
230: PetscReal rs;
231: MatScalar d;
233: PetscFunctionBegin;
234: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
235: PetscCall(MatSeqAIJGetArrayWrite(B, &ba));
236: /* MatPivotSetUp(): initialize shift context sctx */
237: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
239: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
240: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &ddiag, NULL));
241: sctx.shift_top = info->zeropivot;
242: for (i = 0; i < n; i++) {
243: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
244: d = aa[ddiag[i]];
245: rs = -PetscAbsScalar(d) - PetscRealPart(d);
246: v = aa + ai[i];
247: nz = ai[i + 1] - ai[i];
248: for (j = 0; j < nz; j++) rs += PetscAbsScalar(v[j]);
249: if (rs > sctx.shift_top) sctx.shift_top = rs;
250: }
251: sctx.shift_top *= 1.1;
252: sctx.nshift_max = 5;
253: sctx.shift_lo = 0.;
254: sctx.shift_hi = 1.;
255: }
257: PetscCall(ISGetIndices(isrow, &r));
258: PetscCall(ISGetIndices(isicol, &ic));
259: PetscCall(PetscMalloc1(n + 1, &rtmp));
260: ics = ic;
262: do {
263: sctx.newshift = PETSC_FALSE;
264: for (i = 0; i < n; i++) {
265: /* zero rtmp */
266: /* L part */
267: nz = bi[i + 1] - bi[i];
268: bjtmp = bj + bi[i];
269: for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;
271: /* U part */
272: nz = bdiag[i] - bdiag[i + 1];
273: bjtmp = bj + bdiag[i + 1] + 1;
274: for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;
276: /* load in initial (unfactored row) */
277: nz = ai[r[i] + 1] - ai[r[i]];
278: ajtmp = aj + ai[r[i]];
279: v = aa + ai[r[i]];
280: for (j = 0; j < nz; j++) rtmp[ics[ajtmp[j]]] = v[j];
281: /* ZeropivotApply() */
282: rtmp[i] += sctx.shift_amount; /* shift the diagonal of the matrix */
284: /* elimination */
285: bjtmp = bj + bi[i];
286: row = *bjtmp++;
287: nzL = bi[i + 1] - bi[i];
288: for (k = 0; k < nzL; k++) {
289: pc = rtmp + row;
290: if (*pc != 0.0) {
291: pv = ba + bdiag[row];
292: multiplier = *pc * (*pv);
293: *pc = multiplier;
295: pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
296: pv = ba + bdiag[row + 1] + 1;
297: nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries in U(row,:) excluding diag */
299: for (j = 0; j < nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
300: PetscCall(PetscLogFlops(1 + 2.0 * nz));
301: }
302: row = *bjtmp++;
303: }
305: /* finished row so stick it into b->a */
306: rs = 0.0;
307: /* L part */
308: pv = ba + bi[i];
309: pj = b->j + bi[i];
310: nz = bi[i + 1] - bi[i];
311: for (j = 0; j < nz; j++) {
312: pv[j] = rtmp[pj[j]];
313: rs += PetscAbsScalar(pv[j]);
314: }
316: /* U part */
317: pv = ba + bdiag[i + 1] + 1;
318: pj = b->j + bdiag[i + 1] + 1;
319: nz = bdiag[i] - bdiag[i + 1] - 1;
320: for (j = 0; j < nz; j++) {
321: pv[j] = rtmp[pj[j]];
322: rs += PetscAbsScalar(pv[j]);
323: }
325: sctx.rs = rs;
326: sctx.pv = rtmp[i];
327: PetscCall(MatPivotCheck(B, A, info, &sctx, i));
328: if (sctx.newshift) break; /* break for-loop */
329: rtmp[i] = sctx.pv; /* sctx.pv might be updated in the case of MAT_SHIFT_INBLOCKS */
331: /* Mark diagonal and invert diagonal for simpler triangular solves */
332: pv = ba + bdiag[i];
333: *pv = 1.0 / rtmp[i];
335: } /* endof for (i=0; i<n; i++) { */
337: /* MatPivotRefine() */
338: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction > 0 && sctx.nshift < sctx.nshift_max) {
339: /*
340: * if no shift in this attempt & shifting & started shifting & can refine,
341: * then try lower shift
342: */
343: sctx.shift_hi = sctx.shift_fraction;
344: sctx.shift_fraction = (sctx.shift_hi + sctx.shift_lo) / 2.;
345: sctx.shift_amount = sctx.shift_fraction * sctx.shift_top;
346: sctx.newshift = PETSC_TRUE;
347: sctx.nshift++;
348: }
349: } while (sctx.newshift);
351: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
352: PetscCall(MatSeqAIJRestoreArrayWrite(B, &ba));
354: PetscCall(PetscFree(rtmp));
355: PetscCall(ISRestoreIndices(isicol, &ic));
356: PetscCall(ISRestoreIndices(isrow, &r));
358: PetscCall(ISIdentity(isrow, &row_identity));
359: PetscCall(ISIdentity(isicol, &col_identity));
360: if (b->inode.size_csr) {
361: C->ops->solve = MatSolve_SeqAIJ_Inode;
362: } else if (row_identity && col_identity) {
363: C->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
364: } else {
365: C->ops->solve = MatSolve_SeqAIJ;
366: }
367: C->ops->solveadd = MatSolveAdd_SeqAIJ;
368: C->ops->solvetranspose = MatSolveTranspose_SeqAIJ;
369: C->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ;
370: C->ops->matsolve = MatMatSolve_SeqAIJ;
371: C->ops->matsolvetranspose = MatMatSolveTranspose_SeqAIJ;
372: C->assembled = PETSC_TRUE;
373: C->preallocated = PETSC_TRUE;
375: PetscCall(PetscLogFlops(C->cmap->n));
377: /* MatShiftView(A,info,&sctx) */
378: if (sctx.nshift) {
379: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
380: 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));
381: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
382: PetscCall(PetscInfo(A, "number of shift_nz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
383: } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
384: PetscCall(PetscInfo(A, "number of shift_inblocks applied %" PetscInt_FMT ", each shift_amount %g\n", sctx.nshift, (double)info->shiftamount));
385: }
386: }
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: static PetscErrorCode MatMatSolve_SeqAIJ_inplace(Mat, Mat, Mat);
391: static PetscErrorCode MatSolve_SeqAIJ_NaturalOrdering_inplace(Mat, Vec, Vec);
392: static PetscErrorCode MatSolveAdd_SeqAIJ_inplace(Mat, Vec, Vec, Vec);
394: PetscErrorCode MatLUFactorNumeric_SeqAIJ_inplace(Mat B, Mat A, const MatFactorInfo *info)
395: {
396: Mat C = B;
397: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)C->data;
398: IS isrow = b->row, isicol = b->icol;
399: const PetscInt *r, *ic, *ics;
400: PetscInt nz, row, i, j, n = A->rmap->n, diag;
401: const PetscInt *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
402: const PetscInt *ajtmp, *bjtmp, *ddiag, *pj;
403: MatScalar *pv, *rtmp, *pc, multiplier, d;
404: const MatScalar *v, *aa;
405: MatScalar *ba;
406: PetscReal rs = 0.0;
407: FactorShiftCtx sctx;
408: PetscBool row_identity, col_identity;
410: PetscFunctionBegin;
411: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &ddiag, NULL));
413: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
414: PetscCall(MatSeqAIJGetArrayWrite(B, &ba));
415: /* MatPivotSetUp(): initialize shift context sctx */
416: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
418: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
419: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &ddiag, NULL));
420: sctx.shift_top = info->zeropivot;
421: for (i = 0; i < n; i++) {
422: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
423: d = aa[ddiag[i]];
424: rs = -PetscAbsScalar(d) - PetscRealPart(d);
425: v = aa + ai[i];
426: nz = ai[i + 1] - ai[i];
427: for (j = 0; j < nz; j++) rs += PetscAbsScalar(v[j]);
428: if (rs > sctx.shift_top) sctx.shift_top = rs;
429: }
430: sctx.shift_top *= 1.1;
431: sctx.nshift_max = 5;
432: sctx.shift_lo = 0.;
433: sctx.shift_hi = 1.;
434: }
436: PetscCall(ISGetIndices(isrow, &r));
437: PetscCall(ISGetIndices(isicol, &ic));
438: PetscCall(PetscMalloc1(n + 1, &rtmp));
439: ics = ic;
441: do {
442: sctx.newshift = PETSC_FALSE;
443: for (i = 0; i < n; i++) {
444: nz = bi[i + 1] - bi[i];
445: bjtmp = bj + bi[i];
446: for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;
448: /* load in initial (unfactored row) */
449: nz = ai[r[i] + 1] - ai[r[i]];
450: ajtmp = aj + ai[r[i]];
451: v = aa + ai[r[i]];
452: for (j = 0; j < nz; j++) rtmp[ics[ajtmp[j]]] = v[j];
453: rtmp[ics[r[i]]] += sctx.shift_amount; /* shift the diagonal of the matrix */
455: row = *bjtmp++;
456: while (row < i) {
457: pc = rtmp + row;
458: if (*pc != 0.0) {
459: pv = ba + ddiag[row];
460: pj = b->j + ddiag[row] + 1;
461: multiplier = *pc / *pv++;
462: *pc = multiplier;
463: nz = bi[row + 1] - ddiag[row] - 1;
464: for (j = 0; j < nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
465: PetscCall(PetscLogFlops(1 + 2.0 * nz));
466: }
467: row = *bjtmp++;
468: }
469: /* finished row so stick it into b->a */
470: pv = ba + bi[i];
471: pj = b->j + bi[i];
472: nz = bi[i + 1] - bi[i];
473: diag = ddiag[i] - bi[i];
474: rs = 0.0;
475: for (j = 0; j < nz; j++) {
476: pv[j] = rtmp[pj[j]];
477: rs += PetscAbsScalar(pv[j]);
478: }
479: rs -= PetscAbsScalar(pv[diag]);
481: sctx.rs = rs;
482: sctx.pv = pv[diag];
483: PetscCall(MatPivotCheck(B, A, info, &sctx, i));
484: if (sctx.newshift) break;
485: pv[diag] = sctx.pv;
486: }
488: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction > 0 && sctx.nshift < sctx.nshift_max) {
489: /*
490: * if no shift in this attempt & shifting & started shifting & can refine,
491: * then try lower shift
492: */
493: sctx.shift_hi = sctx.shift_fraction;
494: sctx.shift_fraction = (sctx.shift_hi + sctx.shift_lo) / 2.;
495: sctx.shift_amount = sctx.shift_fraction * sctx.shift_top;
496: sctx.newshift = PETSC_TRUE;
497: sctx.nshift++;
498: }
499: } while (sctx.newshift);
501: /* invert diagonal entries for simpler triangular solves */
502: for (i = 0; i < n; i++) ba[ddiag[i]] = 1.0 / ba[ddiag[i]];
503: PetscCall(PetscFree(rtmp));
504: PetscCall(ISRestoreIndices(isicol, &ic));
505: PetscCall(ISRestoreIndices(isrow, &r));
506: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
507: PetscCall(MatSeqAIJRestoreArrayWrite(B, &ba));
509: PetscCall(ISIdentity(isrow, &row_identity));
510: PetscCall(ISIdentity(isicol, &col_identity));
511: if (row_identity && col_identity) {
512: C->ops->solve = MatSolve_SeqAIJ_NaturalOrdering_inplace;
513: } else {
514: C->ops->solve = MatSolve_SeqAIJ_inplace;
515: }
516: C->ops->solveadd = MatSolveAdd_SeqAIJ_inplace;
517: C->ops->solvetranspose = MatSolveTranspose_SeqAIJ_inplace;
518: C->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ_inplace;
519: C->ops->matsolve = MatMatSolve_SeqAIJ_inplace;
520: C->ops->matsolvetranspose = NULL;
522: C->assembled = PETSC_TRUE;
523: C->preallocated = PETSC_TRUE;
525: PetscCall(PetscLogFlops(C->cmap->n));
526: if (sctx.nshift) {
527: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
528: 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));
529: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
530: PetscCall(PetscInfo(A, "number of shift_nz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
531: }
532: }
533: C->ops->solve = MatSolve_SeqAIJ_inplace;
534: C->ops->solvetranspose = MatSolveTranspose_SeqAIJ_inplace;
536: PetscCall(MatSeqAIJCheckInode(C));
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: static PetscErrorCode MatSolve_SeqAIJ_InplaceWithPerm(Mat, Vec, Vec);
542: /*
543: This routine implements inplace ILU(0) with row or/and column permutations.
544: Input:
545: A - original matrix
546: Output;
547: A - a->i (rowptr) is same as original rowptr, but factored i-the row is stored in rowperm[i]
548: a->j (col index) is permuted by the inverse of colperm, then sorted
549: a->a reordered accordingly with a->j
550: a->diag (ptr to diagonal elements) is updated.
551: */
552: PetscErrorCode MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(Mat B, Mat A, const MatFactorInfo *info)
553: {
554: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
555: IS isrow = a->row, isicol = a->icol;
556: const PetscInt *r, *ic, *ics;
557: PetscInt i, j, n = A->rmap->n, *ai = a->i, *aj = a->j;
558: PetscInt *ajtmp, nz, row;
559: PetscInt nbdiag, *pj;
560: PetscScalar *rtmp, *pc, multiplier, d;
561: MatScalar *pv, *v;
562: PetscReal rs;
563: FactorShiftCtx sctx;
564: MatScalar *aa, *vtmp;
565: PetscInt *diag;
567: PetscFunctionBegin;
568: PetscCheck(A == B, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "input and output matrix must have same address");
570: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, (const PetscInt **)&diag, NULL));
571: PetscCall(MatSeqAIJGetArray(A, &aa));
572: /* MatPivotSetUp(): initialize shift context sctx */
573: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
575: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
576: const PetscInt *ddiag;
578: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &ddiag, NULL));
579: sctx.shift_top = info->zeropivot;
580: for (i = 0; i < n; i++) {
581: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
582: d = aa[ddiag[i]];
583: rs = -PetscAbsScalar(d) - PetscRealPart(d);
584: vtmp = aa + ai[i];
585: nz = ai[i + 1] - ai[i];
586: for (j = 0; j < nz; j++) rs += PetscAbsScalar(vtmp[j]);
587: if (rs > sctx.shift_top) sctx.shift_top = rs;
588: }
589: sctx.shift_top *= 1.1;
590: sctx.nshift_max = 5;
591: sctx.shift_lo = 0.;
592: sctx.shift_hi = 1.;
593: }
595: PetscCall(ISGetIndices(isrow, &r));
596: PetscCall(ISGetIndices(isicol, &ic));
597: PetscCall(PetscMalloc1(n + 1, &rtmp));
598: PetscCall(PetscArrayzero(rtmp, n + 1));
599: ics = ic;
601: #if defined(MV)
602: sctx.shift_top = 0.;
603: sctx.nshift_max = 0;
604: sctx.shift_lo = 0.;
605: sctx.shift_hi = 0.;
606: sctx.shift_fraction = 0.;
608: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
609: sctx.shift_top = 0.;
610: for (i = 0; i < n; i++) {
611: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
612: d = aa[diag[i]];
613: rs = -PetscAbsScalar(d) - PetscRealPart(d);
614: v = aa + ai[i];
615: nz = ai[i + 1] - ai[i];
616: for (j = 0; j < nz; j++) rs += PetscAbsScalar(v[j]);
617: if (rs > sctx.shift_top) sctx.shift_top = rs;
618: }
619: if (sctx.shift_top < info->zeropivot) sctx.shift_top = info->zeropivot;
620: sctx.shift_top *= 1.1;
621: sctx.nshift_max = 5;
622: sctx.shift_lo = 0.;
623: sctx.shift_hi = 1.;
624: }
626: sctx.shift_amount = 0.;
627: sctx.nshift = 0;
628: #endif
630: do {
631: sctx.newshift = PETSC_FALSE;
632: for (i = 0; i < n; i++) {
633: /* load in initial unfactored row */
634: nz = ai[r[i] + 1] - ai[r[i]];
635: ajtmp = aj + ai[r[i]];
636: v = aa + ai[r[i]];
637: /* sort permuted ajtmp and values v accordingly */
638: for (j = 0; j < nz; j++) ajtmp[j] = ics[ajtmp[j]];
639: PetscCall(PetscSortIntWithScalarArray(nz, ajtmp, v));
641: diag[r[i]] = ai[r[i]];
642: for (j = 0; j < nz; j++) {
643: rtmp[ajtmp[j]] = v[j];
644: if (ajtmp[j] < i) diag[r[i]]++; /* update a->diag */
645: }
646: rtmp[r[i]] += sctx.shift_amount; /* shift the diagonal of the matrix */
648: row = *ajtmp++;
649: while (row < i) {
650: pc = rtmp + row;
651: if (*pc != 0.0) {
652: pv = aa + diag[r[row]];
653: pj = aj + diag[r[row]] + 1;
655: multiplier = *pc / *pv++;
656: *pc = multiplier;
657: nz = ai[r[row] + 1] - diag[r[row]] - 1;
658: for (j = 0; j < nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
659: PetscCall(PetscLogFlops(1 + 2.0 * nz));
660: }
661: row = *ajtmp++;
662: }
663: /* finished row so overwrite it onto aa */
664: pv = aa + ai[r[i]];
665: pj = aj + ai[r[i]];
666: nz = ai[r[i] + 1] - ai[r[i]];
667: nbdiag = diag[r[i]] - ai[r[i]]; /* num of entries before the diagonal */
669: rs = 0.0;
670: for (j = 0; j < nz; j++) {
671: pv[j] = rtmp[pj[j]];
672: if (j != nbdiag) rs += PetscAbsScalar(pv[j]);
673: }
675: sctx.rs = rs;
676: sctx.pv = pv[nbdiag];
677: PetscCall(MatPivotCheck(B, A, info, &sctx, i));
678: if (sctx.newshift) break;
679: pv[nbdiag] = sctx.pv;
680: }
682: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction > 0 && sctx.nshift < sctx.nshift_max) {
683: /*
684: * if no shift in this attempt & shifting & started shifting & can refine,
685: * then try lower shift
686: */
687: sctx.shift_hi = sctx.shift_fraction;
688: sctx.shift_fraction = (sctx.shift_hi + sctx.shift_lo) / 2.;
689: sctx.shift_amount = sctx.shift_fraction * sctx.shift_top;
690: sctx.newshift = PETSC_TRUE;
691: sctx.nshift++;
692: }
693: } while (sctx.newshift);
695: /* invert diagonal entries for simpler triangular solves */
696: for (i = 0; i < n; i++) aa[diag[r[i]]] = 1.0 / aa[diag[r[i]]];
698: PetscCall(MatSeqAIJRestoreArray(A, &aa));
699: PetscCall(PetscFree(rtmp));
700: PetscCall(ISRestoreIndices(isicol, &ic));
701: PetscCall(ISRestoreIndices(isrow, &r));
703: A->ops->solve = MatSolve_SeqAIJ_InplaceWithPerm;
704: A->ops->solveadd = MatSolveAdd_SeqAIJ_inplace;
705: A->ops->solvetranspose = MatSolveTranspose_SeqAIJ_inplace;
706: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ_inplace;
708: A->assembled = PETSC_TRUE;
709: A->preallocated = PETSC_TRUE;
711: PetscCall(PetscLogFlops(A->cmap->n));
712: if (sctx.nshift) {
713: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
714: 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));
715: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
716: PetscCall(PetscInfo(A, "number of shift_nz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
717: }
718: }
719: PetscFunctionReturn(PETSC_SUCCESS);
720: }
722: PetscErrorCode MatLUFactor_SeqAIJ(Mat A, IS row, IS col, const MatFactorInfo *info)
723: {
724: Mat C;
726: PetscFunctionBegin;
727: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &C));
728: PetscCall(MatLUFactorSymbolic(C, A, row, col, info));
729: PetscCall(MatLUFactorNumeric(C, A, info));
731: A->ops->solve = C->ops->solve;
732: A->ops->solvetranspose = C->ops->solvetranspose;
734: PetscCall(MatHeaderMerge(A, &C));
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }
738: PetscErrorCode MatSolve_SeqAIJ_inplace(Mat A, Vec bb, Vec xx)
739: {
740: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
741: IS iscol = a->col, isrow = a->row;
742: PetscInt i, n = A->rmap->n, *vi, *ai = a->i, *aj = a->j;
743: PetscInt nz;
744: const PetscInt *rout, *cout, *r, *c;
745: PetscScalar *x, *tmp, *tmps, sum;
746: const PetscScalar *b;
747: const MatScalar *aa, *v;
748: const PetscInt *adiag;
750: PetscFunctionBegin;
751: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
753: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
754: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
755: PetscCall(VecGetArrayRead(bb, &b));
756: PetscCall(VecGetArrayWrite(xx, &x));
757: tmp = a->solve_work;
759: PetscCall(ISGetIndices(isrow, &rout));
760: r = rout;
761: PetscCall(ISGetIndices(iscol, &cout));
762: c = cout + (n - 1);
764: /* forward solve the lower triangular */
765: tmp[0] = b[*r++];
766: tmps = tmp;
767: for (i = 1; i < n; i++) {
768: v = aa + ai[i];
769: vi = aj + ai[i];
770: nz = adiag[i] - ai[i];
771: sum = b[*r++];
772: PetscSparseDenseMinusDot(sum, tmps, v, vi, nz);
773: tmp[i] = sum;
774: }
776: /* backward solve the upper triangular */
777: for (i = n - 1; i >= 0; i--) {
778: v = aa + adiag[i] + 1;
779: vi = aj + adiag[i] + 1;
780: nz = ai[i + 1] - adiag[i] - 1;
781: sum = tmp[i];
782: PetscSparseDenseMinusDot(sum, tmps, v, vi, nz);
783: x[*c--] = tmp[i] = sum * aa[adiag[i]];
784: }
785: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
786: PetscCall(ISRestoreIndices(isrow, &rout));
787: PetscCall(ISRestoreIndices(iscol, &cout));
788: PetscCall(VecRestoreArrayRead(bb, &b));
789: PetscCall(VecRestoreArrayWrite(xx, &x));
790: PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: static PetscErrorCode MatMatSolve_SeqAIJ_inplace(Mat A, Mat B, Mat X)
795: {
796: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
797: IS iscol = a->col, isrow = a->row;
798: PetscInt i, n = A->rmap->n, *vi, *ai = a->i, *aj = a->j;
799: PetscInt nz, neq, ldb, ldx;
800: const PetscInt *rout, *cout, *r, *c;
801: PetscScalar *x, *tmp = a->solve_work, *tmps, sum;
802: const PetscScalar *b, *aa, *v;
803: PetscBool isdense;
804: const PetscInt *adiag;
806: PetscFunctionBegin;
807: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
808: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQDENSE, &isdense));
809: PetscCheck(isdense, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "B matrix must be a SeqDense matrix");
810: if (X != B) {
811: PetscCall(PetscObjectTypeCompare((PetscObject)X, MATSEQDENSE, &isdense));
812: PetscCheck(isdense, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "X matrix must be a SeqDense matrix");
813: }
814: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
815: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
816: PetscCall(MatDenseGetArrayRead(B, &b));
817: PetscCall(MatDenseGetLDA(B, &ldb));
818: PetscCall(MatDenseGetArray(X, &x));
819: PetscCall(MatDenseGetLDA(X, &ldx));
820: PetscCall(ISGetIndices(isrow, &rout));
821: r = rout;
822: PetscCall(ISGetIndices(iscol, &cout));
823: c = cout;
824: for (neq = 0; neq < B->cmap->n; neq++) {
825: /* forward solve the lower triangular */
826: tmp[0] = b[r[0]];
827: tmps = tmp;
828: for (i = 1; i < n; i++) {
829: v = aa + ai[i];
830: vi = aj + ai[i];
831: nz = adiag[i] - ai[i];
832: sum = b[r[i]];
833: PetscSparseDenseMinusDot(sum, tmps, v, vi, nz);
834: tmp[i] = sum;
835: }
836: /* backward solve the upper triangular */
837: for (i = n - 1; i >= 0; i--) {
838: v = aa + adiag[i] + 1;
839: vi = aj + adiag[i] + 1;
840: nz = ai[i + 1] - adiag[i] - 1;
841: sum = tmp[i];
842: PetscSparseDenseMinusDot(sum, tmps, v, vi, nz);
843: x[c[i]] = tmp[i] = sum * aa[adiag[i]];
844: }
845: b += ldb;
846: x += ldx;
847: }
848: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
849: PetscCall(ISRestoreIndices(isrow, &rout));
850: PetscCall(ISRestoreIndices(iscol, &cout));
851: PetscCall(MatDenseRestoreArrayRead(B, &b));
852: PetscCall(MatDenseRestoreArray(X, &x));
853: PetscCall(PetscLogFlops(B->cmap->n * (2.0 * a->nz - n)));
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: PetscErrorCode MatMatSolve_SeqAIJ(Mat A, Mat B, Mat X)
858: {
859: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
860: IS iscol = a->col, isrow = a->row;
861: PetscInt i, n = A->rmap->n, *vi, *ai = a->i, *aj = a->j;
862: const PetscInt *adiag;
863: PetscInt nz, neq, ldb, ldx;
864: const PetscInt *rout, *cout, *r, *c;
865: PetscScalar *x, *tmp = a->solve_work, sum;
866: const PetscScalar *b, *aa, *v;
867: PetscBool isdense;
869: PetscFunctionBegin;
870: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
871: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQDENSE, &isdense));
872: PetscCheck(isdense, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "B matrix must be a SeqDense matrix");
873: if (X != B) {
874: PetscCall(PetscObjectTypeCompare((PetscObject)X, MATSEQDENSE, &isdense));
875: PetscCheck(isdense, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "X matrix must be a SeqDense matrix");
876: }
877: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
878: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
879: PetscCall(MatDenseGetArrayRead(B, &b));
880: PetscCall(MatDenseGetLDA(B, &ldb));
881: PetscCall(MatDenseGetArray(X, &x));
882: PetscCall(MatDenseGetLDA(X, &ldx));
883: PetscCall(ISGetIndices(isrow, &rout));
884: r = rout;
885: PetscCall(ISGetIndices(iscol, &cout));
886: c = cout;
887: for (neq = 0; neq < B->cmap->n; neq++) {
888: /* forward solve the lower triangular */
889: tmp[0] = b[r[0]];
890: v = aa;
891: vi = aj;
892: for (i = 1; i < n; i++) {
893: nz = ai[i + 1] - ai[i];
894: sum = b[r[i]];
895: PetscSparseDenseMinusDot(sum, tmp, v, vi, nz);
896: tmp[i] = sum;
897: v += nz;
898: vi += nz;
899: }
900: /* backward solve the upper triangular */
901: for (i = n - 1; i >= 0; i--) {
902: v = aa + adiag[i + 1] + 1;
903: vi = aj + adiag[i + 1] + 1;
904: nz = adiag[i] - adiag[i + 1] - 1;
905: sum = tmp[i];
906: PetscSparseDenseMinusDot(sum, tmp, v, vi, nz);
907: x[c[i]] = tmp[i] = sum * v[nz]; /* v[nz] = aa[adiag[i]] */
908: }
909: b += ldb;
910: x += ldx;
911: }
912: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
913: PetscCall(ISRestoreIndices(isrow, &rout));
914: PetscCall(ISRestoreIndices(iscol, &cout));
915: PetscCall(MatDenseRestoreArrayRead(B, &b));
916: PetscCall(MatDenseRestoreArray(X, &x));
917: PetscCall(PetscLogFlops(B->cmap->n * (2.0 * a->nz - n)));
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
921: PetscErrorCode MatMatSolveTranspose_SeqAIJ(Mat A, Mat B, Mat X)
922: {
923: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
924: IS iscol = a->col, isrow = a->row;
925: PetscInt i, n = A->rmap->n, *vi, *ai = a->i, *aj = a->j, j;
926: const PetscInt *adiag = a->diag;
927: PetscInt nz, neq, ldb, ldx;
928: const PetscInt *rout, *cout, *r, *c;
929: PetscScalar *x, *tmp = a->solve_work, s1;
930: const PetscScalar *b, *aa, *v;
931: PetscBool isdense;
933: PetscFunctionBegin;
934: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
935: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQDENSE, &isdense));
936: PetscCheck(isdense, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "B matrix must be a SeqDense matrix");
937: if (X != B) {
938: PetscCall(PetscObjectTypeCompare((PetscObject)X, MATSEQDENSE, &isdense));
939: PetscCheck(isdense, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "X matrix must be a SeqDense matrix");
940: }
941: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
942: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
943: PetscCall(MatDenseGetArrayRead(B, &b));
944: PetscCall(MatDenseGetLDA(B, &ldb));
945: PetscCall(MatDenseGetArray(X, &x));
946: PetscCall(MatDenseGetLDA(X, &ldx));
947: PetscCall(ISGetIndices(isrow, &rout));
948: r = rout;
949: PetscCall(ISGetIndices(iscol, &cout));
950: c = cout;
951: for (neq = 0; neq < B->cmap->n; neq++) {
952: /* copy the b into temp work space according to permutation */
953: for (i = 0; i < n; i++) tmp[i] = b[c[i]];
955: /* forward solve the U^T */
956: for (i = 0; i < n; i++) {
957: v = aa + adiag[i + 1] + 1;
958: vi = aj + adiag[i + 1] + 1;
959: nz = adiag[i] - adiag[i + 1] - 1;
960: s1 = tmp[i];
961: s1 *= v[nz]; /* multiply by inverse of diagonal entry */
962: for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j];
963: tmp[i] = s1;
964: }
966: /* backward solve the L^T */
967: for (i = n - 1; i >= 0; i--) {
968: v = aa + ai[i];
969: vi = aj + ai[i];
970: nz = ai[i + 1] - ai[i];
971: s1 = tmp[i];
972: for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j];
973: }
975: /* copy tmp into x according to permutation */
976: for (i = 0; i < n; i++) x[r[i]] = tmp[i];
977: b += ldb;
978: x += ldx;
979: }
980: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
981: PetscCall(ISRestoreIndices(isrow, &rout));
982: PetscCall(ISRestoreIndices(iscol, &cout));
983: PetscCall(MatDenseRestoreArrayRead(B, &b));
984: PetscCall(MatDenseRestoreArray(X, &x));
985: PetscCall(PetscLogFlops(B->cmap->n * (2.0 * a->nz - n)));
986: PetscFunctionReturn(PETSC_SUCCESS);
987: }
989: static PetscErrorCode MatSolve_SeqAIJ_InplaceWithPerm(Mat A, Vec bb, Vec xx)
990: {
991: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
992: IS iscol = a->col, isrow = a->row;
993: const PetscInt *r, *c, *rout, *cout, *adiag;
994: PetscInt i, n = A->rmap->n, *vi, *ai = a->i, *aj = a->j;
995: PetscInt nz;
996: PetscScalar *x, *tmp, *tmps, sum;
997: const PetscScalar *b;
998: const MatScalar *aa, *v;
1000: PetscFunctionBegin;
1001: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
1003: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
1004: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1005: PetscCall(VecGetArrayRead(bb, &b));
1006: PetscCall(VecGetArrayWrite(xx, &x));
1007: tmp = a->solve_work;
1009: PetscCall(ISGetIndices(isrow, &rout));
1010: r = rout;
1011: PetscCall(ISGetIndices(iscol, &cout));
1012: c = cout + (n - 1);
1014: /* forward solve the lower triangular */
1015: tmp[0] = b[*r++];
1016: tmps = tmp;
1017: for (PetscInt row = 1; row < n; row++) {
1018: i = rout[row]; /* permuted row */
1019: v = aa + ai[i];
1020: vi = aj + ai[i];
1021: nz = adiag[i] - ai[i];
1022: sum = b[*r++];
1023: PetscSparseDenseMinusDot(sum, tmps, v, vi, nz);
1024: tmp[row] = sum;
1025: }
1027: /* backward solve the upper triangular */
1028: for (PetscInt row = n - 1; row >= 0; row--) {
1029: i = rout[row]; /* permuted row */
1030: v = aa + adiag[i] + 1;
1031: vi = aj + adiag[i] + 1;
1032: nz = ai[i + 1] - adiag[i] - 1;
1033: sum = tmp[row];
1034: PetscSparseDenseMinusDot(sum, tmps, v, vi, nz);
1035: x[*c--] = tmp[row] = sum * aa[adiag[i]];
1036: }
1037: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
1038: PetscCall(ISRestoreIndices(isrow, &rout));
1039: PetscCall(ISRestoreIndices(iscol, &cout));
1040: PetscCall(VecRestoreArrayRead(bb, &b));
1041: PetscCall(VecRestoreArrayWrite(xx, &x));
1042: PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
1043: PetscFunctionReturn(PETSC_SUCCESS);
1044: }
1046: #include <../src/mat/impls/aij/seq/ftn-kernels/fsolve.h>
1047: static PetscErrorCode MatSolve_SeqAIJ_NaturalOrdering_inplace(Mat A, Vec bb, Vec xx)
1048: {
1049: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1050: PetscInt n = A->rmap->n;
1051: const PetscInt *ai = a->i, *aj = a->j, *adiag;
1052: PetscScalar *x;
1053: const PetscScalar *b;
1054: const MatScalar *aa;
1055: #if !PetscDefined(USE_FORTRAN_KERNEL_SOLVEAIJ)
1056: PetscInt adiag_i, i, nz, ai_i;
1057: const PetscInt *vi;
1058: const MatScalar *v;
1059: PetscScalar sum;
1060: #endif
1062: PetscFunctionBegin;
1063: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
1064: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
1065: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1066: PetscCall(VecGetArrayRead(bb, &b));
1067: PetscCall(VecGetArrayWrite(xx, &x));
1069: #if PetscDefined(USE_FORTRAN_KERNEL_SOLVEAIJ)
1070: fortransolveaij_(&n, x, ai, aj, adiag, aa, b);
1071: #else
1072: /* forward solve the lower triangular */
1073: x[0] = b[0];
1074: for (i = 1; i < n; i++) {
1075: ai_i = ai[i];
1076: v = aa + ai_i;
1077: vi = aj + ai_i;
1078: nz = adiag[i] - ai_i;
1079: sum = b[i];
1080: PetscSparseDenseMinusDot(sum, x, v, vi, nz);
1081: x[i] = sum;
1082: }
1084: /* backward solve the upper triangular */
1085: for (i = n - 1; i >= 0; i--) {
1086: adiag_i = adiag[i];
1087: v = aa + adiag_i + 1;
1088: vi = aj + adiag_i + 1;
1089: nz = ai[i + 1] - adiag_i - 1;
1090: sum = x[i];
1091: PetscSparseDenseMinusDot(sum, x, v, vi, nz);
1092: x[i] = sum * aa[adiag_i];
1093: }
1094: #endif
1095: PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
1096: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1097: PetscCall(VecRestoreArrayRead(bb, &b));
1098: PetscCall(VecRestoreArrayWrite(xx, &x));
1099: PetscFunctionReturn(PETSC_SUCCESS);
1100: }
1102: static PetscErrorCode MatSolveAdd_SeqAIJ_inplace(Mat A, Vec bb, Vec yy, Vec xx)
1103: {
1104: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1105: IS iscol = a->col, isrow = a->row;
1106: PetscInt i, n = A->rmap->n, j;
1107: PetscInt nz;
1108: const PetscInt *rout, *cout, *r, *c, *vi, *ai = a->i, *aj = a->j, *adiag;
1109: PetscScalar *x, *tmp, sum;
1110: const PetscScalar *b;
1111: const MatScalar *aa, *v;
1113: PetscFunctionBegin;
1114: if (yy != xx) PetscCall(VecCopy(yy, xx));
1116: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
1117: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1118: PetscCall(VecGetArrayRead(bb, &b));
1119: PetscCall(VecGetArray(xx, &x));
1120: tmp = a->solve_work;
1122: PetscCall(ISGetIndices(isrow, &rout));
1123: r = rout;
1124: PetscCall(ISGetIndices(iscol, &cout));
1125: c = cout + (n - 1);
1127: /* forward solve the lower triangular */
1128: tmp[0] = b[*r++];
1129: for (i = 1; i < n; i++) {
1130: v = aa + ai[i];
1131: vi = aj + ai[i];
1132: nz = adiag[i] - ai[i];
1133: sum = b[*r++];
1134: for (j = 0; j < nz; j++) sum -= v[j] * tmp[vi[j]];
1135: tmp[i] = sum;
1136: }
1138: /* backward solve the upper triangular */
1139: for (i = n - 1; i >= 0; i--) {
1140: v = aa + adiag[i] + 1;
1141: vi = aj + adiag[i] + 1;
1142: nz = ai[i + 1] - adiag[i] - 1;
1143: sum = tmp[i];
1144: for (j = 0; j < nz; j++) sum -= v[j] * tmp[vi[j]];
1145: tmp[i] = sum * aa[adiag[i]];
1146: x[*c--] += tmp[i];
1147: }
1149: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
1150: PetscCall(ISRestoreIndices(isrow, &rout));
1151: PetscCall(ISRestoreIndices(iscol, &cout));
1152: PetscCall(VecRestoreArrayRead(bb, &b));
1153: PetscCall(VecRestoreArray(xx, &x));
1154: PetscCall(PetscLogFlops(2.0 * a->nz));
1155: PetscFunctionReturn(PETSC_SUCCESS);
1156: }
1158: PetscErrorCode MatSolveAdd_SeqAIJ(Mat A, Vec bb, Vec yy, Vec xx)
1159: {
1160: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1161: IS iscol = a->col, isrow = a->row;
1162: PetscInt i, n = A->rmap->n, j;
1163: PetscInt nz;
1164: const PetscInt *rout, *cout, *r, *c, *vi, *ai = a->i, *aj = a->j, *adiag;
1165: PetscScalar *x, *tmp, sum;
1166: const PetscScalar *b;
1167: const MatScalar *aa, *v;
1169: PetscFunctionBegin;
1170: if (yy != xx) PetscCall(VecCopy(yy, xx));
1172: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
1173: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1174: PetscCall(VecGetArrayRead(bb, &b));
1175: PetscCall(VecGetArray(xx, &x));
1176: tmp = a->solve_work;
1178: PetscCall(ISGetIndices(isrow, &rout));
1179: r = rout;
1180: PetscCall(ISGetIndices(iscol, &cout));
1181: c = cout;
1183: /* forward solve the lower triangular */
1184: tmp[0] = b[r[0]];
1185: v = aa;
1186: vi = aj;
1187: for (i = 1; i < n; i++) {
1188: nz = ai[i + 1] - ai[i];
1189: sum = b[r[i]];
1190: for (j = 0; j < nz; j++) sum -= v[j] * tmp[vi[j]];
1191: tmp[i] = sum;
1192: v += nz;
1193: vi += nz;
1194: }
1196: /* backward solve the upper triangular */
1197: v = aa + adiag[n - 1];
1198: vi = aj + adiag[n - 1];
1199: for (i = n - 1; i >= 0; i--) {
1200: nz = adiag[i] - adiag[i + 1] - 1;
1201: sum = tmp[i];
1202: for (j = 0; j < nz; j++) sum -= v[j] * tmp[vi[j]];
1203: tmp[i] = sum * v[nz];
1204: x[c[i]] += tmp[i];
1205: v += nz + 1;
1206: vi += nz + 1;
1207: }
1209: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
1210: PetscCall(ISRestoreIndices(isrow, &rout));
1211: PetscCall(ISRestoreIndices(iscol, &cout));
1212: PetscCall(VecRestoreArrayRead(bb, &b));
1213: PetscCall(VecRestoreArray(xx, &x));
1214: PetscCall(PetscLogFlops(2.0 * a->nz));
1215: PetscFunctionReturn(PETSC_SUCCESS);
1216: }
1218: PetscErrorCode MatSolveTranspose_SeqAIJ_inplace(Mat A, Vec bb, Vec xx)
1219: {
1220: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1221: IS iscol = a->col, isrow = a->row;
1222: const PetscInt *rout, *cout, *r, *c, *diag = a->diag, *ai = a->i, *aj = a->j, *vi;
1223: PetscInt i, n = A->rmap->n, j;
1224: PetscInt nz;
1225: PetscScalar *x, *tmp, s1;
1226: const MatScalar *aa, *v;
1227: const PetscScalar *b;
1229: PetscFunctionBegin;
1230: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1231: PetscCall(VecGetArrayRead(bb, &b));
1232: PetscCall(VecGetArrayWrite(xx, &x));
1233: tmp = a->solve_work;
1235: PetscCall(ISGetIndices(isrow, &rout));
1236: r = rout;
1237: PetscCall(ISGetIndices(iscol, &cout));
1238: c = cout;
1240: /* copy the b into temp work space according to permutation */
1241: for (i = 0; i < n; i++) tmp[i] = b[c[i]];
1243: /* forward solve the U^T */
1244: for (i = 0; i < n; i++) {
1245: v = aa + diag[i];
1246: vi = aj + diag[i] + 1;
1247: nz = ai[i + 1] - diag[i] - 1;
1248: s1 = tmp[i];
1249: s1 *= (*v++); /* multiply by inverse of diagonal entry */
1250: for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j];
1251: tmp[i] = s1;
1252: }
1254: /* backward solve the L^T */
1255: for (i = n - 1; i >= 0; i--) {
1256: v = aa + diag[i] - 1;
1257: vi = aj + diag[i] - 1;
1258: nz = diag[i] - ai[i];
1259: s1 = tmp[i];
1260: for (j = 0; j > -nz; j--) tmp[vi[j]] -= s1 * v[j];
1261: }
1263: /* copy tmp into x according to permutation */
1264: for (i = 0; i < n; i++) x[r[i]] = tmp[i];
1266: PetscCall(ISRestoreIndices(isrow, &rout));
1267: PetscCall(ISRestoreIndices(iscol, &cout));
1268: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
1269: PetscCall(VecRestoreArrayRead(bb, &b));
1270: PetscCall(VecRestoreArrayWrite(xx, &x));
1272: PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
1273: PetscFunctionReturn(PETSC_SUCCESS);
1274: }
1276: PetscErrorCode MatSolveTranspose_SeqAIJ(Mat A, Vec bb, Vec xx)
1277: {
1278: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1279: IS iscol = a->col, isrow = a->row;
1280: const PetscInt *rout, *cout, *r, *c, *adiag, *ai = a->i, *aj = a->j, *vi;
1281: PetscInt i, n = A->rmap->n, j;
1282: PetscInt nz;
1283: PetscScalar *x, *tmp, s1;
1284: const MatScalar *aa, *v;
1285: const PetscScalar *b;
1287: PetscFunctionBegin;
1288: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
1289: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1290: PetscCall(VecGetArrayRead(bb, &b));
1291: PetscCall(VecGetArrayWrite(xx, &x));
1292: tmp = a->solve_work;
1294: PetscCall(ISGetIndices(isrow, &rout));
1295: r = rout;
1296: PetscCall(ISGetIndices(iscol, &cout));
1297: c = cout;
1299: /* copy the b into temp work space according to permutation */
1300: for (i = 0; i < n; i++) tmp[i] = b[c[i]];
1302: /* forward solve the U^T */
1303: for (i = 0; i < n; i++) {
1304: v = aa + adiag[i + 1] + 1;
1305: vi = aj + adiag[i + 1] + 1;
1306: nz = adiag[i] - adiag[i + 1] - 1;
1307: s1 = tmp[i];
1308: s1 *= v[nz]; /* multiply by inverse of diagonal entry */
1309: for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j];
1310: tmp[i] = s1;
1311: }
1313: /* backward solve the L^T */
1314: for (i = n - 1; i >= 0; i--) {
1315: v = aa + ai[i];
1316: vi = aj + ai[i];
1317: nz = ai[i + 1] - ai[i];
1318: s1 = tmp[i];
1319: for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j];
1320: }
1322: /* copy tmp into x according to permutation */
1323: for (i = 0; i < n; i++) x[r[i]] = tmp[i];
1325: PetscCall(ISRestoreIndices(isrow, &rout));
1326: PetscCall(ISRestoreIndices(iscol, &cout));
1327: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
1328: PetscCall(VecRestoreArrayRead(bb, &b));
1329: PetscCall(VecRestoreArrayWrite(xx, &x));
1331: PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
1332: PetscFunctionReturn(PETSC_SUCCESS);
1333: }
1335: PetscErrorCode MatSolveTransposeAdd_SeqAIJ_inplace(Mat A, Vec bb, Vec zz, Vec xx)
1336: {
1337: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1338: IS iscol = a->col, isrow = a->row;
1339: const PetscInt *rout, *cout, *r, *c, *diag = a->diag, *ai = a->i, *aj = a->j, *vi;
1340: PetscInt i, n = A->rmap->n, j;
1341: PetscInt nz;
1342: PetscScalar *x, *tmp, s1;
1343: const MatScalar *aa, *v;
1344: const PetscScalar *b;
1346: PetscFunctionBegin;
1347: if (zz != xx) PetscCall(VecCopy(zz, xx));
1348: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1349: PetscCall(VecGetArrayRead(bb, &b));
1350: PetscCall(VecGetArray(xx, &x));
1351: tmp = a->solve_work;
1353: PetscCall(ISGetIndices(isrow, &rout));
1354: r = rout;
1355: PetscCall(ISGetIndices(iscol, &cout));
1356: c = cout;
1358: /* copy the b into temp work space according to permutation */
1359: for (i = 0; i < n; i++) tmp[i] = b[c[i]];
1361: /* forward solve the U^T */
1362: for (i = 0; i < n; i++) {
1363: v = aa + diag[i];
1364: vi = aj + diag[i] + 1;
1365: nz = ai[i + 1] - diag[i] - 1;
1366: s1 = tmp[i];
1367: s1 *= (*v++); /* multiply by inverse of diagonal entry */
1368: for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j];
1369: tmp[i] = s1;
1370: }
1372: /* backward solve the L^T */
1373: for (i = n - 1; i >= 0; i--) {
1374: v = aa + diag[i] - 1;
1375: vi = aj + diag[i] - 1;
1376: nz = diag[i] - ai[i];
1377: s1 = tmp[i];
1378: for (j = 0; j > -nz; j--) tmp[vi[j]] -= s1 * v[j];
1379: }
1381: /* copy tmp into x according to permutation */
1382: for (i = 0; i < n; i++) x[r[i]] += tmp[i];
1384: PetscCall(ISRestoreIndices(isrow, &rout));
1385: PetscCall(ISRestoreIndices(iscol, &cout));
1386: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
1387: PetscCall(VecRestoreArrayRead(bb, &b));
1388: PetscCall(VecRestoreArray(xx, &x));
1390: PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
1391: PetscFunctionReturn(PETSC_SUCCESS);
1392: }
1394: PetscErrorCode MatSolveTransposeAdd_SeqAIJ(Mat A, Vec bb, Vec zz, Vec xx)
1395: {
1396: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1397: IS iscol = a->col, isrow = a->row;
1398: const PetscInt *rout, *cout, *r, *c, *adiag, *ai = a->i, *aj = a->j, *vi;
1399: PetscInt i, n = A->rmap->n, j;
1400: PetscInt nz;
1401: PetscScalar *x, *tmp, s1;
1402: const MatScalar *aa, *v;
1403: const PetscScalar *b;
1405: PetscFunctionBegin;
1406: if (zz != xx) PetscCall(VecCopy(zz, xx));
1407: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
1408: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1409: PetscCall(VecGetArrayRead(bb, &b));
1410: PetscCall(VecGetArray(xx, &x));
1411: tmp = a->solve_work;
1413: PetscCall(ISGetIndices(isrow, &rout));
1414: r = rout;
1415: PetscCall(ISGetIndices(iscol, &cout));
1416: c = cout;
1418: /* copy the b into temp work space according to permutation */
1419: for (i = 0; i < n; i++) tmp[i] = b[c[i]];
1421: /* forward solve the U^T */
1422: for (i = 0; i < n; i++) {
1423: v = aa + adiag[i + 1] + 1;
1424: vi = aj + adiag[i + 1] + 1;
1425: nz = adiag[i] - adiag[i + 1] - 1;
1426: s1 = tmp[i];
1427: s1 *= v[nz]; /* multiply by inverse of diagonal entry */
1428: for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j];
1429: tmp[i] = s1;
1430: }
1432: /* backward solve the L^T */
1433: for (i = n - 1; i >= 0; i--) {
1434: v = aa + ai[i];
1435: vi = aj + ai[i];
1436: nz = ai[i + 1] - ai[i];
1437: s1 = tmp[i];
1438: for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j];
1439: }
1441: /* copy tmp into x according to permutation */
1442: for (i = 0; i < n; i++) x[r[i]] += tmp[i];
1444: PetscCall(ISRestoreIndices(isrow, &rout));
1445: PetscCall(ISRestoreIndices(iscol, &cout));
1446: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
1447: PetscCall(VecRestoreArrayRead(bb, &b));
1448: PetscCall(VecRestoreArray(xx, &x));
1450: PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
1451: PetscFunctionReturn(PETSC_SUCCESS);
1452: }
1454: /*
1455: ilu() under revised new data structure.
1456: Factored arrays bj and ba are stored as
1457: L(0,:), L(1,:), ...,L(n-1,:), U(n-1,:),...,U(i,:),U(i-1,:),...,U(0,:)
1459: bi=fact->i is an array of size n+1, in which
1460: bi[i]: points to 1st entry of L(i,:),i=0,...,n-1
1461: bi[n]: points to L(n-1,n-1)+1
1463: bdiag=fact->diag is an array of size n+1,in which
1464: bdiag[i]: points to diagonal of U(i,:), i=0,...,n-1
1465: bdiag[n]: points to entry of U(n-1,0)-1
1467: U(i,:) contains bdiag[i] as its last entry, i.e.,
1468: U(i,:) = (u[i,i+1],...,u[i,n-1],diag[i])
1469: */
1470: PetscErrorCode MatILUFactorSymbolic_SeqAIJ_ilu0(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
1471: {
1472: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b;
1473: const PetscInt n = A->rmap->n, *ai = a->i, *aj, *adiag;
1474: PetscInt i, j, k = 0, nz, *bi, *bj, *bdiag;
1475: IS isicol;
1477: PetscFunctionBegin;
1478: PetscCall(ISInvertPermutation(iscol, PETSC_DECIDE, &isicol));
1479: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
1480: PetscCall(MatDuplicateNoCreate_SeqAIJ(fact, A, MAT_DO_NOT_COPY_VALUES, PETSC_FALSE));
1481: b = (Mat_SeqAIJ *)fact->data;
1483: /* allocate matrix arrays for new data structure */
1484: PetscCall(PetscShmgetAllocateArray(ai[n], sizeof(PetscScalar), (void **)&b->a));
1485: PetscCall(PetscShmgetAllocateArray(ai[n], sizeof(PetscInt), (void **)&b->j));
1486: PetscCall(PetscShmgetAllocateArray(n + 1, sizeof(PetscInt), (void **)&b->i));
1487: if (n > 0) PetscCall(PetscArrayzero(b->a, ai[n]));
1488: b->free_a = PETSC_TRUE;
1489: b->free_ij = PETSC_TRUE;
1491: if (!b->diag) PetscCall(PetscMalloc1(n + 1, &b->diag));
1492: bdiag = b->diag;
1494: /* set bi and bj with new data structure */
1495: bi = b->i;
1496: bj = b->j;
1498: /* L part */
1499: bi[0] = 0;
1500: for (i = 0; i < n; i++) {
1501: nz = adiag[i] - ai[i];
1502: bi[i + 1] = bi[i] + nz;
1503: aj = a->j + ai[i];
1504: for (j = 0; j < nz; j++) bj[k++] = aj[j];
1505: }
1507: /* U part */
1508: bdiag[n] = bi[n] - 1;
1509: for (i = n - 1; i >= 0; i--) {
1510: nz = ai[i + 1] - adiag[i] - 1;
1511: aj = a->j + adiag[i] + 1;
1512: for (j = 0; j < nz; j++) bj[k++] = aj[j];
1513: /* diag[i] */
1514: bj[k++] = i;
1515: bdiag[i] = bdiag[i + 1] + nz + 1;
1516: }
1518: fact->factortype = MAT_FACTOR_ILU;
1519: fact->info.factor_mallocs = 0;
1520: fact->info.fill_ratio_given = info->fill;
1521: fact->info.fill_ratio_needed = 1.0;
1522: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ;
1523: PetscCall(MatSeqAIJCheckInode_FactorLU(fact));
1525: b = (Mat_SeqAIJ *)fact->data;
1526: b->row = isrow;
1527: b->col = iscol;
1528: b->icol = isicol;
1529: PetscCall(PetscMalloc1(fact->rmap->n, &b->solve_work));
1530: PetscCall(PetscObjectReference((PetscObject)isrow));
1531: PetscCall(PetscObjectReference((PetscObject)iscol));
1532: PetscFunctionReturn(PETSC_SUCCESS);
1533: }
1535: PetscErrorCode MatILUFactorSymbolic_SeqAIJ(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
1536: {
1537: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b;
1538: IS isicol;
1539: const PetscInt *r, *ic;
1540: PetscInt n = A->rmap->n, *ai = a->i, *aj = a->j;
1541: PetscInt *bi, *cols, nnz, *cols_lvl;
1542: PetscInt *bdiag, prow, fm, nzbd, reallocs = 0, dcount = 0;
1543: PetscInt i, levels = (PetscInt)info->levels, diagonal_fill;
1544: PetscBool col_identity, row_identity;
1545: PetscReal f;
1546: PetscInt nlnk, *lnk, *lnk_lvl = NULL;
1547: PetscBT lnkbt;
1548: PetscInt nzi, *bj, **bj_ptr, **bjlvl_ptr;
1549: PetscFreeSpaceList free_space = NULL, current_space = NULL;
1550: PetscFreeSpaceList free_space_lvl = NULL, current_space_lvl = NULL;
1551: PetscBool diagDense;
1553: PetscFunctionBegin;
1554: 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);
1555: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, NULL, &diagDense));
1556: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entries");
1558: PetscCall(ISIdentity(isrow, &row_identity));
1559: PetscCall(ISIdentity(iscol, &col_identity));
1560: if (!levels && row_identity && col_identity) {
1561: /* special case: ilu(0) with natural ordering */
1562: PetscCall(MatILUFactorSymbolic_SeqAIJ_ilu0(fact, A, isrow, iscol, info));
1563: if (a->inode.size_csr) fact->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
1564: PetscFunctionReturn(PETSC_SUCCESS);
1565: }
1567: PetscCall(ISInvertPermutation(iscol, PETSC_DECIDE, &isicol));
1568: PetscCall(ISGetIndices(isrow, &r));
1569: PetscCall(ISGetIndices(isicol, &ic));
1571: /* get new row and diagonal pointers, must be allocated separately because they will be given to the Mat_SeqAIJ and freed separately */
1572: PetscCall(PetscShmgetAllocateArray(n + 1, sizeof(PetscInt), (void **)&bi));
1573: PetscCall(PetscMalloc1(n + 1, &bdiag));
1574: bi[0] = bdiag[0] = 0;
1575: PetscCall(PetscMalloc2(n, &bj_ptr, n, &bjlvl_ptr));
1577: /* create a linked list for storing column indices of the active row */
1578: nlnk = n + 1;
1579: PetscCall(PetscIncompleteLLCreate(n, n, nlnk, lnk, lnk_lvl, lnkbt));
1581: /* initial FreeSpace size is f*(ai[n]+1) */
1582: f = info->fill;
1583: diagonal_fill = (PetscInt)info->diagonal_fill;
1584: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space));
1585: current_space = free_space;
1586: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space_lvl));
1587: current_space_lvl = free_space_lvl;
1588: for (i = 0; i < n; i++) {
1589: nzi = 0;
1590: /* copy current row into linked list */
1591: nnz = ai[r[i] + 1] - ai[r[i]];
1592: 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);
1593: cols = aj + ai[r[i]];
1594: lnk[i] = -1; /* marker to indicate if diagonal exists */
1595: PetscCall(PetscIncompleteLLInit(nnz, cols, n, ic, &nlnk, lnk, lnk_lvl, lnkbt));
1596: nzi += nlnk;
1598: /* make sure diagonal entry is included */
1599: if (diagonal_fill && lnk[i] == -1) {
1600: fm = n;
1601: while (lnk[fm] < i) fm = lnk[fm];
1602: lnk[i] = lnk[fm]; /* insert diagonal into linked list */
1603: lnk[fm] = i;
1604: lnk_lvl[i] = 0;
1605: nzi++;
1606: dcount++;
1607: }
1609: /* add pivot rows into the active row */
1610: nzbd = 0;
1611: prow = lnk[n];
1612: while (prow < i) {
1613: nnz = bdiag[prow];
1614: cols = bj_ptr[prow] + nnz + 1;
1615: cols_lvl = bjlvl_ptr[prow] + nnz + 1;
1616: nnz = bi[prow + 1] - bi[prow] - nnz - 1;
1617: PetscCall(PetscILULLAddSorted(nnz, cols, levels, cols_lvl, prow, &nlnk, lnk, lnk_lvl, lnkbt, prow));
1618: nzi += nlnk;
1619: prow = lnk[prow];
1620: nzbd++;
1621: }
1622: bdiag[i] = nzbd;
1623: bi[i + 1] = bi[i] + nzi;
1624: /* if free space is not available, make more free space */
1625: if (current_space->local_remaining < nzi) {
1626: nnz = PetscIntMultTruncate(2, PetscIntMultTruncate(nzi, n - i)); /* estimated and max additional space needed */
1627: PetscCall(PetscFreeSpaceGet(nnz, ¤t_space));
1628: PetscCall(PetscFreeSpaceGet(nnz, ¤t_space_lvl));
1629: reallocs++;
1630: }
1632: /* copy data into free_space and free_space_lvl, then initialize lnk */
1633: PetscCall(PetscIncompleteLLClean(n, n, nzi, lnk, lnk_lvl, current_space->array, current_space_lvl->array, lnkbt));
1634: bj_ptr[i] = current_space->array;
1635: bjlvl_ptr[i] = current_space_lvl->array;
1637: /* make sure the active row i has diagonal entry */
1638: 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);
1640: current_space->array += nzi;
1641: current_space->local_used += nzi;
1642: current_space->local_remaining -= nzi;
1643: current_space_lvl->array += nzi;
1644: current_space_lvl->local_used += nzi;
1645: current_space_lvl->local_remaining -= nzi;
1646: }
1648: PetscCall(ISRestoreIndices(isrow, &r));
1649: PetscCall(ISRestoreIndices(isicol, &ic));
1650: /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
1651: PetscCall(PetscShmgetAllocateArray(bi[n], sizeof(PetscInt), (void **)&bj));
1652: PetscCall(PetscFreeSpaceContiguous_LU(&free_space, bj, n, bi, bdiag));
1654: PetscCall(PetscIncompleteLLDestroy(lnk, lnkbt));
1655: PetscCall(PetscFreeSpaceDestroy(free_space_lvl));
1656: PetscCall(PetscFree2(bj_ptr, bjlvl_ptr));
1658: #if PetscDefined(USE_INFO)
1659: {
1660: PetscReal af = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);
1661: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)f, (double)af));
1662: PetscCall(PetscInfo(A, "Run with -[sub_]pc_factor_fill %g or use \n", (double)af));
1663: PetscCall(PetscInfo(A, "PCFactorSetFill([sub]pc,%g);\n", (double)af));
1664: PetscCall(PetscInfo(A, "for best performance.\n"));
1665: if (diagonal_fill) PetscCall(PetscInfo(A, "Detected and replaced %" PetscInt_FMT " missing diagonals\n", dcount));
1666: }
1667: #endif
1668: /* put together the new matrix */
1669: PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(fact, MAT_SKIP_ALLOCATION, NULL));
1670: b = (Mat_SeqAIJ *)fact->data;
1671: b->free_ij = PETSC_TRUE;
1672: PetscCall(PetscShmgetAllocateArray(bdiag[0] + 1, sizeof(PetscScalar), (void **)&b->a));
1673: b->free_a = PETSC_TRUE;
1674: b->j = bj;
1675: b->i = bi;
1676: b->diag = bdiag;
1677: b->ilen = NULL;
1678: b->imax = NULL;
1679: b->row = isrow;
1680: b->col = iscol;
1681: PetscCall(PetscObjectReference((PetscObject)isrow));
1682: PetscCall(PetscObjectReference((PetscObject)iscol));
1683: b->icol = isicol;
1685: PetscCall(PetscMalloc1(n, &b->solve_work));
1686: /* In b structure: Free imax, ilen, old a, old j.
1687: Allocate bdiag, solve_work, new a, new j */
1688: b->maxnz = b->nz = bdiag[0] + 1;
1690: fact->info.factor_mallocs = reallocs;
1691: fact->info.fill_ratio_given = f;
1692: fact->info.fill_ratio_needed = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);
1693: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ;
1694: if (a->inode.size_csr) fact->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
1695: PetscCall(MatSeqAIJCheckInode_FactorLU(fact));
1696: PetscFunctionReturn(PETSC_SUCCESS);
1697: }
1699: PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ(Mat B, Mat A, const MatFactorInfo *info)
1700: {
1701: Mat C = B;
1702: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1703: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)C->data;
1704: IS ip = b->row, iip = b->icol;
1705: const PetscInt *rip, *riip;
1706: PetscInt i, j, mbs = A->rmap->n, *bi = b->i, *bj = b->j, *bdiag = b->diag, *bjtmp;
1707: PetscInt *ai = a->i, *aj = a->j;
1708: PetscInt k, jmin, jmax, *c2r, *il, col, nexti, ili, nz;
1709: MatScalar *rtmp, *ba = b->a, *bval, dk, uikdi;
1710: PetscBool perm_identity;
1711: FactorShiftCtx sctx;
1712: PetscReal rs;
1713: const MatScalar *aa, *v;
1714: MatScalar d;
1715: const PetscInt *adiag;
1717: PetscFunctionBegin;
1718: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
1719: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1720: /* MatPivotSetUp(): initialize shift context sctx */
1721: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
1723: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
1724: sctx.shift_top = info->zeropivot;
1725: for (i = 0; i < mbs; i++) {
1726: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
1727: d = aa[adiag[i]];
1728: rs = -PetscAbsScalar(d) - PetscRealPart(d);
1729: v = aa + ai[i];
1730: nz = ai[i + 1] - ai[i];
1731: for (j = 0; j < nz; j++) rs += PetscAbsScalar(v[j]);
1732: if (rs > sctx.shift_top) sctx.shift_top = rs;
1733: }
1734: sctx.shift_top *= 1.1;
1735: sctx.nshift_max = 5;
1736: sctx.shift_lo = 0.;
1737: sctx.shift_hi = 1.;
1738: }
1740: PetscCall(ISGetIndices(ip, &rip));
1741: PetscCall(ISGetIndices(iip, &riip));
1743: /* allocate working arrays
1744: c2r: linked list, keep track of pivot rows for a given column. c2r[col]: head of the list for a given col
1745: il: for active k row, il[i] gives the index of the 1st nonzero entry in U[i,k:n-1] in bj and ba arrays
1746: */
1747: PetscCall(PetscMalloc3(mbs, &rtmp, mbs, &il, mbs, &c2r));
1749: do {
1750: sctx.newshift = PETSC_FALSE;
1752: for (i = 0; i < mbs; i++) c2r[i] = mbs;
1753: if (mbs) il[0] = 0;
1755: for (k = 0; k < mbs; k++) {
1756: /* zero rtmp */
1757: nz = bi[k + 1] - bi[k];
1758: bjtmp = bj + bi[k];
1759: for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;
1761: /* load in initial unfactored row */
1762: bval = ba + bi[k];
1763: jmin = ai[rip[k]];
1764: jmax = ai[rip[k] + 1];
1765: for (j = jmin; j < jmax; j++) {
1766: col = riip[aj[j]];
1767: if (col >= k) { /* only take upper triangular entry */
1768: rtmp[col] = aa[j];
1769: *bval++ = 0.0; /* for in-place factorization */
1770: }
1771: }
1772: /* shift the diagonal of the matrix: ZeropivotApply() */
1773: rtmp[k] += sctx.shift_amount; /* shift the diagonal of the matrix */
1775: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1776: dk = rtmp[k];
1777: i = c2r[k]; /* first row to be added to k_th row */
1779: while (i < k) {
1780: nexti = c2r[i]; /* next row to be added to k_th row */
1782: /* compute multiplier, update diag(k) and U(i,k) */
1783: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1784: uikdi = -ba[ili] * ba[bdiag[i]]; /* diagonal(k) */
1785: dk += uikdi * ba[ili]; /* update diag[k] */
1786: ba[ili] = uikdi; /* -U(i,k) */
1788: /* add multiple of row i to k-th row */
1789: jmin = ili + 1;
1790: jmax = bi[i + 1];
1791: if (jmin < jmax) {
1792: for (j = jmin; j < jmax; j++) rtmp[bj[j]] += uikdi * ba[j];
1793: /* update il and c2r for row i */
1794: il[i] = jmin;
1795: j = bj[jmin];
1796: c2r[i] = c2r[j];
1797: c2r[j] = i;
1798: }
1799: i = nexti;
1800: }
1802: /* copy data into U(k,:) */
1803: rs = 0.0;
1804: jmin = bi[k];
1805: jmax = bi[k + 1] - 1;
1806: if (jmin < jmax) {
1807: for (j = jmin; j < jmax; j++) {
1808: col = bj[j];
1809: ba[j] = rtmp[col];
1810: rs += PetscAbsScalar(ba[j]);
1811: }
1812: /* add the k-th row into il and c2r */
1813: il[k] = jmin;
1814: i = bj[jmin];
1815: c2r[k] = c2r[i];
1816: c2r[i] = k;
1817: }
1819: /* MatPivotCheck() */
1820: sctx.rs = rs;
1821: sctx.pv = dk;
1822: PetscCall(MatPivotCheck(B, A, info, &sctx, i));
1823: if (sctx.newshift) break;
1824: dk = sctx.pv;
1826: ba[bdiag[k]] = 1.0 / dk; /* U(k,k) */
1827: }
1828: } while (sctx.newshift);
1830: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
1831: PetscCall(PetscFree3(rtmp, il, c2r));
1832: PetscCall(ISRestoreIndices(ip, &rip));
1833: PetscCall(ISRestoreIndices(iip, &riip));
1835: PetscCall(ISIdentity(ip, &perm_identity));
1836: if (perm_identity) {
1837: B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1838: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1839: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
1840: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;
1841: } else {
1842: B->ops->solve = MatSolve_SeqSBAIJ_1;
1843: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1;
1844: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1;
1845: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1;
1846: }
1848: C->assembled = PETSC_TRUE;
1849: C->preallocated = PETSC_TRUE;
1851: PetscCall(PetscLogFlops(C->rmap->n));
1853: /* MatPivotView() */
1854: if (sctx.nshift) {
1855: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1856: 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));
1857: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1858: PetscCall(PetscInfo(A, "number of shift_nz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1859: } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
1860: PetscCall(PetscInfo(A, "number of shift_inblocks applied %" PetscInt_FMT ", each shift_amount %g\n", sctx.nshift, (double)info->shiftamount));
1861: }
1862: }
1863: PetscFunctionReturn(PETSC_SUCCESS);
1864: }
1866: PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ_inplace(Mat B, Mat A, const MatFactorInfo *info)
1867: {
1868: Mat C = B;
1869: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1870: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)C->data;
1871: IS ip = b->row, iip = b->icol;
1872: const PetscInt *rip, *riip;
1873: PetscInt i, j, mbs = A->rmap->n, *bi = b->i, *bj = b->j, *bcol, *bjtmp;
1874: PetscInt *ai = a->i, *aj = a->j;
1875: PetscInt k, jmin, jmax, *jl, *il, col, nexti, ili, nz;
1876: MatScalar *rtmp, *ba = b->a, *bval, dk, uikdi;
1877: const MatScalar *aa, *v;
1878: PetscBool perm_identity;
1879: FactorShiftCtx sctx;
1880: PetscReal rs;
1881: MatScalar d;
1882: const PetscInt *adiag;
1884: PetscFunctionBegin;
1885: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
1886: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1887: /* MatPivotSetUp(): initialize shift context sctx */
1888: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
1890: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
1891: sctx.shift_top = info->zeropivot;
1892: for (i = 0; i < mbs; i++) {
1893: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
1894: d = aa[adiag[i]];
1895: rs = -PetscAbsScalar(d) - PetscRealPart(d);
1896: v = aa + ai[i];
1897: nz = ai[i + 1] - ai[i];
1898: for (j = 0; j < nz; j++) rs += PetscAbsScalar(v[j]);
1899: if (rs > sctx.shift_top) sctx.shift_top = rs;
1900: }
1901: sctx.shift_top *= 1.1;
1902: sctx.nshift_max = 5;
1903: sctx.shift_lo = 0.;
1904: sctx.shift_hi = 1.;
1905: }
1907: PetscCall(ISGetIndices(ip, &rip));
1908: PetscCall(ISGetIndices(iip, &riip));
1910: /* initialization */
1911: PetscCall(PetscMalloc3(mbs, &rtmp, mbs, &il, mbs, &jl));
1913: do {
1914: sctx.newshift = PETSC_FALSE;
1916: for (i = 0; i < mbs; i++) jl[i] = mbs;
1917: il[0] = 0;
1919: for (k = 0; k < mbs; k++) {
1920: /* zero rtmp */
1921: nz = bi[k + 1] - bi[k];
1922: bjtmp = bj + bi[k];
1923: for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;
1925: bval = ba + bi[k];
1926: /* initialize k-th row by the perm[k]-th row of A */
1927: jmin = ai[rip[k]];
1928: jmax = ai[rip[k] + 1];
1929: for (j = jmin; j < jmax; j++) {
1930: col = riip[aj[j]];
1931: if (col >= k) { /* only take upper triangular entry */
1932: rtmp[col] = aa[j];
1933: *bval++ = 0.0; /* for in-place factorization */
1934: }
1935: }
1936: /* shift the diagonal of the matrix */
1937: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
1939: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1940: dk = rtmp[k];
1941: i = jl[k]; /* first row to be added to k_th row */
1943: while (i < k) {
1944: nexti = jl[i]; /* next row to be added to k_th row */
1946: /* compute multiplier, update diag(k) and U(i,k) */
1947: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1948: uikdi = -ba[ili] * ba[bi[i]]; /* diagonal(k) */
1949: dk += uikdi * ba[ili];
1950: ba[ili] = uikdi; /* -U(i,k) */
1952: /* add multiple of row i to k-th row */
1953: jmin = ili + 1;
1954: jmax = bi[i + 1];
1955: if (jmin < jmax) {
1956: for (j = jmin; j < jmax; j++) rtmp[bj[j]] += uikdi * ba[j];
1957: /* update il and jl for row i */
1958: il[i] = jmin;
1959: j = bj[jmin];
1960: jl[i] = jl[j];
1961: jl[j] = i;
1962: }
1963: i = nexti;
1964: }
1966: /* shift the diagonals when zero pivot is detected */
1967: /* compute rs=sum of abs(off-diagonal) */
1968: rs = 0.0;
1969: jmin = bi[k] + 1;
1970: nz = bi[k + 1] - jmin;
1971: bcol = bj + jmin;
1972: for (j = 0; j < nz; j++) rs += PetscAbsScalar(rtmp[bcol[j]]);
1974: sctx.rs = rs;
1975: sctx.pv = dk;
1976: PetscCall(MatPivotCheck(B, A, info, &sctx, k));
1977: if (sctx.newshift) break;
1978: dk = sctx.pv;
1980: /* copy data into U(k,:) */
1981: ba[bi[k]] = 1.0 / dk; /* U(k,k) */
1982: jmin = bi[k] + 1;
1983: jmax = bi[k + 1];
1984: if (jmin < jmax) {
1985: for (j = jmin; j < jmax; j++) {
1986: col = bj[j];
1987: ba[j] = rtmp[col];
1988: }
1989: /* add the k-th row into il and jl */
1990: il[k] = jmin;
1991: i = bj[jmin];
1992: jl[k] = jl[i];
1993: jl[i] = k;
1994: }
1995: }
1996: } while (sctx.newshift);
1998: PetscCall(PetscFree3(rtmp, il, jl));
1999: PetscCall(ISRestoreIndices(ip, &rip));
2000: PetscCall(ISRestoreIndices(iip, &riip));
2001: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
2003: PetscCall(ISIdentity(ip, &perm_identity));
2004: if (perm_identity) {
2005: B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2006: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2007: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2008: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2009: } else {
2010: B->ops->solve = MatSolve_SeqSBAIJ_1_inplace;
2011: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
2012: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_inplace;
2013: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_inplace;
2014: }
2016: C->assembled = PETSC_TRUE;
2017: C->preallocated = PETSC_TRUE;
2019: PetscCall(PetscLogFlops(C->rmap->n));
2020: if (sctx.nshift) {
2021: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
2022: PetscCall(PetscInfo(A, "number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
2023: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
2024: PetscCall(PetscInfo(A, "number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
2025: }
2026: }
2027: PetscFunctionReturn(PETSC_SUCCESS);
2028: }
2030: /*
2031: icc() under revised new data structure.
2032: Factored arrays bj and ba are stored as
2033: U(0,:),...,U(i,:),U(n-1,:)
2035: ui=fact->i is an array of size n+1, in which
2036: ui+
2037: ui[i]: points to 1st entry of U(i,:),i=0,...,n-1
2038: ui[n]: points to U(n-1,n-1)+1
2040: udiag=fact->diag is an array of size n,in which
2041: udiag[i]: points to diagonal of U(i,:), i=0,...,n-1
2043: U(i,:) contains udiag[i] as its last entry, i.e.,
2044: U(i,:) = (u[i,i+1],...,u[i,n-1],diag[i])
2045: */
2047: PetscErrorCode MatICCFactorSymbolic_SeqAIJ(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
2048: {
2049: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
2050: Mat_SeqSBAIJ *b;
2051: PetscBool perm_identity;
2052: PetscInt reallocs = 0, i, *ai = a->i, *aj = a->j, am = A->rmap->n, *ui, *udiag, levels = (PetscInt)info->levels;
2053: const PetscInt *rip, *riip, *adiag;
2054: PetscInt jmin, jmax, nzk, k, j, *jl, prow, *il, nextprow;
2055: PetscInt nlnk, *lnk, *lnk_lvl = NULL;
2056: PetscInt ncols, ncols_upper, *cols, *ajtmp, *uj, **uj_ptr, **uj_lvl_ptr;
2057: PetscReal fill = info->fill;
2058: PetscFreeSpaceList free_space = NULL, current_space = NULL;
2059: PetscFreeSpaceList free_space_lvl = NULL, current_space_lvl = NULL;
2060: PetscBT lnkbt;
2061: IS iperm;
2062: PetscBool diagDense;
2064: PetscFunctionBegin;
2065: 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);
2066: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, &diagDense));
2067: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entries");
2068: PetscCall(ISIdentity(perm, &perm_identity));
2069: PetscCall(ISInvertPermutation(perm, PETSC_DECIDE, &iperm));
2071: PetscCall(PetscShmgetAllocateArray(am + 1, sizeof(PetscInt), (void **)&ui));
2072: PetscCall(PetscMalloc1(am + 1, &udiag));
2073: ui[0] = 0;
2075: /* ICC(0) without matrix ordering: simply rearrange column indices */
2076: if (!levels && perm_identity) {
2077: for (i = 0; i < am; i++) {
2078: ncols = ai[i + 1] - adiag[i];
2079: ui[i + 1] = ui[i] + ncols;
2080: udiag[i] = ui[i + 1] - 1; /* points to the last entry of U(i,:) */
2081: }
2082: PetscCall(PetscMalloc1(ui[am] + 1, &uj));
2083: cols = uj;
2084: for (i = 0; i < am; i++) {
2085: aj = a->j + adiag[i] + 1; /* 1st entry of U(i,:) without diagonal */
2086: ncols = ai[i + 1] - adiag[i] - 1;
2087: for (j = 0; j < ncols; j++) *cols++ = aj[j];
2088: *cols++ = i; /* diagonal is located as the last entry of U(i,:) */
2089: }
2090: } else { /* case: levels>0 || (levels=0 && !perm_identity) */
2091: PetscCall(ISGetIndices(iperm, &riip));
2092: PetscCall(ISGetIndices(perm, &rip));
2094: /* initialization */
2095: PetscCall(PetscMalloc1(am + 1, &ajtmp));
2097: /* jl: linked list for storing indices of the pivot rows
2098: il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2099: PetscCall(PetscMalloc4(am, &uj_ptr, am, &uj_lvl_ptr, am, &jl, am, &il));
2100: for (i = 0; i < am; i++) {
2101: jl[i] = am;
2102: il[i] = 0;
2103: }
2105: /* create and initialize a linked list for storing column indices of the active row k */
2106: nlnk = am + 1;
2107: PetscCall(PetscIncompleteLLCreate(am, am, nlnk, lnk, lnk_lvl, lnkbt));
2109: /* initial FreeSpace size is fill*(ai[am]+am)/2 */
2110: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, (ai[am] + am) / 2), &free_space));
2111: current_space = free_space;
2112: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, (ai[am] + am) / 2), &free_space_lvl));
2113: current_space_lvl = free_space_lvl;
2115: for (k = 0; k < am; k++) { /* for each active row k */
2116: /* initialize lnk by the column indices of row rip[k] of A */
2117: nzk = 0;
2118: ncols = ai[rip[k] + 1] - ai[rip[k]];
2119: PetscCheck(ncols, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Empty row in matrix: row in original ordering %" PetscInt_FMT " in permuted ordering %" PetscInt_FMT, rip[k], k);
2120: ncols_upper = 0;
2121: for (j = 0; j < ncols; j++) {
2122: i = *(aj + ai[rip[k]] + j); /* unpermuted column index */
2123: if (riip[i] >= k) { /* only take upper triangular entry */
2124: ajtmp[ncols_upper] = i;
2125: ncols_upper++;
2126: }
2127: }
2128: PetscCall(PetscIncompleteLLInit(ncols_upper, ajtmp, am, riip, &nlnk, lnk, lnk_lvl, lnkbt));
2129: nzk += nlnk;
2131: /* update lnk by computing fill-in for each pivot row to be merged in */
2132: prow = jl[k]; /* 1st pivot row */
2134: while (prow < k) {
2135: nextprow = jl[prow];
2137: /* merge prow into k-th row */
2138: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
2139: jmax = ui[prow + 1];
2140: ncols = jmax - jmin;
2141: i = jmin - ui[prow];
2142: cols = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
2143: uj = uj_lvl_ptr[prow] + i; /* levels of cols */
2144: j = *(uj - 1);
2145: PetscCall(PetscICCLLAddSorted(ncols, cols, levels, uj, am, &nlnk, lnk, lnk_lvl, lnkbt, j));
2146: nzk += nlnk;
2148: /* update il and jl for prow */
2149: if (jmin < jmax) {
2150: il[prow] = jmin;
2151: j = *cols;
2152: jl[prow] = jl[j];
2153: jl[j] = prow;
2154: }
2155: prow = nextprow;
2156: }
2158: /* if free space is not available, make more free space */
2159: if (current_space->local_remaining < nzk) {
2160: i = am - k + 1; /* num of unfactored rows */
2161: i = PetscIntMultTruncate(i, PetscMin(nzk, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
2162: PetscCall(PetscFreeSpaceGet(i, ¤t_space));
2163: PetscCall(PetscFreeSpaceGet(i, ¤t_space_lvl));
2164: reallocs++;
2165: }
2167: /* copy data into free_space and free_space_lvl, then initialize lnk */
2168: PetscCheck(nzk != 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Empty row %" PetscInt_FMT " in ICC matrix factor", k);
2169: PetscCall(PetscIncompleteLLClean(am, am, nzk, lnk, lnk_lvl, current_space->array, current_space_lvl->array, lnkbt));
2171: /* add the k-th row into il and jl */
2172: if (nzk > 1) {
2173: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
2174: jl[k] = jl[i];
2175: jl[i] = k;
2176: il[k] = ui[k] + 1;
2177: }
2178: uj_ptr[k] = current_space->array;
2179: uj_lvl_ptr[k] = current_space_lvl->array;
2181: current_space->array += nzk;
2182: current_space->local_used += nzk;
2183: current_space->local_remaining -= nzk;
2185: current_space_lvl->array += nzk;
2186: current_space_lvl->local_used += nzk;
2187: current_space_lvl->local_remaining -= nzk;
2189: ui[k + 1] = ui[k] + nzk;
2190: }
2192: PetscCall(ISRestoreIndices(perm, &rip));
2193: PetscCall(ISRestoreIndices(iperm, &riip));
2194: PetscCall(PetscFree4(uj_ptr, uj_lvl_ptr, jl, il));
2195: PetscCall(PetscFree(ajtmp));
2197: /* copy free_space into uj and free free_space; set ui, uj, udiag in new datastructure; */
2198: PetscCall(PetscShmgetAllocateArray(ui[am] + 1, sizeof(PetscInt), (void **)&uj));
2199: PetscCall(PetscFreeSpaceContiguous_Cholesky(&free_space, uj, am, ui, udiag)); /* store matrix factor */
2200: PetscCall(PetscIncompleteLLDestroy(lnk, lnkbt));
2201: PetscCall(PetscFreeSpaceDestroy(free_space_lvl));
2203: } /* end of case: levels>0 || (levels=0 && !perm_identity) */
2205: /* put together the new matrix in MATSEQSBAIJ format */
2206: b = (Mat_SeqSBAIJ *)fact->data;
2207: b->free_ij = PETSC_TRUE;
2208: PetscCall(PetscShmgetAllocateArray(ui[am], sizeof(PetscScalar), (void **)&b->a));
2209: b->free_a = PETSC_TRUE;
2210: b->j = uj;
2211: b->i = ui;
2212: b->diag = udiag;
2213: b->ilen = NULL;
2214: b->imax = NULL;
2215: b->row = perm;
2216: b->col = perm;
2217: PetscCall(PetscObjectReference((PetscObject)perm));
2218: PetscCall(PetscObjectReference((PetscObject)perm));
2219: b->icol = iperm;
2220: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
2222: PetscCall(PetscMalloc1(am, &b->solve_work));
2224: b->maxnz = b->nz = ui[am];
2226: fact->info.factor_mallocs = reallocs;
2227: fact->info.fill_ratio_given = fill;
2228: if (ai[am] != 0) {
2229: /* nonzeros in lower triangular part of A (including diagonals) = (ai[am]+am)/2 */
2230: fact->info.fill_ratio_needed = ((PetscReal)2 * ui[am]) / (ai[am] + am);
2231: } else {
2232: fact->info.fill_ratio_needed = 0.0;
2233: }
2234: #if PetscDefined(USE_INFO)
2235: if (ai[am] != 0) {
2236: PetscReal af = fact->info.fill_ratio_needed;
2237: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
2238: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
2239: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
2240: } else PetscCall(PetscInfo(A, "Empty matrix\n"));
2241: #endif
2242: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ;
2243: PetscFunctionReturn(PETSC_SUCCESS);
2244: }
2246: PetscErrorCode MatCholeskyFactorSymbolic_SeqAIJ(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
2247: {
2248: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
2249: Mat_SeqSBAIJ *b;
2250: PetscBool perm_identity;
2251: PetscReal fill = info->fill;
2252: const PetscInt *rip, *riip;
2253: PetscInt i, am = A->rmap->n, *ai = a->i, *aj = a->j, reallocs = 0, prow;
2254: PetscInt *jl, jmin, jmax, nzk, *ui, k, j, *il, nextprow;
2255: PetscInt nlnk, *lnk, ncols, ncols_upper, *cols, *uj, **ui_ptr, *uj_ptr, *udiag;
2256: PetscFreeSpaceList free_space = NULL, current_space = NULL;
2257: PetscBT lnkbt;
2258: IS iperm;
2259: PetscBool diagDense;
2261: PetscFunctionBegin;
2262: 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);
2263: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, NULL, &diagDense));
2264: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entries");
2266: /* check whether perm is the identity mapping */
2267: PetscCall(ISIdentity(perm, &perm_identity));
2268: PetscCall(ISInvertPermutation(perm, PETSC_DECIDE, &iperm));
2269: PetscCall(ISGetIndices(iperm, &riip));
2270: PetscCall(ISGetIndices(perm, &rip));
2272: /* initialization */
2273: PetscCall(PetscShmgetAllocateArray(am + 1, sizeof(PetscInt), (void **)&ui));
2274: PetscCall(PetscMalloc1(am + 1, &udiag));
2275: ui[0] = 0;
2277: /* jl: linked list for storing indices of the pivot rows
2278: il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2279: PetscCall(PetscMalloc4(am, &ui_ptr, am, &jl, am, &il, am, &cols));
2280: for (i = 0; i < am; i++) {
2281: jl[i] = am;
2282: il[i] = 0;
2283: }
2285: /* create and initialize a linked list for storing column indices of the active row k */
2286: nlnk = am + 1;
2287: PetscCall(PetscLLCreate(am, am, nlnk, lnk, lnkbt));
2289: /* initial FreeSpace size is fill*(ai[am]+am)/2 */
2290: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, (ai[am] + am) / 2), &free_space));
2291: current_space = free_space;
2293: for (k = 0; k < am; k++) { /* for each active row k */
2294: /* initialize lnk by the column indices of row rip[k] of A */
2295: nzk = 0;
2296: ncols = ai[rip[k] + 1] - ai[rip[k]];
2297: PetscCheck(ncols, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Empty row in matrix: row in original ordering %" PetscInt_FMT " in permuted ordering %" PetscInt_FMT, rip[k], k);
2298: ncols_upper = 0;
2299: for (j = 0; j < ncols; j++) {
2300: i = riip[*(aj + ai[rip[k]] + j)];
2301: if (i >= k) { /* only take upper triangular entry */
2302: cols[ncols_upper] = i;
2303: ncols_upper++;
2304: }
2305: }
2306: PetscCall(PetscLLAdd(ncols_upper, cols, am, &nlnk, lnk, lnkbt));
2307: nzk += nlnk;
2309: /* update lnk by computing fill-in for each pivot row to be merged in */
2310: prow = jl[k]; /* 1st pivot row */
2312: while (prow < k) {
2313: nextprow = jl[prow];
2314: /* merge prow into k-th row */
2315: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
2316: jmax = ui[prow + 1];
2317: ncols = jmax - jmin;
2318: uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:am-1) */
2319: PetscCall(PetscLLAddSorted(ncols, uj_ptr, am, &nlnk, lnk, lnkbt));
2320: nzk += nlnk;
2322: /* update il and jl for prow */
2323: if (jmin < jmax) {
2324: il[prow] = jmin;
2325: j = *uj_ptr;
2326: jl[prow] = jl[j];
2327: jl[j] = prow;
2328: }
2329: prow = nextprow;
2330: }
2332: /* if free space is not available, make more free space */
2333: if (current_space->local_remaining < nzk) {
2334: i = am - k + 1; /* num of unfactored rows */
2335: i = PetscIntMultTruncate(i, PetscMin(nzk, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
2336: PetscCall(PetscFreeSpaceGet(i, ¤t_space));
2337: reallocs++;
2338: }
2340: /* copy data into free space, then initialize lnk */
2341: PetscCall(PetscLLClean(am, am, nzk, lnk, current_space->array, lnkbt));
2343: /* add the k-th row into il and jl */
2344: if (nzk > 1) {
2345: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
2346: jl[k] = jl[i];
2347: jl[i] = k;
2348: il[k] = ui[k] + 1;
2349: }
2350: ui_ptr[k] = current_space->array;
2352: current_space->array += nzk;
2353: current_space->local_used += nzk;
2354: current_space->local_remaining -= nzk;
2356: ui[k + 1] = ui[k] + nzk;
2357: }
2359: PetscCall(ISRestoreIndices(perm, &rip));
2360: PetscCall(ISRestoreIndices(iperm, &riip));
2361: PetscCall(PetscFree4(ui_ptr, jl, il, cols));
2363: /* copy free_space into uj and free free_space; set ui, uj, udiag in new datastructure; */
2364: PetscCall(PetscShmgetAllocateArray(ui[am], sizeof(PetscInt), (void **)&uj));
2365: PetscCall(PetscFreeSpaceContiguous_Cholesky(&free_space, uj, am, ui, udiag)); /* store matrix factor */
2366: PetscCall(PetscLLDestroy(lnk, lnkbt));
2368: /* put together the new matrix in MATSEQSBAIJ format */
2369: b = (Mat_SeqSBAIJ *)fact->data;
2370: b->free_ij = PETSC_TRUE;
2371: PetscCall(PetscShmgetAllocateArray(ui[am], sizeof(PetscScalar), (void **)&b->a));
2372: b->free_a = PETSC_TRUE;
2373: b->j = uj;
2374: b->i = ui;
2375: b->diag = udiag;
2376: b->ilen = NULL;
2377: b->imax = NULL;
2378: b->row = perm;
2379: b->col = perm;
2381: PetscCall(PetscObjectReference((PetscObject)perm));
2382: PetscCall(PetscObjectReference((PetscObject)perm));
2384: b->icol = iperm;
2385: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
2387: PetscCall(PetscMalloc1(am, &b->solve_work));
2389: b->maxnz = b->nz = ui[am];
2391: fact->info.factor_mallocs = reallocs;
2392: fact->info.fill_ratio_given = fill;
2393: if (ai[am] != 0) {
2394: /* nonzeros in lower triangular part of A (including diagonals) = (ai[am]+am)/2 */
2395: fact->info.fill_ratio_needed = ((PetscReal)2 * ui[am]) / (ai[am] + am);
2396: } else {
2397: fact->info.fill_ratio_needed = 0.0;
2398: }
2399: #if PetscDefined(USE_INFO)
2400: if (ai[am] != 0) {
2401: PetscReal af = fact->info.fill_ratio_needed;
2402: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
2403: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
2404: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
2405: } else PetscCall(PetscInfo(A, "Empty matrix\n"));
2406: #endif
2407: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ;
2408: PetscFunctionReturn(PETSC_SUCCESS);
2409: }
2411: PetscErrorCode MatSolve_SeqAIJ_NaturalOrdering(Mat A, Vec bb, Vec xx)
2412: {
2413: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
2414: PetscInt n = A->rmap->n;
2415: const PetscInt *ai = a->i, *aj = a->j, *adiag = a->diag, *vi;
2416: PetscScalar *x, sum;
2417: const PetscScalar *b;
2418: const MatScalar *aa, *v;
2419: PetscInt i, nz;
2421: PetscFunctionBegin;
2422: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
2424: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
2425: PetscCall(VecGetArrayRead(bb, &b));
2426: PetscCall(VecGetArrayWrite(xx, &x));
2428: /* forward solve the lower triangular */
2429: x[0] = b[0];
2430: v = aa;
2431: vi = aj;
2432: for (i = 1; i < n; i++) {
2433: nz = ai[i + 1] - ai[i];
2434: sum = b[i];
2435: PetscSparseDenseMinusDot(sum, x, v, vi, nz);
2436: v += nz;
2437: vi += nz;
2438: x[i] = sum;
2439: }
2441: /* backward solve the upper triangular */
2442: for (i = n - 1; i >= 0; i--) {
2443: v = aa + adiag[i + 1] + 1;
2444: vi = aj + adiag[i + 1] + 1;
2445: nz = adiag[i] - adiag[i + 1] - 1;
2446: sum = x[i];
2447: PetscSparseDenseMinusDot(sum, x, v, vi, nz);
2448: x[i] = sum * v[nz]; /* x[i]=aa[adiag[i]]*sum; v++; */
2449: }
2451: PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
2452: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
2453: PetscCall(VecRestoreArrayRead(bb, &b));
2454: PetscCall(VecRestoreArrayWrite(xx, &x));
2455: PetscFunctionReturn(PETSC_SUCCESS);
2456: }
2458: PetscErrorCode MatSolve_SeqAIJ(Mat A, Vec bb, Vec xx)
2459: {
2460: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
2461: IS iscol = a->col, isrow = a->row;
2462: PetscInt i, n = A->rmap->n, *vi, *ai = a->i, *aj = a->j, *adiag = a->diag, nz;
2463: const PetscInt *rout, *cout, *r, *c;
2464: PetscScalar *x, *tmp, sum;
2465: const PetscScalar *b;
2466: const MatScalar *aa, *v;
2468: PetscFunctionBegin;
2469: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
2471: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
2472: PetscCall(VecGetArrayRead(bb, &b));
2473: PetscCall(VecGetArrayWrite(xx, &x));
2474: tmp = a->solve_work;
2476: PetscCall(ISGetIndices(isrow, &rout));
2477: r = rout;
2478: PetscCall(ISGetIndices(iscol, &cout));
2479: c = cout;
2481: /* forward solve the lower triangular */
2482: tmp[0] = b[r[0]];
2483: v = aa;
2484: vi = aj;
2485: for (i = 1; i < n; i++) {
2486: nz = ai[i + 1] - ai[i];
2487: sum = b[r[i]];
2488: PetscSparseDenseMinusDot(sum, tmp, v, vi, nz);
2489: tmp[i] = sum;
2490: v += nz;
2491: vi += nz;
2492: }
2494: /* backward solve the upper triangular */
2495: for (i = n - 1; i >= 0; i--) {
2496: v = aa + adiag[i + 1] + 1;
2497: vi = aj + adiag[i + 1] + 1;
2498: nz = adiag[i] - adiag[i + 1] - 1;
2499: sum = tmp[i];
2500: PetscSparseDenseMinusDot(sum, tmp, v, vi, nz);
2501: x[c[i]] = tmp[i] = sum * v[nz]; /* v[nz] = aa[adiag[i]] */
2502: }
2504: PetscCall(ISRestoreIndices(isrow, &rout));
2505: PetscCall(ISRestoreIndices(iscol, &cout));
2506: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
2507: PetscCall(VecRestoreArrayRead(bb, &b));
2508: PetscCall(VecRestoreArrayWrite(xx, &x));
2509: PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
2510: PetscFunctionReturn(PETSC_SUCCESS);
2511: }