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