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, 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: levels = (PetscInt)info->levels;
1560: PetscCall(ISIdentity(isrow, &row_identity));
1561: PetscCall(ISIdentity(iscol, &col_identity));
1562: if (!levels && row_identity && col_identity) {
1563: /* special case: ilu(0) with natural ordering */
1564: PetscCall(MatILUFactorSymbolic_SeqAIJ_ilu0(fact, A, isrow, iscol, info));
1565: if (a->inode.size_csr) fact->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
1566: PetscFunctionReturn(PETSC_SUCCESS);
1567: }
1569: PetscCall(ISInvertPermutation(iscol, PETSC_DECIDE, &isicol));
1570: PetscCall(ISGetIndices(isrow, &r));
1571: PetscCall(ISGetIndices(isicol, &ic));
1573: /* get new row and diagonal pointers, must be allocated separately because they will be given to the Mat_SeqAIJ and freed separately */
1574: PetscCall(PetscShmgetAllocateArray(n + 1, sizeof(PetscInt), (void **)&bi));
1575: PetscCall(PetscMalloc1(n + 1, &bdiag));
1576: bi[0] = bdiag[0] = 0;
1577: PetscCall(PetscMalloc2(n, &bj_ptr, n, &bjlvl_ptr));
1579: /* create a linked list for storing column indices of the active row */
1580: nlnk = n + 1;
1581: PetscCall(PetscIncompleteLLCreate(n, n, nlnk, lnk, lnk_lvl, lnkbt));
1583: /* initial FreeSpace size is f*(ai[n]+1) */
1584: f = info->fill;
1585: diagonal_fill = (PetscInt)info->diagonal_fill;
1586: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space));
1587: current_space = free_space;
1588: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space_lvl));
1589: current_space_lvl = free_space_lvl;
1590: for (i = 0; i < n; i++) {
1591: nzi = 0;
1592: /* copy current row into linked list */
1593: nnz = ai[r[i] + 1] - ai[r[i]];
1594: 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);
1595: cols = aj + ai[r[i]];
1596: lnk[i] = -1; /* marker to indicate if diagonal exists */
1597: PetscCall(PetscIncompleteLLInit(nnz, cols, n, ic, &nlnk, lnk, lnk_lvl, lnkbt));
1598: nzi += nlnk;
1600: /* make sure diagonal entry is included */
1601: if (diagonal_fill && lnk[i] == -1) {
1602: fm = n;
1603: while (lnk[fm] < i) fm = lnk[fm];
1604: lnk[i] = lnk[fm]; /* insert diagonal into linked list */
1605: lnk[fm] = i;
1606: lnk_lvl[i] = 0;
1607: nzi++;
1608: dcount++;
1609: }
1611: /* add pivot rows into the active row */
1612: nzbd = 0;
1613: prow = lnk[n];
1614: while (prow < i) {
1615: nnz = bdiag[prow];
1616: cols = bj_ptr[prow] + nnz + 1;
1617: cols_lvl = bjlvl_ptr[prow] + nnz + 1;
1618: nnz = bi[prow + 1] - bi[prow] - nnz - 1;
1619: PetscCall(PetscILULLAddSorted(nnz, cols, levels, cols_lvl, prow, &nlnk, lnk, lnk_lvl, lnkbt, prow));
1620: nzi += nlnk;
1621: prow = lnk[prow];
1622: nzbd++;
1623: }
1624: bdiag[i] = nzbd;
1625: bi[i + 1] = bi[i] + nzi;
1626: /* if free space is not available, make more free space */
1627: if (current_space->local_remaining < nzi) {
1628: nnz = PetscIntMultTruncate(2, PetscIntMultTruncate(nzi, n - i)); /* estimated and max additional space needed */
1629: PetscCall(PetscFreeSpaceGet(nnz, ¤t_space));
1630: PetscCall(PetscFreeSpaceGet(nnz, ¤t_space_lvl));
1631: reallocs++;
1632: }
1634: /* copy data into free_space and free_space_lvl, then initialize lnk */
1635: PetscCall(PetscIncompleteLLClean(n, n, nzi, lnk, lnk_lvl, current_space->array, current_space_lvl->array, lnkbt));
1636: bj_ptr[i] = current_space->array;
1637: bjlvl_ptr[i] = current_space_lvl->array;
1639: /* make sure the active row i has diagonal entry */
1640: 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);
1642: current_space->array += nzi;
1643: current_space->local_used += nzi;
1644: current_space->local_remaining -= nzi;
1645: current_space_lvl->array += nzi;
1646: current_space_lvl->local_used += nzi;
1647: current_space_lvl->local_remaining -= nzi;
1648: }
1650: PetscCall(ISRestoreIndices(isrow, &r));
1651: PetscCall(ISRestoreIndices(isicol, &ic));
1652: /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
1653: PetscCall(PetscShmgetAllocateArray(bi[n], sizeof(PetscInt), (void **)&bj));
1654: PetscCall(PetscFreeSpaceContiguous_LU(&free_space, bj, n, bi, bdiag));
1656: PetscCall(PetscIncompleteLLDestroy(lnk, lnkbt));
1657: PetscCall(PetscFreeSpaceDestroy(free_space_lvl));
1658: PetscCall(PetscFree2(bj_ptr, bjlvl_ptr));
1660: #if PetscDefined(USE_INFO)
1661: {
1662: PetscReal af = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);
1663: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)f, (double)af));
1664: PetscCall(PetscInfo(A, "Run with -[sub_]pc_factor_fill %g or use \n", (double)af));
1665: PetscCall(PetscInfo(A, "PCFactorSetFill([sub]pc,%g);\n", (double)af));
1666: PetscCall(PetscInfo(A, "for best performance.\n"));
1667: if (diagonal_fill) PetscCall(PetscInfo(A, "Detected and replaced %" PetscInt_FMT " missing diagonals\n", dcount));
1668: }
1669: #endif
1670: /* put together the new matrix */
1671: PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(fact, MAT_SKIP_ALLOCATION, NULL));
1672: b = (Mat_SeqAIJ *)fact->data;
1673: b->free_ij = PETSC_TRUE;
1674: PetscCall(PetscShmgetAllocateArray(bdiag[0] + 1, sizeof(PetscScalar), (void **)&b->a));
1675: b->free_a = PETSC_TRUE;
1676: b->j = bj;
1677: b->i = bi;
1678: b->diag = bdiag;
1679: b->ilen = NULL;
1680: b->imax = NULL;
1681: b->row = isrow;
1682: b->col = iscol;
1683: PetscCall(PetscObjectReference((PetscObject)isrow));
1684: PetscCall(PetscObjectReference((PetscObject)iscol));
1685: b->icol = isicol;
1687: PetscCall(PetscMalloc1(n, &b->solve_work));
1688: /* In b structure: Free imax, ilen, old a, old j.
1689: Allocate bdiag, solve_work, new a, new j */
1690: b->maxnz = b->nz = bdiag[0] + 1;
1692: fact->info.factor_mallocs = reallocs;
1693: fact->info.fill_ratio_given = f;
1694: fact->info.fill_ratio_needed = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);
1695: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ;
1696: if (a->inode.size_csr) fact->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
1697: PetscCall(MatSeqAIJCheckInode_FactorLU(fact));
1698: PetscFunctionReturn(PETSC_SUCCESS);
1699: }
1701: PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ(Mat B, Mat A, const MatFactorInfo *info)
1702: {
1703: Mat C = B;
1704: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1705: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)C->data;
1706: IS ip = b->row, iip = b->icol;
1707: const PetscInt *rip, *riip;
1708: PetscInt i, j, mbs = A->rmap->n, *bi = b->i, *bj = b->j, *bdiag = b->diag, *bjtmp;
1709: PetscInt *ai = a->i, *aj = a->j;
1710: PetscInt k, jmin, jmax, *c2r, *il, col, nexti, ili, nz;
1711: MatScalar *rtmp, *ba = b->a, *bval, dk, uikdi;
1712: PetscBool perm_identity;
1713: FactorShiftCtx sctx;
1714: PetscReal rs;
1715: const MatScalar *aa, *v;
1716: MatScalar d;
1717: const PetscInt *adiag;
1719: PetscFunctionBegin;
1720: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
1721: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1722: /* MatPivotSetUp(): initialize shift context sctx */
1723: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
1725: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
1726: sctx.shift_top = info->zeropivot;
1727: for (i = 0; i < mbs; i++) {
1728: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
1729: d = aa[adiag[i]];
1730: rs = -PetscAbsScalar(d) - PetscRealPart(d);
1731: v = aa + ai[i];
1732: nz = ai[i + 1] - ai[i];
1733: for (j = 0; j < nz; j++) rs += PetscAbsScalar(v[j]);
1734: if (rs > sctx.shift_top) sctx.shift_top = rs;
1735: }
1736: sctx.shift_top *= 1.1;
1737: sctx.nshift_max = 5;
1738: sctx.shift_lo = 0.;
1739: sctx.shift_hi = 1.;
1740: }
1742: PetscCall(ISGetIndices(ip, &rip));
1743: PetscCall(ISGetIndices(iip, &riip));
1745: /* allocate working arrays
1746: c2r: linked list, keep track of pivot rows for a given column. c2r[col]: head of the list for a given col
1747: 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
1748: */
1749: PetscCall(PetscMalloc3(mbs, &rtmp, mbs, &il, mbs, &c2r));
1751: do {
1752: sctx.newshift = PETSC_FALSE;
1754: for (i = 0; i < mbs; i++) c2r[i] = mbs;
1755: if (mbs) il[0] = 0;
1757: for (k = 0; k < mbs; k++) {
1758: /* zero rtmp */
1759: nz = bi[k + 1] - bi[k];
1760: bjtmp = bj + bi[k];
1761: for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;
1763: /* load in initial unfactored row */
1764: bval = ba + bi[k];
1765: jmin = ai[rip[k]];
1766: jmax = ai[rip[k] + 1];
1767: for (j = jmin; j < jmax; j++) {
1768: col = riip[aj[j]];
1769: if (col >= k) { /* only take upper triangular entry */
1770: rtmp[col] = aa[j];
1771: *bval++ = 0.0; /* for in-place factorization */
1772: }
1773: }
1774: /* shift the diagonal of the matrix: ZeropivotApply() */
1775: rtmp[k] += sctx.shift_amount; /* shift the diagonal of the matrix */
1777: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1778: dk = rtmp[k];
1779: i = c2r[k]; /* first row to be added to k_th row */
1781: while (i < k) {
1782: nexti = c2r[i]; /* next row to be added to k_th row */
1784: /* compute multiplier, update diag(k) and U(i,k) */
1785: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1786: uikdi = -ba[ili] * ba[bdiag[i]]; /* diagonal(k) */
1787: dk += uikdi * ba[ili]; /* update diag[k] */
1788: ba[ili] = uikdi; /* -U(i,k) */
1790: /* add multiple of row i to k-th row */
1791: jmin = ili + 1;
1792: jmax = bi[i + 1];
1793: if (jmin < jmax) {
1794: for (j = jmin; j < jmax; j++) rtmp[bj[j]] += uikdi * ba[j];
1795: /* update il and c2r for row i */
1796: il[i] = jmin;
1797: j = bj[jmin];
1798: c2r[i] = c2r[j];
1799: c2r[j] = i;
1800: }
1801: i = nexti;
1802: }
1804: /* copy data into U(k,:) */
1805: rs = 0.0;
1806: jmin = bi[k];
1807: jmax = bi[k + 1] - 1;
1808: if (jmin < jmax) {
1809: for (j = jmin; j < jmax; j++) {
1810: col = bj[j];
1811: ba[j] = rtmp[col];
1812: rs += PetscAbsScalar(ba[j]);
1813: }
1814: /* add the k-th row into il and c2r */
1815: il[k] = jmin;
1816: i = bj[jmin];
1817: c2r[k] = c2r[i];
1818: c2r[i] = k;
1819: }
1821: /* MatPivotCheck() */
1822: sctx.rs = rs;
1823: sctx.pv = dk;
1824: PetscCall(MatPivotCheck(B, A, info, &sctx, i));
1825: if (sctx.newshift) break;
1826: dk = sctx.pv;
1828: ba[bdiag[k]] = 1.0 / dk; /* U(k,k) */
1829: }
1830: } while (sctx.newshift);
1832: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
1833: PetscCall(PetscFree3(rtmp, il, c2r));
1834: PetscCall(ISRestoreIndices(ip, &rip));
1835: PetscCall(ISRestoreIndices(iip, &riip));
1837: PetscCall(ISIdentity(ip, &perm_identity));
1838: if (perm_identity) {
1839: B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1840: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1841: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
1842: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;
1843: } else {
1844: B->ops->solve = MatSolve_SeqSBAIJ_1;
1845: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1;
1846: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1;
1847: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1;
1848: }
1850: C->assembled = PETSC_TRUE;
1851: C->preallocated = PETSC_TRUE;
1853: PetscCall(PetscLogFlops(C->rmap->n));
1855: /* MatPivotView() */
1856: if (sctx.nshift) {
1857: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1858: 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));
1859: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1860: PetscCall(PetscInfo(A, "number of shift_nz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1861: } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
1862: PetscCall(PetscInfo(A, "number of shift_inblocks applied %" PetscInt_FMT ", each shift_amount %g\n", sctx.nshift, (double)info->shiftamount));
1863: }
1864: }
1865: PetscFunctionReturn(PETSC_SUCCESS);
1866: }
1868: PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ_inplace(Mat B, Mat A, const MatFactorInfo *info)
1869: {
1870: Mat C = B;
1871: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1872: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)C->data;
1873: IS ip = b->row, iip = b->icol;
1874: const PetscInt *rip, *riip;
1875: PetscInt i, j, mbs = A->rmap->n, *bi = b->i, *bj = b->j, *bcol, *bjtmp;
1876: PetscInt *ai = a->i, *aj = a->j;
1877: PetscInt k, jmin, jmax, *jl, *il, col, nexti, ili, nz;
1878: MatScalar *rtmp, *ba = b->a, *bval, dk, uikdi;
1879: const MatScalar *aa, *v;
1880: PetscBool perm_identity;
1881: FactorShiftCtx sctx;
1882: PetscReal rs;
1883: MatScalar d;
1884: const PetscInt *adiag;
1886: PetscFunctionBegin;
1887: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
1888: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
1889: /* MatPivotSetUp(): initialize shift context sctx */
1890: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
1892: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
1893: sctx.shift_top = info->zeropivot;
1894: for (i = 0; i < mbs; i++) {
1895: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
1896: d = aa[adiag[i]];
1897: rs = -PetscAbsScalar(d) - PetscRealPart(d);
1898: v = aa + ai[i];
1899: nz = ai[i + 1] - ai[i];
1900: for (j = 0; j < nz; j++) rs += PetscAbsScalar(v[j]);
1901: if (rs > sctx.shift_top) sctx.shift_top = rs;
1902: }
1903: sctx.shift_top *= 1.1;
1904: sctx.nshift_max = 5;
1905: sctx.shift_lo = 0.;
1906: sctx.shift_hi = 1.;
1907: }
1909: PetscCall(ISGetIndices(ip, &rip));
1910: PetscCall(ISGetIndices(iip, &riip));
1912: /* initialization */
1913: PetscCall(PetscMalloc3(mbs, &rtmp, mbs, &il, mbs, &jl));
1915: do {
1916: sctx.newshift = PETSC_FALSE;
1918: for (i = 0; i < mbs; i++) jl[i] = mbs;
1919: il[0] = 0;
1921: for (k = 0; k < mbs; k++) {
1922: /* zero rtmp */
1923: nz = bi[k + 1] - bi[k];
1924: bjtmp = bj + bi[k];
1925: for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;
1927: bval = ba + bi[k];
1928: /* initialize k-th row by the perm[k]-th row of A */
1929: jmin = ai[rip[k]];
1930: jmax = ai[rip[k] + 1];
1931: for (j = jmin; j < jmax; j++) {
1932: col = riip[aj[j]];
1933: if (col >= k) { /* only take upper triangular entry */
1934: rtmp[col] = aa[j];
1935: *bval++ = 0.0; /* for in-place factorization */
1936: }
1937: }
1938: /* shift the diagonal of the matrix */
1939: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
1941: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1942: dk = rtmp[k];
1943: i = jl[k]; /* first row to be added to k_th row */
1945: while (i < k) {
1946: nexti = jl[i]; /* next row to be added to k_th row */
1948: /* compute multiplier, update diag(k) and U(i,k) */
1949: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1950: uikdi = -ba[ili] * ba[bi[i]]; /* diagonal(k) */
1951: dk += uikdi * ba[ili];
1952: ba[ili] = uikdi; /* -U(i,k) */
1954: /* add multiple of row i to k-th row */
1955: jmin = ili + 1;
1956: jmax = bi[i + 1];
1957: if (jmin < jmax) {
1958: for (j = jmin; j < jmax; j++) rtmp[bj[j]] += uikdi * ba[j];
1959: /* update il and jl for row i */
1960: il[i] = jmin;
1961: j = bj[jmin];
1962: jl[i] = jl[j];
1963: jl[j] = i;
1964: }
1965: i = nexti;
1966: }
1968: /* shift the diagonals when zero pivot is detected */
1969: /* compute rs=sum of abs(off-diagonal) */
1970: rs = 0.0;
1971: jmin = bi[k] + 1;
1972: nz = bi[k + 1] - jmin;
1973: bcol = bj + jmin;
1974: for (j = 0; j < nz; j++) rs += PetscAbsScalar(rtmp[bcol[j]]);
1976: sctx.rs = rs;
1977: sctx.pv = dk;
1978: PetscCall(MatPivotCheck(B, A, info, &sctx, k));
1979: if (sctx.newshift) break;
1980: dk = sctx.pv;
1982: /* copy data into U(k,:) */
1983: ba[bi[k]] = 1.0 / dk; /* U(k,k) */
1984: jmin = bi[k] + 1;
1985: jmax = bi[k + 1];
1986: if (jmin < jmax) {
1987: for (j = jmin; j < jmax; j++) {
1988: col = bj[j];
1989: ba[j] = rtmp[col];
1990: }
1991: /* add the k-th row into il and jl */
1992: il[k] = jmin;
1993: i = bj[jmin];
1994: jl[k] = jl[i];
1995: jl[i] = k;
1996: }
1997: }
1998: } while (sctx.newshift);
2000: PetscCall(PetscFree3(rtmp, il, jl));
2001: PetscCall(ISRestoreIndices(ip, &rip));
2002: PetscCall(ISRestoreIndices(iip, &riip));
2003: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
2005: PetscCall(ISIdentity(ip, &perm_identity));
2006: if (perm_identity) {
2007: B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2008: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2009: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2010: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2011: } else {
2012: B->ops->solve = MatSolve_SeqSBAIJ_1_inplace;
2013: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
2014: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_inplace;
2015: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_inplace;
2016: }
2018: C->assembled = PETSC_TRUE;
2019: C->preallocated = PETSC_TRUE;
2021: PetscCall(PetscLogFlops(C->rmap->n));
2022: if (sctx.nshift) {
2023: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
2024: PetscCall(PetscInfo(A, "number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
2025: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
2026: PetscCall(PetscInfo(A, "number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
2027: }
2028: }
2029: PetscFunctionReturn(PETSC_SUCCESS);
2030: }
2032: /*
2033: icc() under revised new data structure.
2034: Factored arrays bj and ba are stored as
2035: U(0,:),...,U(i,:),U(n-1,:)
2037: ui=fact->i is an array of size n+1, in which
2038: ui+
2039: ui[i]: points to 1st entry of U(i,:),i=0,...,n-1
2040: ui[n]: points to U(n-1,n-1)+1
2042: udiag=fact->diag is an array of size n,in which
2043: udiag[i]: points to diagonal of U(i,:), i=0,...,n-1
2045: U(i,:) contains udiag[i] as its last entry, i.e.,
2046: U(i,:) = (u[i,i+1],...,u[i,n-1],diag[i])
2047: */
2049: PetscErrorCode MatICCFactorSymbolic_SeqAIJ(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
2050: {
2051: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
2052: Mat_SeqSBAIJ *b;
2053: PetscBool perm_identity;
2054: PetscInt reallocs = 0, i, *ai = a->i, *aj = a->j, am = A->rmap->n, *ui, *udiag;
2055: const PetscInt *rip, *riip, *adiag;
2056: PetscInt jmin, jmax, nzk, k, j, *jl, prow, *il, nextprow;
2057: PetscInt nlnk, *lnk, *lnk_lvl = NULL;
2058: PetscInt ncols, ncols_upper, *cols, *ajtmp, *uj, **uj_ptr, **uj_lvl_ptr;
2059: PetscReal fill = info->fill, levels = info->levels;
2060: PetscFreeSpaceList free_space = NULL, current_space = NULL;
2061: PetscFreeSpaceList free_space_lvl = NULL, current_space_lvl = NULL;
2062: PetscBT lnkbt;
2063: IS iperm;
2064: PetscBool diagDense;
2066: PetscFunctionBegin;
2067: 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);
2068: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, &diagDense));
2069: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entries");
2070: PetscCall(ISIdentity(perm, &perm_identity));
2071: PetscCall(ISInvertPermutation(perm, PETSC_DECIDE, &iperm));
2073: PetscCall(PetscShmgetAllocateArray(am + 1, sizeof(PetscInt), (void **)&ui));
2074: PetscCall(PetscMalloc1(am + 1, &udiag));
2075: ui[0] = 0;
2077: /* ICC(0) without matrix ordering: simply rearrange column indices */
2078: if (!levels && perm_identity) {
2079: for (i = 0; i < am; i++) {
2080: ncols = ai[i + 1] - adiag[i];
2081: ui[i + 1] = ui[i] + ncols;
2082: udiag[i] = ui[i + 1] - 1; /* points to the last entry of U(i,:) */
2083: }
2084: PetscCall(PetscMalloc1(ui[am] + 1, &uj));
2085: cols = uj;
2086: for (i = 0; i < am; i++) {
2087: aj = a->j + adiag[i] + 1; /* 1st entry of U(i,:) without diagonal */
2088: ncols = ai[i + 1] - adiag[i] - 1;
2089: for (j = 0; j < ncols; j++) *cols++ = aj[j];
2090: *cols++ = i; /* diagonal is located as the last entry of U(i,:) */
2091: }
2092: } else { /* case: levels>0 || (levels=0 && !perm_identity) */
2093: PetscCall(ISGetIndices(iperm, &riip));
2094: PetscCall(ISGetIndices(perm, &rip));
2096: /* initialization */
2097: PetscCall(PetscMalloc1(am + 1, &ajtmp));
2099: /* jl: linked list for storing indices of the pivot rows
2100: il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2101: PetscCall(PetscMalloc4(am, &uj_ptr, am, &uj_lvl_ptr, am, &jl, am, &il));
2102: for (i = 0; i < am; i++) {
2103: jl[i] = am;
2104: il[i] = 0;
2105: }
2107: /* create and initialize a linked list for storing column indices of the active row k */
2108: nlnk = am + 1;
2109: PetscCall(PetscIncompleteLLCreate(am, am, nlnk, lnk, lnk_lvl, lnkbt));
2111: /* initial FreeSpace size is fill*(ai[am]+am)/2 */
2112: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, (ai[am] + am) / 2), &free_space));
2113: current_space = free_space;
2114: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, (ai[am] + am) / 2), &free_space_lvl));
2115: current_space_lvl = free_space_lvl;
2117: for (k = 0; k < am; k++) { /* for each active row k */
2118: /* initialize lnk by the column indices of row rip[k] of A */
2119: nzk = 0;
2120: ncols = ai[rip[k] + 1] - ai[rip[k]];
2121: 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);
2122: ncols_upper = 0;
2123: for (j = 0; j < ncols; j++) {
2124: i = *(aj + ai[rip[k]] + j); /* unpermuted column index */
2125: if (riip[i] >= k) { /* only take upper triangular entry */
2126: ajtmp[ncols_upper] = i;
2127: ncols_upper++;
2128: }
2129: }
2130: PetscCall(PetscIncompleteLLInit(ncols_upper, ajtmp, am, riip, &nlnk, lnk, lnk_lvl, lnkbt));
2131: nzk += nlnk;
2133: /* update lnk by computing fill-in for each pivot row to be merged in */
2134: prow = jl[k]; /* 1st pivot row */
2136: while (prow < k) {
2137: nextprow = jl[prow];
2139: /* merge prow into k-th row */
2140: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
2141: jmax = ui[prow + 1];
2142: ncols = jmax - jmin;
2143: i = jmin - ui[prow];
2144: cols = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
2145: uj = uj_lvl_ptr[prow] + i; /* levels of cols */
2146: j = *(uj - 1);
2147: PetscCall(PetscICCLLAddSorted(ncols, cols, levels, uj, am, &nlnk, lnk, lnk_lvl, lnkbt, j));
2148: nzk += nlnk;
2150: /* update il and jl for prow */
2151: if (jmin < jmax) {
2152: il[prow] = jmin;
2153: j = *cols;
2154: jl[prow] = jl[j];
2155: jl[j] = prow;
2156: }
2157: prow = nextprow;
2158: }
2160: /* if free space is not available, make more free space */
2161: if (current_space->local_remaining < nzk) {
2162: i = am - k + 1; /* num of unfactored rows */
2163: i = PetscIntMultTruncate(i, PetscMin(nzk, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
2164: PetscCall(PetscFreeSpaceGet(i, ¤t_space));
2165: PetscCall(PetscFreeSpaceGet(i, ¤t_space_lvl));
2166: reallocs++;
2167: }
2169: /* copy data into free_space and free_space_lvl, then initialize lnk */
2170: PetscCheck(nzk != 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Empty row %" PetscInt_FMT " in ICC matrix factor", k);
2171: PetscCall(PetscIncompleteLLClean(am, am, nzk, lnk, lnk_lvl, current_space->array, current_space_lvl->array, lnkbt));
2173: /* add the k-th row into il and jl */
2174: if (nzk > 1) {
2175: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
2176: jl[k] = jl[i];
2177: jl[i] = k;
2178: il[k] = ui[k] + 1;
2179: }
2180: uj_ptr[k] = current_space->array;
2181: uj_lvl_ptr[k] = current_space_lvl->array;
2183: current_space->array += nzk;
2184: current_space->local_used += nzk;
2185: current_space->local_remaining -= nzk;
2187: current_space_lvl->array += nzk;
2188: current_space_lvl->local_used += nzk;
2189: current_space_lvl->local_remaining -= nzk;
2191: ui[k + 1] = ui[k] + nzk;
2192: }
2194: PetscCall(ISRestoreIndices(perm, &rip));
2195: PetscCall(ISRestoreIndices(iperm, &riip));
2196: PetscCall(PetscFree4(uj_ptr, uj_lvl_ptr, jl, il));
2197: PetscCall(PetscFree(ajtmp));
2199: /* copy free_space into uj and free free_space; set ui, uj, udiag in new datastructure; */
2200: PetscCall(PetscShmgetAllocateArray(ui[am] + 1, sizeof(PetscInt), (void **)&uj));
2201: PetscCall(PetscFreeSpaceContiguous_Cholesky(&free_space, uj, am, ui, udiag)); /* store matrix factor */
2202: PetscCall(PetscIncompleteLLDestroy(lnk, lnkbt));
2203: PetscCall(PetscFreeSpaceDestroy(free_space_lvl));
2205: } /* end of case: levels>0 || (levels=0 && !perm_identity) */
2207: /* put together the new matrix in MATSEQSBAIJ format */
2208: b = (Mat_SeqSBAIJ *)fact->data;
2209: b->free_ij = PETSC_TRUE;
2210: PetscCall(PetscShmgetAllocateArray(ui[am], sizeof(PetscScalar), (void **)&b->a));
2211: b->free_a = PETSC_TRUE;
2212: b->j = uj;
2213: b->i = ui;
2214: b->diag = udiag;
2215: b->ilen = NULL;
2216: b->imax = NULL;
2217: b->row = perm;
2218: b->col = perm;
2219: PetscCall(PetscObjectReference((PetscObject)perm));
2220: PetscCall(PetscObjectReference((PetscObject)perm));
2221: b->icol = iperm;
2222: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
2224: PetscCall(PetscMalloc1(am, &b->solve_work));
2226: b->maxnz = b->nz = ui[am];
2228: fact->info.factor_mallocs = reallocs;
2229: fact->info.fill_ratio_given = fill;
2230: if (ai[am] != 0) {
2231: /* nonzeros in lower triangular part of A (including diagonals) = (ai[am]+am)/2 */
2232: fact->info.fill_ratio_needed = ((PetscReal)2 * ui[am]) / (ai[am] + am);
2233: } else {
2234: fact->info.fill_ratio_needed = 0.0;
2235: }
2236: #if PetscDefined(USE_INFO)
2237: if (ai[am] != 0) {
2238: PetscReal af = fact->info.fill_ratio_needed;
2239: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
2240: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
2241: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
2242: } else PetscCall(PetscInfo(A, "Empty matrix\n"));
2243: #endif
2244: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ;
2245: PetscFunctionReturn(PETSC_SUCCESS);
2246: }
2248: PetscErrorCode MatCholeskyFactorSymbolic_SeqAIJ(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
2249: {
2250: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
2251: Mat_SeqSBAIJ *b;
2252: PetscBool perm_identity;
2253: PetscReal fill = info->fill;
2254: const PetscInt *rip, *riip;
2255: PetscInt i, am = A->rmap->n, *ai = a->i, *aj = a->j, reallocs = 0, prow;
2256: PetscInt *jl, jmin, jmax, nzk, *ui, k, j, *il, nextprow;
2257: PetscInt nlnk, *lnk, ncols, ncols_upper, *cols, *uj, **ui_ptr, *uj_ptr, *udiag;
2258: PetscFreeSpaceList free_space = NULL, current_space = NULL;
2259: PetscBT lnkbt;
2260: IS iperm;
2261: PetscBool diagDense;
2263: PetscFunctionBegin;
2264: 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);
2265: PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, NULL, &diagDense));
2266: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entries");
2268: /* check whether perm is the identity mapping */
2269: PetscCall(ISIdentity(perm, &perm_identity));
2270: PetscCall(ISInvertPermutation(perm, PETSC_DECIDE, &iperm));
2271: PetscCall(ISGetIndices(iperm, &riip));
2272: PetscCall(ISGetIndices(perm, &rip));
2274: /* initialization */
2275: PetscCall(PetscShmgetAllocateArray(am + 1, sizeof(PetscInt), (void **)&ui));
2276: PetscCall(PetscMalloc1(am + 1, &udiag));
2277: ui[0] = 0;
2279: /* jl: linked list for storing indices of the pivot rows
2280: il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2281: PetscCall(PetscMalloc4(am, &ui_ptr, am, &jl, am, &il, am, &cols));
2282: for (i = 0; i < am; i++) {
2283: jl[i] = am;
2284: il[i] = 0;
2285: }
2287: /* create and initialize a linked list for storing column indices of the active row k */
2288: nlnk = am + 1;
2289: PetscCall(PetscLLCreate(am, am, nlnk, lnk, lnkbt));
2291: /* initial FreeSpace size is fill*(ai[am]+am)/2 */
2292: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, (ai[am] + am) / 2), &free_space));
2293: current_space = free_space;
2295: for (k = 0; k < am; k++) { /* for each active row k */
2296: /* initialize lnk by the column indices of row rip[k] of A */
2297: nzk = 0;
2298: ncols = ai[rip[k] + 1] - ai[rip[k]];
2299: 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);
2300: ncols_upper = 0;
2301: for (j = 0; j < ncols; j++) {
2302: i = riip[*(aj + ai[rip[k]] + j)];
2303: if (i >= k) { /* only take upper triangular entry */
2304: cols[ncols_upper] = i;
2305: ncols_upper++;
2306: }
2307: }
2308: PetscCall(PetscLLAdd(ncols_upper, cols, am, &nlnk, lnk, lnkbt));
2309: nzk += nlnk;
2311: /* update lnk by computing fill-in for each pivot row to be merged in */
2312: prow = jl[k]; /* 1st pivot row */
2314: while (prow < k) {
2315: nextprow = jl[prow];
2316: /* merge prow into k-th row */
2317: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
2318: jmax = ui[prow + 1];
2319: ncols = jmax - jmin;
2320: uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:am-1) */
2321: PetscCall(PetscLLAddSorted(ncols, uj_ptr, am, &nlnk, lnk, lnkbt));
2322: nzk += nlnk;
2324: /* update il and jl for prow */
2325: if (jmin < jmax) {
2326: il[prow] = jmin;
2327: j = *uj_ptr;
2328: jl[prow] = jl[j];
2329: jl[j] = prow;
2330: }
2331: prow = nextprow;
2332: }
2334: /* if free space is not available, make more free space */
2335: if (current_space->local_remaining < nzk) {
2336: i = am - k + 1; /* num of unfactored rows */
2337: i = PetscIntMultTruncate(i, PetscMin(nzk, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
2338: PetscCall(PetscFreeSpaceGet(i, ¤t_space));
2339: reallocs++;
2340: }
2342: /* copy data into free space, then initialize lnk */
2343: PetscCall(PetscLLClean(am, am, nzk, lnk, current_space->array, lnkbt));
2345: /* add the k-th row into il and jl */
2346: if (nzk > 1) {
2347: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
2348: jl[k] = jl[i];
2349: jl[i] = k;
2350: il[k] = ui[k] + 1;
2351: }
2352: ui_ptr[k] = current_space->array;
2354: current_space->array += nzk;
2355: current_space->local_used += nzk;
2356: current_space->local_remaining -= nzk;
2358: ui[k + 1] = ui[k] + nzk;
2359: }
2361: PetscCall(ISRestoreIndices(perm, &rip));
2362: PetscCall(ISRestoreIndices(iperm, &riip));
2363: PetscCall(PetscFree4(ui_ptr, jl, il, cols));
2365: /* copy free_space into uj and free free_space; set ui, uj, udiag in new datastructure; */
2366: PetscCall(PetscShmgetAllocateArray(ui[am], sizeof(PetscInt), (void **)&uj));
2367: PetscCall(PetscFreeSpaceContiguous_Cholesky(&free_space, uj, am, ui, udiag)); /* store matrix factor */
2368: PetscCall(PetscLLDestroy(lnk, lnkbt));
2370: /* put together the new matrix in MATSEQSBAIJ format */
2371: b = (Mat_SeqSBAIJ *)fact->data;
2372: b->free_ij = PETSC_TRUE;
2373: PetscCall(PetscShmgetAllocateArray(ui[am], sizeof(PetscScalar), (void **)&b->a));
2374: b->free_a = PETSC_TRUE;
2375: b->j = uj;
2376: b->i = ui;
2377: b->diag = udiag;
2378: b->ilen = NULL;
2379: b->imax = NULL;
2380: b->row = perm;
2381: b->col = perm;
2383: PetscCall(PetscObjectReference((PetscObject)perm));
2384: PetscCall(PetscObjectReference((PetscObject)perm));
2386: b->icol = iperm;
2387: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
2389: PetscCall(PetscMalloc1(am, &b->solve_work));
2391: b->maxnz = b->nz = ui[am];
2393: fact->info.factor_mallocs = reallocs;
2394: fact->info.fill_ratio_given = fill;
2395: if (ai[am] != 0) {
2396: /* nonzeros in lower triangular part of A (including diagonals) = (ai[am]+am)/2 */
2397: fact->info.fill_ratio_needed = ((PetscReal)2 * ui[am]) / (ai[am] + am);
2398: } else {
2399: fact->info.fill_ratio_needed = 0.0;
2400: }
2401: #if PetscDefined(USE_INFO)
2402: if (ai[am] != 0) {
2403: PetscReal af = fact->info.fill_ratio_needed;
2404: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
2405: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
2406: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
2407: } else PetscCall(PetscInfo(A, "Empty matrix\n"));
2408: #endif
2409: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ;
2410: PetscFunctionReturn(PETSC_SUCCESS);
2411: }
2413: PetscErrorCode MatSolve_SeqAIJ_NaturalOrdering(Mat A, Vec bb, Vec xx)
2414: {
2415: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
2416: PetscInt n = A->rmap->n;
2417: const PetscInt *ai = a->i, *aj = a->j, *adiag = a->diag, *vi;
2418: PetscScalar *x, sum;
2419: const PetscScalar *b;
2420: const MatScalar *aa, *v;
2421: PetscInt i, nz;
2423: PetscFunctionBegin;
2424: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
2426: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
2427: PetscCall(VecGetArrayRead(bb, &b));
2428: PetscCall(VecGetArrayWrite(xx, &x));
2430: /* forward solve the lower triangular */
2431: x[0] = b[0];
2432: v = aa;
2433: vi = aj;
2434: for (i = 1; i < n; i++) {
2435: nz = ai[i + 1] - ai[i];
2436: sum = b[i];
2437: PetscSparseDenseMinusDot(sum, x, v, vi, nz);
2438: v += nz;
2439: vi += nz;
2440: x[i] = sum;
2441: }
2443: /* backward solve the upper triangular */
2444: for (i = n - 1; i >= 0; i--) {
2445: v = aa + adiag[i + 1] + 1;
2446: vi = aj + adiag[i + 1] + 1;
2447: nz = adiag[i] - adiag[i + 1] - 1;
2448: sum = x[i];
2449: PetscSparseDenseMinusDot(sum, x, v, vi, nz);
2450: x[i] = sum * v[nz]; /* x[i]=aa[adiag[i]]*sum; v++; */
2451: }
2453: PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
2454: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
2455: PetscCall(VecRestoreArrayRead(bb, &b));
2456: PetscCall(VecRestoreArrayWrite(xx, &x));
2457: PetscFunctionReturn(PETSC_SUCCESS);
2458: }
2460: PetscErrorCode MatSolve_SeqAIJ(Mat A, Vec bb, Vec xx)
2461: {
2462: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
2463: IS iscol = a->col, isrow = a->row;
2464: PetscInt i, n = A->rmap->n, *vi, *ai = a->i, *aj = a->j, *adiag = a->diag, nz;
2465: const PetscInt *rout, *cout, *r, *c;
2466: PetscScalar *x, *tmp, sum;
2467: const PetscScalar *b;
2468: const MatScalar *aa, *v;
2470: PetscFunctionBegin;
2471: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
2473: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
2474: PetscCall(VecGetArrayRead(bb, &b));
2475: PetscCall(VecGetArrayWrite(xx, &x));
2476: tmp = a->solve_work;
2478: PetscCall(ISGetIndices(isrow, &rout));
2479: r = rout;
2480: PetscCall(ISGetIndices(iscol, &cout));
2481: c = cout;
2483: /* forward solve the lower triangular */
2484: tmp[0] = b[r[0]];
2485: v = aa;
2486: vi = aj;
2487: for (i = 1; i < n; i++) {
2488: nz = ai[i + 1] - ai[i];
2489: sum = b[r[i]];
2490: PetscSparseDenseMinusDot(sum, tmp, v, vi, nz);
2491: tmp[i] = sum;
2492: v += nz;
2493: vi += nz;
2494: }
2496: /* backward solve the upper triangular */
2497: for (i = n - 1; i >= 0; i--) {
2498: v = aa + adiag[i + 1] + 1;
2499: vi = aj + adiag[i + 1] + 1;
2500: nz = adiag[i] - adiag[i + 1] - 1;
2501: sum = tmp[i];
2502: PetscSparseDenseMinusDot(sum, tmp, v, vi, nz);
2503: x[c[i]] = tmp[i] = sum * v[nz]; /* v[nz] = aa[adiag[i]] */
2504: }
2506: PetscCall(ISRestoreIndices(isrow, &rout));
2507: PetscCall(ISRestoreIndices(iscol, &cout));
2508: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
2509: PetscCall(VecRestoreArrayRead(bb, &b));
2510: PetscCall(VecRestoreArrayWrite(xx, &x));
2511: PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
2512: PetscFunctionReturn(PETSC_SUCCESS);
2513: }