Actual source code: sbaijfact.c
1: #include <../src/mat/impls/baij/seq/baij.h>
2: #include <../src/mat/impls/sbaij/seq/sbaij.h>
3: #include <petsc/private/kernels/blockinvert.h>
4: #include <petscis.h>
6: PetscErrorCode MatGetInertia_SeqSBAIJ(Mat F, PetscInt *nneg, PetscInt *nzero, PetscInt *npos)
7: {
8: Mat_SeqSBAIJ *fact = (Mat_SeqSBAIJ *)F->data;
9: MatScalar *dd = fact->a;
10: PetscInt mbs = fact->mbs, bs = F->rmap->bs, i, nneg_tmp, npos_tmp, *fi = fact->diag;
12: PetscFunctionBegin;
13: PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for bs: %" PetscInt_FMT " >1 yet", bs);
15: nneg_tmp = 0;
16: npos_tmp = 0;
17: if (fi) {
18: for (i = 0; i < mbs; i++) {
19: if (PetscRealPart(dd[*fi]) > 0.0) npos_tmp++;
20: else if (PetscRealPart(dd[*fi]) < 0.0) nneg_tmp++;
21: fi++;
22: }
23: } else {
24: for (i = 0; i < mbs; i++) {
25: if (PetscRealPart(dd[fact->i[i]]) > 0.0) npos_tmp++;
26: else if (PetscRealPart(dd[fact->i[i]]) < 0.0) nneg_tmp++;
27: }
28: }
29: if (nneg) *nneg = nneg_tmp;
30: if (npos) *npos = npos_tmp;
31: if (nzero) *nzero = mbs - nneg_tmp - npos_tmp;
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: /*
36: Symbolic U^T*D*U factorization for SBAIJ format. Modified from SSF of YSMP.
37: Use Modified Sparse Row (MSR) storage for u and ju. See page 85, "Iterative Methods ..." by Saad.
38: */
39: static PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ_MSR(Mat F, Mat A, IS perm, const MatFactorInfo *info)
40: {
41: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b;
42: const PetscInt *rip, *ai, *aj;
43: PetscInt i, mbs = a->mbs, *jutmp, bs = A->rmap->bs, bs2 = a->bs2;
44: PetscInt m, reallocs = 0, prow;
45: PetscInt *jl, *q, jmin, jmax, juidx, nzk, qm, *iu, *ju, k, j, vj, umax, maxadd;
46: PetscReal f = info->fill;
47: PetscBool perm_identity;
49: PetscFunctionBegin;
50: /* check whether perm is the identity mapping */
51: PetscCall(ISIdentity(perm, &perm_identity));
52: PetscCall(ISGetIndices(perm, &rip));
54: if (perm_identity) { /* without permutation */
55: a->permute = PETSC_FALSE;
57: ai = a->i;
58: aj = a->j;
59: } else { /* non-trivial permutation */
60: a->permute = PETSC_TRUE;
62: PetscCall(MatReorderingSeqSBAIJ(A, perm));
64: ai = a->inew;
65: aj = a->jnew;
66: }
68: /* initialization */
69: PetscCall(PetscShmgetAllocateArray(mbs + 1, sizeof(PetscInt), (void **)&iu));
70: umax = (PetscInt)(f * ai[mbs] + 1);
71: umax += mbs + 1;
72: PetscCall(PetscShmgetAllocateArray(umax, sizeof(PetscInt), (void **)&ju));
73: iu[0] = mbs + 1;
74: juidx = mbs + 1; /* index for ju */
75: /* jl linked list for pivot row -- linked list for col index */
76: PetscCall(PetscMalloc2(mbs, &jl, mbs, &q));
77: for (i = 0; i < mbs; i++) {
78: jl[i] = mbs;
79: q[i] = 0;
80: }
82: /* for each row k */
83: for (k = 0; k < mbs; k++) {
84: for (i = 0; i < mbs; i++) q[i] = 0; /* to be removed! */
85: nzk = 0; /* num. of nz blocks in k-th block row with diagonal block excluded */
86: q[k] = mbs;
87: /* initialize nonzero structure of k-th row to row rip[k] of A */
88: jmin = ai[rip[k]] + 1; /* exclude diag[k] */
89: jmax = ai[rip[k] + 1];
90: for (j = jmin; j < jmax; j++) {
91: vj = rip[aj[j]]; /* col. value */
92: if (vj > k) {
93: qm = k;
94: do {
95: m = qm;
96: qm = q[m];
97: } while (qm < vj);
98: PetscCheck(qm != vj, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Duplicate entry in A");
99: nzk++;
100: q[m] = vj;
101: q[vj] = qm;
102: } /* if (vj > k) */
103: } /* for (j=jmin; j<jmax; j++) */
105: /* modify nonzero structure of k-th row by computing fill-in
106: for each row i to be merged in */
107: prow = k;
108: prow = jl[prow]; /* next pivot row (== mbs for symbolic factorization) */
110: while (prow < k) {
111: /* merge row prow into k-th row */
112: jmin = iu[prow] + 1;
113: jmax = iu[prow + 1];
114: qm = k;
115: for (j = jmin; j < jmax; j++) {
116: vj = ju[j];
117: do {
118: m = qm;
119: qm = q[m];
120: } while (qm < vj);
121: if (qm != vj) {
122: nzk++;
123: q[m] = vj;
124: q[vj] = qm;
125: qm = vj;
126: }
127: }
128: prow = jl[prow]; /* next pivot row */
129: }
131: /* add k to row list for first nonzero element in k-th row */
132: if (nzk > 0) {
133: i = q[k]; /* col value of first nonzero element in U(k, k+1:mbs-1) */
134: jl[k] = jl[i];
135: jl[i] = k;
136: }
137: iu[k + 1] = iu[k] + nzk;
139: /* allocate more space to ju if needed */
140: if (iu[k + 1] > umax) {
141: /* estimate how much additional space we will need */
142: /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
143: /* just double the memory each time */
144: maxadd = umax;
145: if (maxadd < nzk) maxadd = (mbs - k) * (nzk + 1) / 2;
146: umax += maxadd;
148: /* allocate a longer ju */
149: PetscCall(PetscShmgetAllocateArray(umax, sizeof(PetscInt), (void **)&jutmp));
150: PetscCall(PetscArraycpy(jutmp, ju, iu[k]));
151: PetscCall(PetscShmgetDeallocateArray((void **)&ju));
152: ju = jutmp;
153: reallocs++; /* count how many times we realloc */
154: }
156: /* save nonzero structure of k-th row in ju */
157: i = k;
158: while (nzk--) {
159: i = q[i];
160: ju[juidx++] = i;
161: }
162: }
164: #if PetscDefined(USE_INFO)
165: if (ai[mbs] != 0) {
166: PetscReal af = ((PetscReal)iu[mbs]) / ((PetscReal)ai[mbs]);
167: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)f, (double)af));
168: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
169: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g);\n", (double)af));
170: PetscCall(PetscInfo(A, "for best performance.\n"));
171: } else PetscCall(PetscInfo(A, "Empty matrix\n"));
172: #endif
174: PetscCall(ISRestoreIndices(perm, &rip));
175: PetscCall(PetscFree2(jl, q));
177: /* put together the new matrix */
178: PetscCall(MatSeqSBAIJSetPreallocation(F, bs, MAT_SKIP_ALLOCATION, NULL));
180: b = (Mat_SeqSBAIJ *)F->data;
181: b->free_ij = PETSC_TRUE;
182: PetscCall(PetscShmgetAllocateArray((iu[mbs] + 1) * bs2, sizeof(PetscScalar), (void **)&b->a));
183: b->free_a = PETSC_TRUE;
184: b->j = ju;
185: b->i = iu;
186: b->diag = NULL;
187: b->ilen = NULL;
188: b->imax = NULL;
189: b->row = perm;
191: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
193: PetscCall(PetscObjectReference((PetscObject)perm));
195: b->icol = perm;
196: PetscCall(PetscObjectReference((PetscObject)perm));
197: PetscCall(PetscMalloc1(bs * mbs + bs, &b->solve_work));
198: /* In b structure: Free imax, ilen, old a, old j.
199: Allocate idnew, solve_work, new a, new j */
200: b->maxnz = b->nz = iu[mbs];
202: F->info.factor_mallocs = reallocs;
203: F->info.fill_ratio_given = f;
204: if (ai[mbs] != 0) {
205: F->info.fill_ratio_needed = ((PetscReal)iu[mbs]) / ((PetscReal)ai[mbs]);
206: } else {
207: F->info.fill_ratio_needed = 0.0;
208: }
209: PetscCall(MatSeqSBAIJSetNumericFactorization_inplace(F, perm_identity));
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
212: /*
213: Symbolic U^T*D*U factorization for SBAIJ format.
214: See MatICCFactorSymbolic_SeqAIJ() for description of its data structure.
215: */
216: #include <../src/mat/utils/freespace.h>
217: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
218: {
219: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
220: Mat_SeqSBAIJ *b;
221: PetscBool perm_identity;
222: PetscReal fill = info->fill;
223: const PetscInt *rip, *ai = a->i, *aj = a->j;
224: PetscInt i, mbs = a->mbs, bs = A->rmap->bs, reallocs = 0, prow;
225: PetscInt *jl, jmin, jmax, nzk, *ui, k, j, *il, nextprow;
226: PetscInt nlnk, *lnk, ncols, *cols, *uj, **ui_ptr, *uj_ptr, *udiag;
227: PetscFreeSpaceList free_space = NULL, current_space = NULL;
228: PetscBT lnkbt;
229: PetscBool diagDense;
231: PetscFunctionBegin;
232: 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);
233: PetscCall(MatGetDiagonalMarkers_SeqSBAIJ(A, NULL, &diagDense));
234: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry");
235: if (bs > 1) {
236: PetscCall(MatCholeskyFactorSymbolic_SeqSBAIJ_inplace(fact, A, perm, info));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /* check whether perm is the identity mapping */
241: PetscCall(ISIdentity(perm, &perm_identity));
242: PetscCheck(perm_identity, PETSC_COMM_SELF, PETSC_ERR_SUP, "Matrix reordering is not supported for sbaij matrix. Use aij format");
243: a->permute = PETSC_FALSE;
244: PetscCall(ISGetIndices(perm, &rip));
246: /* initialization */
247: PetscCall(PetscShmgetAllocateArray(mbs + 1, sizeof(PetscInt), (void **)&ui));
248: PetscCall(PetscMalloc1(mbs + 1, &udiag));
249: ui[0] = 0;
251: /* jl: linked list for storing indices of the pivot rows
252: il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
253: PetscCall(PetscMalloc4(mbs, &ui_ptr, mbs, &il, mbs, &jl, mbs, &cols));
254: for (i = 0; i < mbs; i++) {
255: jl[i] = mbs;
256: il[i] = 0;
257: }
259: /* create and initialize a linked list for storing column indices of the active row k */
260: nlnk = mbs + 1;
261: PetscCall(PetscLLCreate(mbs, mbs, nlnk, lnk, lnkbt));
263: /* initial FreeSpace size is fill*(ai[mbs]+1) */
264: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, ai[mbs] + 1), &free_space));
265: current_space = free_space;
267: for (k = 0; k < mbs; k++) { /* for each active row k */
268: /* initialize lnk by the column indices of row rip[k] of A */
269: nzk = 0;
270: ncols = ai[k + 1] - ai[k];
271: PetscCheck(ncols, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Empty row %" PetscInt_FMT " in matrix ", k);
272: for (j = 0; j < ncols; j++) {
273: i = *(aj + ai[k] + j);
274: cols[j] = i;
275: }
276: PetscCall(PetscLLAdd(ncols, cols, mbs, &nlnk, lnk, lnkbt));
277: nzk += nlnk;
279: /* update lnk by computing fill-in for each pivot row to be merged in */
280: prow = jl[k]; /* 1st pivot row */
282: while (prow < k) {
283: nextprow = jl[prow];
284: /* merge prow into k-th row */
285: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
286: jmax = ui[prow + 1];
287: ncols = jmax - jmin;
288: uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
289: PetscCall(PetscLLAddSorted(ncols, uj_ptr, mbs, &nlnk, lnk, lnkbt));
290: nzk += nlnk;
292: /* update il and jl for prow */
293: if (jmin < jmax) {
294: il[prow] = jmin;
295: j = *uj_ptr;
296: jl[prow] = jl[j];
297: jl[j] = prow;
298: }
299: prow = nextprow;
300: }
302: /* if free space is not available, make more free space */
303: if (current_space->local_remaining < nzk) {
304: i = mbs - k + 1; /* num of unfactored rows */
305: i = PetscIntMultTruncate(i, PetscMin(nzk, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
306: PetscCall(PetscFreeSpaceGet(i, ¤t_space));
307: reallocs++;
308: }
310: /* copy data into free space, then initialize lnk */
311: PetscCall(PetscLLClean(mbs, mbs, nzk, lnk, current_space->array, lnkbt));
313: /* add the k-th row into il and jl */
314: if (nzk > 1) {
315: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
316: jl[k] = jl[i];
317: jl[i] = k;
318: il[k] = ui[k] + 1;
319: }
320: ui_ptr[k] = current_space->array;
322: current_space->array += nzk;
323: current_space->local_used += nzk;
324: current_space->local_remaining -= nzk;
326: ui[k + 1] = ui[k] + nzk;
327: }
329: PetscCall(ISRestoreIndices(perm, &rip));
330: PetscCall(PetscFree4(ui_ptr, il, jl, cols));
332: /* destroy list of free space and other temporary array(s) */
333: PetscCall(PetscShmgetAllocateArray(ui[mbs], sizeof(PetscInt), (void **)&uj));
334: PetscCall(PetscFreeSpaceContiguous_Cholesky(&free_space, uj, mbs, ui, udiag)); /* store matrix factor */
335: PetscCall(PetscLLDestroy(lnk, lnkbt));
337: /* put together the new matrix in MATSEQSBAIJ format */
338: PetscCall(MatSeqSBAIJSetPreallocation(fact, bs, MAT_SKIP_ALLOCATION, NULL));
340: b = (Mat_SeqSBAIJ *)fact->data;
341: b->free_ij = PETSC_TRUE;
342: PetscCall(PetscShmgetAllocateArray(ui[mbs], sizeof(PetscScalar), (void **)&b->a));
343: b->free_a = PETSC_TRUE;
344: b->j = uj;
345: b->i = ui;
346: b->diag = udiag;
347: b->ilen = NULL;
348: b->imax = NULL;
349: b->row = perm;
350: b->icol = perm;
352: PetscCall(PetscObjectReference((PetscObject)perm));
353: PetscCall(PetscObjectReference((PetscObject)perm));
355: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
357: PetscCall(PetscMalloc1(mbs + 1, &b->solve_work));
359: b->maxnz = b->nz = ui[mbs];
361: fact->info.factor_mallocs = reallocs;
362: fact->info.fill_ratio_given = fill;
363: if (ai[mbs] != 0) {
364: fact->info.fill_ratio_needed = ((PetscReal)ui[mbs]) / ai[mbs];
365: } else {
366: fact->info.fill_ratio_needed = 0.0;
367: }
368: #if PetscDefined(USE_INFO)
369: if (ai[mbs] != 0) {
370: PetscReal af = fact->info.fill_ratio_needed;
371: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
372: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
373: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
374: } else PetscCall(PetscInfo(A, "Empty matrix\n"));
375: #endif
376: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering;
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ_inplace(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
381: {
382: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
383: Mat_SeqSBAIJ *b;
384: PetscBool perm_identity;
385: PetscReal fill = info->fill;
386: const PetscInt *rip, *ai, *aj;
387: PetscInt i, mbs = a->mbs, bs = A->rmap->bs, reallocs = 0, prow;
388: PetscInt *jl, jmin, jmax, nzk, *ui, k, j, *il, nextprow;
389: PetscInt nlnk, *lnk, ncols, *cols, *uj, **ui_ptr, *uj_ptr;
390: PetscFreeSpaceList free_space = NULL, current_space = NULL;
391: PetscBT lnkbt;
392: PetscBool diagDense;
394: PetscFunctionBegin;
395: PetscCall(MatGetDiagonalMarkers_SeqSBAIJ(A, NULL, &diagDense));
396: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry");
398: /*
399: This code originally uses Modified Sparse Row (MSR) storage
400: (see page 85, "Iterative Methods ..." by Saad) for the output matrix B - bad choice!
401: Then it is rewritten so the factor B takes seqsbaij format. However the associated
402: MatCholeskyFactorNumeric_() have not been modified for the cases of bs>1 or !perm_identity,
403: thus the original code in MSR format is still used for these cases.
404: The code below should replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR() whenever
405: MatCholeskyFactorNumeric_() is modified for using sbaij symbolic factor.
406: */
407: if (bs > 1) {
408: PetscCall(MatCholeskyFactorSymbolic_SeqSBAIJ_MSR(fact, A, perm, info));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: /* check whether perm is the identity mapping */
413: PetscCall(ISIdentity(perm, &perm_identity));
414: PetscCheck(perm_identity, PETSC_COMM_SELF, PETSC_ERR_SUP, "Matrix reordering is not supported for sbaij matrix. Use aij format");
415: a->permute = PETSC_FALSE;
416: ai = a->i;
417: aj = a->j;
418: PetscCall(ISGetIndices(perm, &rip));
420: /* initialization */
421: PetscCall(PetscShmgetAllocateArray(mbs + 1, sizeof(PetscInt), (void **)&ui));
422: ui[0] = 0;
424: /* jl: linked list for storing indices of the pivot rows
425: il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
426: PetscCall(PetscMalloc4(mbs, &ui_ptr, mbs, &il, mbs, &jl, mbs, &cols));
427: for (i = 0; i < mbs; i++) {
428: jl[i] = mbs;
429: il[i] = 0;
430: }
432: /* create and initialize a linked list for storing column indices of the active row k */
433: nlnk = mbs + 1;
434: PetscCall(PetscLLCreate(mbs, mbs, nlnk, lnk, lnkbt));
436: /* initial FreeSpace size is fill*(ai[mbs]+1) */
437: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, ai[mbs] + 1), &free_space));
438: current_space = free_space;
440: for (k = 0; k < mbs; k++) { /* for each active row k */
441: /* initialize lnk by the column indices of row rip[k] of A */
442: nzk = 0;
443: ncols = ai[rip[k] + 1] - ai[rip[k]];
444: for (j = 0; j < ncols; j++) {
445: i = *(aj + ai[rip[k]] + j);
446: cols[j] = rip[i];
447: }
448: PetscCall(PetscLLAdd(ncols, cols, mbs, &nlnk, lnk, lnkbt));
449: nzk += nlnk;
451: /* update lnk by computing fill-in for each pivot row to be merged in */
452: prow = jl[k]; /* 1st pivot row */
454: while (prow < k) {
455: nextprow = jl[prow];
456: /* merge prow into k-th row */
457: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
458: jmax = ui[prow + 1];
459: ncols = jmax - jmin;
460: uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
461: PetscCall(PetscLLAddSorted(ncols, uj_ptr, mbs, &nlnk, lnk, lnkbt));
462: nzk += nlnk;
464: /* update il and jl for prow */
465: if (jmin < jmax) {
466: il[prow] = jmin;
468: j = *uj_ptr;
469: jl[prow] = jl[j];
470: jl[j] = prow;
471: }
472: prow = nextprow;
473: }
475: /* if free space is not available, make more free space */
476: if (current_space->local_remaining < nzk) {
477: i = mbs - k + 1; /* num of unfactored rows */
478: i = PetscMin(PetscIntMultTruncate(i, nzk), PetscIntMultTruncate(i, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
479: PetscCall(PetscFreeSpaceGet(i, ¤t_space));
480: reallocs++;
481: }
483: /* copy data into free space, then initialize lnk */
484: PetscCall(PetscLLClean(mbs, mbs, nzk, lnk, current_space->array, lnkbt));
486: /* add the k-th row into il and jl */
487: if (nzk - 1 > 0) {
488: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
489: jl[k] = jl[i];
490: jl[i] = k;
491: il[k] = ui[k] + 1;
492: }
493: ui_ptr[k] = current_space->array;
495: current_space->array += nzk;
496: current_space->local_used += nzk;
497: current_space->local_remaining -= nzk;
499: ui[k + 1] = ui[k] + nzk;
500: }
502: PetscCall(ISRestoreIndices(perm, &rip));
503: PetscCall(PetscFree4(ui_ptr, il, jl, cols));
505: /* destroy list of free space and other temporary array(s) */
506: PetscCall(PetscShmgetAllocateArray(ui[mbs] + 1, sizeof(PetscInt), (void **)&uj));
507: PetscCall(PetscFreeSpaceContiguous(&free_space, uj));
508: PetscCall(PetscLLDestroy(lnk, lnkbt));
510: /* put together the new matrix in MATSEQSBAIJ format */
511: PetscCall(MatSeqSBAIJSetPreallocation(fact, bs, MAT_SKIP_ALLOCATION, NULL));
513: b = (Mat_SeqSBAIJ *)fact->data;
514: PetscCall(PetscShmgetAllocateArray(ui[mbs] + 1, sizeof(PetscScalar), (void **)&b->a));
515: b->free_a = PETSC_TRUE;
516: b->free_ij = PETSC_TRUE;
517: b->j = uj;
518: b->i = ui;
519: b->diag = NULL;
520: b->ilen = NULL;
521: b->imax = NULL;
522: b->row = perm;
524: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
526: PetscCall(PetscObjectReference((PetscObject)perm));
527: b->icol = perm;
528: PetscCall(PetscObjectReference((PetscObject)perm));
529: PetscCall(PetscMalloc1(mbs + 1, &b->solve_work));
530: b->maxnz = b->nz = ui[mbs];
532: fact->info.factor_mallocs = reallocs;
533: fact->info.fill_ratio_given = fill;
534: if (ai[mbs] != 0) {
535: fact->info.fill_ratio_needed = ((PetscReal)ui[mbs]) / ai[mbs];
536: } else {
537: fact->info.fill_ratio_needed = 0.0;
538: }
539: #if PetscDefined(USE_INFO)
540: if (ai[mbs] != 0) {
541: PetscReal af = fact->info.fill_ratio_needed;
542: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
543: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
544: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
545: } else PetscCall(PetscInfo(A, "Empty matrix\n"));
546: #endif
547: PetscCall(MatSeqSBAIJSetNumericFactorization_inplace(fact, perm_identity));
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat C, Mat A, const MatFactorInfo *info)
552: {
553: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
554: IS perm = b->row;
555: const PetscInt *ai, *aj, *perm_ptr, mbs = a->mbs, *bi = b->i, *bj = b->j;
556: PetscInt i, j;
557: PetscInt *a2anew, k, k1, jmin, jmax, *jl, *il, vj, nexti, ili;
558: PetscInt bs = A->rmap->bs, bs2 = a->bs2;
559: MatScalar *ba = b->a, *aa, *ap, *dk, *uik;
560: MatScalar *u, *diag, *rtmp, *rtmp_ptr;
561: MatScalar *work;
562: PetscInt *pivots;
563: PetscBool allowzeropivot, zeropivotdetected;
565: PetscFunctionBegin;
566: /* initialization */
567: PetscCall(PetscCalloc1(bs2 * mbs, &rtmp));
568: PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
569: allowzeropivot = PetscNot(A->erroriffailure);
571: il[0] = 0;
572: for (i = 0; i < mbs; i++) jl[i] = mbs;
574: PetscCall(PetscMalloc3(bs2, &dk, bs2, &uik, bs, &work));
575: PetscCall(PetscMalloc1(bs, &pivots));
577: PetscCall(ISGetIndices(perm, &perm_ptr));
579: /* check permutation */
580: if (!a->permute) {
581: ai = a->i;
582: aj = a->j;
583: aa = a->a;
584: } else {
585: ai = a->inew;
586: aj = a->jnew;
587: PetscCall(PetscMalloc1(bs2 * ai[mbs], &aa));
588: PetscCall(PetscArraycpy(aa, a->a, bs2 * ai[mbs]));
589: PetscCall(PetscMalloc1(ai[mbs], &a2anew));
590: PetscCall(PetscArraycpy(a2anew, a->a2anew, ai[mbs]));
592: for (i = 0; i < mbs; i++) {
593: jmin = ai[i];
594: jmax = ai[i + 1];
595: for (j = jmin; j < jmax; j++) {
596: while (a2anew[j] != j) {
597: k = a2anew[j];
598: a2anew[j] = a2anew[k];
599: a2anew[k] = k;
600: for (k1 = 0; k1 < bs2; k1++) {
601: dk[k1] = aa[k * bs2 + k1];
602: aa[k * bs2 + k1] = aa[j * bs2 + k1];
603: aa[j * bs2 + k1] = dk[k1];
604: }
605: }
606: /* transform column-oriented blocks that lie in the lower triangle to row-oriented blocks */
607: if (i > aj[j]) {
608: ap = aa + j * bs2; /* ptr to the beginning of j-th block of aa */
609: for (k = 0; k < bs2; k++) dk[k] = ap[k]; /* dk <- j-th block of aa */
610: for (k = 0; k < bs; k++) { /* j-th block of aa <- dk^T */
611: for (k1 = 0; k1 < bs; k1++) *ap++ = dk[k + bs * k1];
612: }
613: }
614: }
615: }
616: PetscCall(PetscFree(a2anew));
617: }
619: /* for each row k */
620: for (k = 0; k < mbs; k++) {
621: /*initialize k-th row with elements nonzero in row perm(k) of A */
622: jmin = ai[perm_ptr[k]];
623: jmax = ai[perm_ptr[k] + 1];
625: ap = aa + jmin * bs2;
626: for (j = jmin; j < jmax; j++) {
627: vj = perm_ptr[aj[j]]; /* block col. index */
628: rtmp_ptr = rtmp + vj * bs2;
629: for (i = 0; i < bs2; i++) *rtmp_ptr++ = *ap++;
630: }
632: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
633: PetscCall(PetscArraycpy(dk, rtmp + k * bs2, bs2));
634: i = jl[k]; /* first row to be added to k_th row */
636: while (i < k) {
637: nexti = jl[i]; /* next row to be added to k_th row */
639: /* compute multiplier */
640: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
642: /* uik = -inv(Di)*U_bar(i,k) */
643: diag = ba + i * bs2;
644: u = ba + ili * bs2;
645: PetscCall(PetscArrayzero(uik, bs2));
646: PetscKernel_A_gets_A_minus_B_times_C(bs, uik, diag, u);
648: /* update D(k) += -U(i,k)^T * U_bar(i,k) */
649: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs, dk, uik, u);
650: PetscCall(PetscLogFlops(4.0 * bs * bs2));
652: /* update -U(i,k) */
653: PetscCall(PetscArraycpy(ba + ili * bs2, uik, bs2));
655: /* add multiple of row i to k-th row ... */
656: jmin = ili + 1;
657: jmax = bi[i + 1];
658: if (jmin < jmax) {
659: for (j = jmin; j < jmax; j++) {
660: /* rtmp += -U(i,k)^T * U_bar(i,j) */
661: rtmp_ptr = rtmp + bj[j] * bs2;
662: u = ba + j * bs2;
663: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs, rtmp_ptr, uik, u);
664: }
665: PetscCall(PetscLogFlops(2.0 * bs * bs2 * (jmax - jmin)));
667: /* ... add i to row list for next nonzero entry */
668: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
669: j = bj[jmin];
670: jl[i] = jl[j];
671: jl[j] = i; /* update jl */
672: }
673: i = nexti;
674: }
676: /* save nonzero entries in k-th row of U ... */
678: /* invert diagonal block */
679: diag = ba + k * bs2;
680: PetscCall(PetscArraycpy(diag, dk, bs2));
682: PetscCall(PetscKernel_A_gets_inverse_A(bs, diag, pivots, work, allowzeropivot, &zeropivotdetected));
683: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
685: jmin = bi[k];
686: jmax = bi[k + 1];
687: if (jmin < jmax) {
688: for (j = jmin; j < jmax; j++) {
689: vj = bj[j]; /* block col. index of U */
690: u = ba + j * bs2;
691: rtmp_ptr = rtmp + vj * bs2;
692: for (k1 = 0; k1 < bs2; k1++) {
693: *u++ = *rtmp_ptr;
694: *rtmp_ptr++ = 0.0;
695: }
696: }
698: /* ... add k to row list for first nonzero entry in k-th row */
699: il[k] = jmin;
700: i = bj[jmin];
701: jl[k] = jl[i];
702: jl[i] = k;
703: }
704: }
706: PetscCall(PetscFree(rtmp));
707: PetscCall(PetscFree2(il, jl));
708: PetscCall(PetscFree3(dk, uik, work));
709: PetscCall(PetscFree(pivots));
710: if (a->permute) PetscCall(PetscFree(aa));
712: PetscCall(ISRestoreIndices(perm, &perm_ptr));
714: C->ops->solve = MatSolve_SeqSBAIJ_N_inplace;
715: C->ops->solvetranspose = MatSolve_SeqSBAIJ_N_inplace;
716: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_N_inplace;
717: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_N_inplace;
719: C->assembled = PETSC_TRUE;
720: C->preallocated = PETSC_TRUE;
722: PetscCall(PetscLogFlops(1.3333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }
726: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering(Mat C, Mat A, const MatFactorInfo *info)
727: {
728: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
729: PetscInt i, j, mbs = a->mbs, *bi = b->i, *bj = b->j;
730: PetscInt *ai, *aj, k, k1, jmin, jmax, *jl, *il, vj, nexti, ili;
731: PetscInt bs = A->rmap->bs, bs2 = a->bs2;
732: MatScalar *ba = b->a, *aa, *ap, *dk, *uik;
733: MatScalar *u, *diag, *rtmp, *rtmp_ptr;
734: MatScalar *work;
735: PetscInt *pivots;
736: PetscBool allowzeropivot, zeropivotdetected;
738: PetscFunctionBegin;
739: PetscCall(PetscCalloc1(bs2 * mbs, &rtmp));
740: PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
741: il[0] = 0;
742: for (i = 0; i < mbs; i++) jl[i] = mbs;
744: PetscCall(PetscMalloc3(bs2, &dk, bs2, &uik, bs, &work));
745: PetscCall(PetscMalloc1(bs, &pivots));
746: allowzeropivot = PetscNot(A->erroriffailure);
748: ai = a->i;
749: aj = a->j;
750: aa = a->a;
752: /* for each row k */
753: for (k = 0; k < mbs; k++) {
754: /*initialize k-th row with elements nonzero in row k of A */
755: jmin = ai[k];
756: jmax = ai[k + 1];
757: ap = aa + jmin * bs2;
758: for (j = jmin; j < jmax; j++) {
759: vj = aj[j]; /* block col. index */
760: rtmp_ptr = rtmp + vj * bs2;
761: for (i = 0; i < bs2; i++) *rtmp_ptr++ = *ap++;
762: }
764: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
765: PetscCall(PetscArraycpy(dk, rtmp + k * bs2, bs2));
766: i = jl[k]; /* first row to be added to k_th row */
768: while (i < k) {
769: nexti = jl[i]; /* next row to be added to k_th row */
771: /* compute multiplier */
772: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
774: /* uik = -inv(Di)*U_bar(i,k) */
775: diag = ba + i * bs2;
776: u = ba + ili * bs2;
777: PetscCall(PetscArrayzero(uik, bs2));
778: PetscKernel_A_gets_A_minus_B_times_C(bs, uik, diag, u);
780: /* update D(k) += -U(i,k)^T * U_bar(i,k) */
781: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs, dk, uik, u);
782: PetscCall(PetscLogFlops(2.0 * bs * bs2));
784: /* update -U(i,k) */
785: PetscCall(PetscArraycpy(ba + ili * bs2, uik, bs2));
787: /* add multiple of row i to k-th row ... */
788: jmin = ili + 1;
789: jmax = bi[i + 1];
790: if (jmin < jmax) {
791: for (j = jmin; j < jmax; j++) {
792: /* rtmp += -U(i,k)^T * U_bar(i,j) */
793: rtmp_ptr = rtmp + bj[j] * bs2;
794: u = ba + j * bs2;
795: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs, rtmp_ptr, uik, u);
796: }
797: PetscCall(PetscLogFlops(2.0 * bs * bs2 * (jmax - jmin)));
799: /* ... add i to row list for next nonzero entry */
800: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
801: j = bj[jmin];
802: jl[i] = jl[j];
803: jl[j] = i; /* update jl */
804: }
805: i = nexti;
806: }
808: /* save nonzero entries in k-th row of U ... */
810: /* invert diagonal block */
811: diag = ba + k * bs2;
812: PetscCall(PetscArraycpy(diag, dk, bs2));
814: PetscCall(PetscKernel_A_gets_inverse_A(bs, diag, pivots, work, allowzeropivot, &zeropivotdetected));
815: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
817: jmin = bi[k];
818: jmax = bi[k + 1];
819: if (jmin < jmax) {
820: for (j = jmin; j < jmax; j++) {
821: vj = bj[j]; /* block col. index of U */
822: u = ba + j * bs2;
823: rtmp_ptr = rtmp + vj * bs2;
824: for (k1 = 0; k1 < bs2; k1++) {
825: *u++ = *rtmp_ptr;
826: *rtmp_ptr++ = 0.0;
827: }
828: }
830: /* ... add k to row list for first nonzero entry in k-th row */
831: il[k] = jmin;
832: i = bj[jmin];
833: jl[k] = jl[i];
834: jl[i] = k;
835: }
836: }
838: PetscCall(PetscFree(rtmp));
839: PetscCall(PetscFree2(il, jl));
840: PetscCall(PetscFree3(dk, uik, work));
841: PetscCall(PetscFree(pivots));
843: C->ops->solve = MatSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
844: C->ops->solvetranspose = MatSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
845: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
846: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
847: C->assembled = PETSC_TRUE;
848: C->preallocated = PETSC_TRUE;
850: PetscCall(PetscLogFlops(1.3333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */
851: PetscFunctionReturn(PETSC_SUCCESS);
852: }
854: /*
855: Numeric U^T*D*U factorization for SBAIJ format. Modified from SNF of YSMP.
856: Version for blocks 2 by 2.
857: */
858: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat C, Mat A, const MatFactorInfo *info)
859: {
860: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
861: IS perm = b->row;
862: const PetscInt *ai, *aj, *perm_ptr;
863: PetscInt i, j, mbs = a->mbs, *bi = b->i, *bj = b->j;
864: PetscInt *a2anew, k, k1, jmin, jmax, *jl, *il, vj, nexti, ili;
865: MatScalar *ba = b->a, *aa, *ap;
866: MatScalar *u, *diag, *rtmp, *rtmp_ptr, dk[4], uik[4];
867: PetscReal shift = info->shiftamount;
868: PetscBool allowzeropivot, zeropivotdetected;
870: PetscFunctionBegin;
871: allowzeropivot = PetscNot(A->erroriffailure);
873: /* initialization */
874: /* il and jl record the first nonzero element in each row of the accessing
875: window U(0:k, k:mbs-1).
876: jl: list of rows to be added to uneliminated rows
877: i>= k: jl(i) is the first row to be added to row i
878: i< k: jl(i) is the row following row i in some list of rows
879: jl(i) = mbs indicates the end of a list
880: il(i): points to the first nonzero element in columns k,...,mbs-1 of
881: row i of U */
882: PetscCall(PetscCalloc1(4 * mbs, &rtmp));
883: PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
884: il[0] = 0;
885: for (i = 0; i < mbs; i++) jl[i] = mbs;
887: PetscCall(ISGetIndices(perm, &perm_ptr));
889: /* check permutation */
890: if (!a->permute) {
891: ai = a->i;
892: aj = a->j;
893: aa = a->a;
894: } else {
895: ai = a->inew;
896: aj = a->jnew;
897: PetscCall(PetscMalloc1(4 * ai[mbs], &aa));
898: PetscCall(PetscArraycpy(aa, a->a, 4 * ai[mbs]));
899: PetscCall(PetscMalloc1(ai[mbs], &a2anew));
900: PetscCall(PetscArraycpy(a2anew, a->a2anew, ai[mbs]));
902: for (i = 0; i < mbs; i++) {
903: jmin = ai[i];
904: jmax = ai[i + 1];
905: for (j = jmin; j < jmax; j++) {
906: while (a2anew[j] != j) {
907: k = a2anew[j];
908: a2anew[j] = a2anew[k];
909: a2anew[k] = k;
910: for (k1 = 0; k1 < 4; k1++) {
911: dk[k1] = aa[k * 4 + k1];
912: aa[k * 4 + k1] = aa[j * 4 + k1];
913: aa[j * 4 + k1] = dk[k1];
914: }
915: }
916: /* transform column-oriented blocks that lie in the lower triangle to row-oriented blocks */
917: if (i > aj[j]) {
918: ap = aa + j * 4; /* ptr to the beginning of the block */
919: dk[1] = ap[1]; /* swap ap[1] and ap[2] */
920: ap[1] = ap[2];
921: ap[2] = dk[1];
922: }
923: }
924: }
925: PetscCall(PetscFree(a2anew));
926: }
928: /* for each row k */
929: for (k = 0; k < mbs; k++) {
930: /*initialize k-th row with elements nonzero in row perm(k) of A */
931: jmin = ai[perm_ptr[k]];
932: jmax = ai[perm_ptr[k] + 1];
933: ap = aa + jmin * 4;
934: for (j = jmin; j < jmax; j++) {
935: vj = perm_ptr[aj[j]]; /* block col. index */
936: rtmp_ptr = rtmp + vj * 4;
937: for (i = 0; i < 4; i++) *rtmp_ptr++ = *ap++;
938: }
940: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
941: PetscCall(PetscArraycpy(dk, rtmp + k * 4, 4));
942: i = jl[k]; /* first row to be added to k_th row */
944: while (i < k) {
945: nexti = jl[i]; /* next row to be added to k_th row */
947: /* compute multiplier */
948: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
950: /* uik = -inv(Di)*U_bar(i,k): - ba[ili]*ba[i] */
951: diag = ba + i * 4;
952: u = ba + ili * 4;
953: uik[0] = -(diag[0] * u[0] + diag[2] * u[1]);
954: uik[1] = -(diag[1] * u[0] + diag[3] * u[1]);
955: uik[2] = -(diag[0] * u[2] + diag[2] * u[3]);
956: uik[3] = -(diag[1] * u[2] + diag[3] * u[3]);
958: /* update D(k) += -U(i,k)^T * U_bar(i,k): dk += uik*ba[ili] */
959: dk[0] += uik[0] * u[0] + uik[1] * u[1];
960: dk[1] += uik[2] * u[0] + uik[3] * u[1];
961: dk[2] += uik[0] * u[2] + uik[1] * u[3];
962: dk[3] += uik[2] * u[2] + uik[3] * u[3];
964: PetscCall(PetscLogFlops(16.0 * 2.0));
966: /* update -U(i,k): ba[ili] = uik */
967: PetscCall(PetscArraycpy(ba + ili * 4, uik, 4));
969: /* add multiple of row i to k-th row ... */
970: jmin = ili + 1;
971: jmax = bi[i + 1];
972: if (jmin < jmax) {
973: for (j = jmin; j < jmax; j++) {
974: /* rtmp += -U(i,k)^T * U_bar(i,j): rtmp[bj[j]] += uik*ba[j]; */
975: rtmp_ptr = rtmp + bj[j] * 4;
976: u = ba + j * 4;
977: rtmp_ptr[0] += uik[0] * u[0] + uik[1] * u[1];
978: rtmp_ptr[1] += uik[2] * u[0] + uik[3] * u[1];
979: rtmp_ptr[2] += uik[0] * u[2] + uik[1] * u[3];
980: rtmp_ptr[3] += uik[2] * u[2] + uik[3] * u[3];
981: }
982: PetscCall(PetscLogFlops(16.0 * (jmax - jmin)));
984: /* ... add i to row list for next nonzero entry */
985: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
986: j = bj[jmin];
987: jl[i] = jl[j];
988: jl[j] = i; /* update jl */
989: }
990: i = nexti;
991: }
993: /* save nonzero entries in k-th row of U ... */
995: /* invert diagonal block */
996: diag = ba + k * 4;
997: PetscCall(PetscArraycpy(diag, dk, 4));
998: PetscCall(PetscKernel_A_gets_inverse_A_2(diag, shift, allowzeropivot, &zeropivotdetected));
999: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1001: jmin = bi[k];
1002: jmax = bi[k + 1];
1003: if (jmin < jmax) {
1004: for (j = jmin; j < jmax; j++) {
1005: vj = bj[j]; /* block col. index of U */
1006: u = ba + j * 4;
1007: rtmp_ptr = rtmp + vj * 4;
1008: for (k1 = 0; k1 < 4; k1++) {
1009: *u++ = *rtmp_ptr;
1010: *rtmp_ptr++ = 0.0;
1011: }
1012: }
1014: /* ... add k to row list for first nonzero entry in k-th row */
1015: il[k] = jmin;
1016: i = bj[jmin];
1017: jl[k] = jl[i];
1018: jl[i] = k;
1019: }
1020: }
1022: PetscCall(PetscFree(rtmp));
1023: PetscCall(PetscFree2(il, jl));
1024: if (a->permute) PetscCall(PetscFree(aa));
1025: PetscCall(ISRestoreIndices(perm, &perm_ptr));
1027: C->ops->solve = MatSolve_SeqSBAIJ_2_inplace;
1028: C->ops->solvetranspose = MatSolve_SeqSBAIJ_2_inplace;
1029: C->assembled = PETSC_TRUE;
1030: C->preallocated = PETSC_TRUE;
1032: PetscCall(PetscLogFlops(1.3333 * 8 * b->mbs)); /* from inverting diagonal blocks */
1033: PetscFunctionReturn(PETSC_SUCCESS);
1034: }
1036: /*
1037: Version for when blocks are 2 by 2 Using natural ordering
1038: */
1039: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering(Mat C, Mat A, const MatFactorInfo *info)
1040: {
1041: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
1042: PetscInt i, j, mbs = a->mbs, *bi = b->i, *bj = b->j;
1043: PetscInt *ai, *aj, k, k1, jmin, jmax, *jl, *il, vj, nexti, ili;
1044: MatScalar *ba = b->a, *aa, *ap, dk[8], uik[8];
1045: MatScalar *u, *diag, *rtmp, *rtmp_ptr;
1046: PetscReal shift = info->shiftamount;
1047: PetscBool allowzeropivot, zeropivotdetected;
1049: PetscFunctionBegin;
1050: allowzeropivot = PetscNot(A->erroriffailure);
1052: /* initialization */
1053: /* il and jl record the first nonzero element in each row of the accessing
1054: window U(0:k, k:mbs-1).
1055: jl: list of rows to be added to uneliminated rows
1056: i>= k: jl(i) is the first row to be added to row i
1057: i< k: jl(i) is the row following row i in some list of rows
1058: jl(i) = mbs indicates the end of a list
1059: il(i): points to the first nonzero element in columns k,...,mbs-1 of
1060: row i of U */
1061: PetscCall(PetscCalloc1(4 * mbs, &rtmp));
1062: PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
1063: il[0] = 0;
1064: for (i = 0; i < mbs; i++) jl[i] = mbs;
1066: ai = a->i;
1067: aj = a->j;
1068: aa = a->a;
1070: /* for each row k */
1071: for (k = 0; k < mbs; k++) {
1072: /*initialize k-th row with elements nonzero in row k of A */
1073: jmin = ai[k];
1074: jmax = ai[k + 1];
1075: ap = aa + jmin * 4;
1076: for (j = jmin; j < jmax; j++) {
1077: vj = aj[j]; /* block col. index */
1078: rtmp_ptr = rtmp + vj * 4;
1079: for (i = 0; i < 4; i++) *rtmp_ptr++ = *ap++;
1080: }
1082: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
1083: PetscCall(PetscArraycpy(dk, rtmp + k * 4, 4));
1084: i = jl[k]; /* first row to be added to k_th row */
1086: while (i < k) {
1087: nexti = jl[i]; /* next row to be added to k_th row */
1089: /* compute multiplier */
1090: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1092: /* uik = -inv(Di)*U_bar(i,k): - ba[ili]*ba[i] */
1093: diag = ba + i * 4;
1094: u = ba + ili * 4;
1095: uik[0] = -(diag[0] * u[0] + diag[2] * u[1]);
1096: uik[1] = -(diag[1] * u[0] + diag[3] * u[1]);
1097: uik[2] = -(diag[0] * u[2] + diag[2] * u[3]);
1098: uik[3] = -(diag[1] * u[2] + diag[3] * u[3]);
1100: /* update D(k) += -U(i,k)^T * U_bar(i,k): dk += uik*ba[ili] */
1101: dk[0] += uik[0] * u[0] + uik[1] * u[1];
1102: dk[1] += uik[2] * u[0] + uik[3] * u[1];
1103: dk[2] += uik[0] * u[2] + uik[1] * u[3];
1104: dk[3] += uik[2] * u[2] + uik[3] * u[3];
1106: PetscCall(PetscLogFlops(16.0 * 2.0));
1108: /* update -U(i,k): ba[ili] = uik */
1109: PetscCall(PetscArraycpy(ba + ili * 4, uik, 4));
1111: /* add multiple of row i to k-th row ... */
1112: jmin = ili + 1;
1113: jmax = bi[i + 1];
1114: if (jmin < jmax) {
1115: for (j = jmin; j < jmax; j++) {
1116: /* rtmp += -U(i,k)^T * U_bar(i,j): rtmp[bj[j]] += uik*ba[j]; */
1117: rtmp_ptr = rtmp + bj[j] * 4;
1118: u = ba + j * 4;
1119: rtmp_ptr[0] += uik[0] * u[0] + uik[1] * u[1];
1120: rtmp_ptr[1] += uik[2] * u[0] + uik[3] * u[1];
1121: rtmp_ptr[2] += uik[0] * u[2] + uik[1] * u[3];
1122: rtmp_ptr[3] += uik[2] * u[2] + uik[3] * u[3];
1123: }
1124: PetscCall(PetscLogFlops(16.0 * (jmax - jmin)));
1126: /* ... add i to row list for next nonzero entry */
1127: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
1128: j = bj[jmin];
1129: jl[i] = jl[j];
1130: jl[j] = i; /* update jl */
1131: }
1132: i = nexti;
1133: }
1135: /* save nonzero entries in k-th row of U ... */
1137: /* invert diagonal block */
1138: diag = ba + k * 4;
1139: PetscCall(PetscArraycpy(diag, dk, 4));
1140: PetscCall(PetscKernel_A_gets_inverse_A_2(diag, shift, allowzeropivot, &zeropivotdetected));
1141: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1143: jmin = bi[k];
1144: jmax = bi[k + 1];
1145: if (jmin < jmax) {
1146: for (j = jmin; j < jmax; j++) {
1147: vj = bj[j]; /* block col. index of U */
1148: u = ba + j * 4;
1149: rtmp_ptr = rtmp + vj * 4;
1150: for (k1 = 0; k1 < 4; k1++) {
1151: *u++ = *rtmp_ptr;
1152: *rtmp_ptr++ = 0.0;
1153: }
1154: }
1156: /* ... add k to row list for first nonzero entry in k-th row */
1157: il[k] = jmin;
1158: i = bj[jmin];
1159: jl[k] = jl[i];
1160: jl[i] = k;
1161: }
1162: }
1164: PetscCall(PetscFree(rtmp));
1165: PetscCall(PetscFree2(il, jl));
1167: C->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1168: C->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1169: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1170: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1171: C->assembled = PETSC_TRUE;
1172: C->preallocated = PETSC_TRUE;
1174: PetscCall(PetscLogFlops(1.3333 * 8 * b->mbs)); /* from inverting diagonal blocks */
1175: PetscFunctionReturn(PETSC_SUCCESS);
1176: }
1178: /*
1179: Numeric U^T*D*U factorization for SBAIJ format.
1180: Version for blocks are 1 by 1.
1181: */
1182: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace(Mat C, Mat A, const MatFactorInfo *info)
1183: {
1184: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
1185: IS ip = b->row;
1186: const PetscInt *ai, *aj, *rip;
1187: PetscInt *a2anew, i, j, mbs = a->mbs, *bi = b->i, *bj = b->j, *bcol;
1188: PetscInt k, jmin, jmax, *jl, *il, col, nexti, ili, nz;
1189: MatScalar *rtmp, *ba = b->a, *bval, *aa, dk, uikdi;
1190: PetscReal rs;
1191: FactorShiftCtx sctx;
1193: PetscFunctionBegin;
1194: /* MatPivotSetUp(): initialize shift context sctx */
1195: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
1197: PetscCall(ISGetIndices(ip, &rip));
1198: if (!a->permute) {
1199: ai = a->i;
1200: aj = a->j;
1201: aa = a->a;
1202: } else {
1203: ai = a->inew;
1204: aj = a->jnew;
1205: nz = ai[mbs];
1206: PetscCall(PetscMalloc1(nz, &aa));
1207: a2anew = a->a2anew;
1208: bval = a->a;
1209: for (j = 0; j < nz; j++) aa[a2anew[j]] = *(bval++);
1210: }
1212: /* initialization */
1213: /* il and jl record the first nonzero element in each row of the accessing
1214: window U(0:k, k:mbs-1).
1215: jl: list of rows to be added to uneliminated rows
1216: i>= k: jl(i) is the first row to be added to row i
1217: i< k: jl(i) is the row following row i in some list of rows
1218: jl(i) = mbs indicates the end of a list
1219: il(i): points to the first nonzero element in columns k,...,mbs-1 of
1220: row i of U */
1221: PetscCall(PetscMalloc3(mbs, &rtmp, mbs, &il, mbs, &jl));
1223: do {
1224: sctx.newshift = PETSC_FALSE;
1225: il[0] = 0;
1226: for (i = 0; i < mbs; i++) {
1227: rtmp[i] = 0.0;
1228: jl[i] = mbs;
1229: }
1231: for (k = 0; k < mbs; k++) {
1232: /*initialize k-th row by the perm[k]-th row of A */
1233: jmin = ai[rip[k]];
1234: jmax = ai[rip[k] + 1];
1235: bval = ba + bi[k];
1236: for (j = jmin; j < jmax; j++) {
1237: col = rip[aj[j]];
1238: rtmp[col] = aa[j];
1239: *bval++ = 0.0; /* for in-place factorization */
1240: }
1242: /* shift the diagonal of the matrix */
1243: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
1245: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1246: dk = rtmp[k];
1247: i = jl[k]; /* first row to be added to k_th row */
1249: while (i < k) {
1250: nexti = jl[i]; /* next row to be added to k_th row */
1252: /* compute multiplier, update diag(k) and U(i,k) */
1253: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1254: uikdi = -ba[ili] * ba[bi[i]]; /* diagonal(k) */
1255: dk += uikdi * ba[ili];
1256: ba[ili] = uikdi; /* -U(i,k) */
1258: /* add multiple of row i to k-th row */
1259: jmin = ili + 1;
1260: jmax = bi[i + 1];
1261: if (jmin < jmax) {
1262: for (j = jmin; j < jmax; j++) rtmp[bj[j]] += uikdi * ba[j];
1263: PetscCall(PetscLogFlops(2.0 * (jmax - jmin)));
1265: /* update il and jl for row i */
1266: il[i] = jmin;
1267: j = bj[jmin];
1268: jl[i] = jl[j];
1269: jl[j] = i;
1270: }
1271: i = nexti;
1272: }
1274: /* shift the diagonals when zero pivot is detected */
1275: /* compute rs=sum of abs(off-diagonal) */
1276: rs = 0.0;
1277: jmin = bi[k] + 1;
1278: nz = bi[k + 1] - jmin;
1279: if (nz) {
1280: bcol = bj + jmin;
1281: while (nz--) {
1282: rs += PetscAbsScalar(rtmp[*bcol]);
1283: bcol++;
1284: }
1285: }
1287: sctx.rs = rs;
1288: sctx.pv = dk;
1289: PetscCall(MatPivotCheck(C, A, info, &sctx, k));
1290: if (sctx.newshift) break; /* sctx.shift_amount is updated */
1291: dk = sctx.pv;
1293: /* copy data into U(k,:) */
1294: ba[bi[k]] = 1.0 / dk; /* U(k,k) */
1295: jmin = bi[k] + 1;
1296: jmax = bi[k + 1];
1297: if (jmin < jmax) {
1298: for (j = jmin; j < jmax; j++) {
1299: col = bj[j];
1300: ba[j] = rtmp[col];
1301: rtmp[col] = 0.0;
1302: }
1303: /* add the k-th row into il and jl */
1304: il[k] = jmin;
1305: i = bj[jmin];
1306: jl[k] = jl[i];
1307: jl[i] = k;
1308: }
1309: }
1310: } while (sctx.newshift);
1311: PetscCall(PetscFree3(rtmp, il, jl));
1312: if (a->permute) PetscCall(PetscFree(aa));
1314: PetscCall(ISRestoreIndices(ip, &rip));
1316: C->ops->solve = MatSolve_SeqSBAIJ_1_inplace;
1317: C->ops->solves = MatSolves_SeqSBAIJ_1_inplace;
1318: C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
1319: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_inplace;
1320: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_inplace;
1321: C->assembled = PETSC_TRUE;
1322: C->preallocated = PETSC_TRUE;
1324: PetscCall(PetscLogFlops(C->rmap->N));
1325: if (sctx.nshift) {
1326: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1327: PetscCall(PetscInfo(A, "number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1328: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1329: PetscCall(PetscInfo(A, "number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1330: }
1331: }
1332: PetscFunctionReturn(PETSC_SUCCESS);
1333: }
1335: /*
1336: Version for when blocks are 1 by 1 Using natural ordering under new datastructure
1337: Modified from MatCholeskyFactorNumeric_SeqAIJ()
1338: */
1339: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info)
1340: {
1341: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1342: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data;
1343: PetscInt i, j, mbs = A->rmap->n, *bi = b->i, *bj = b->j, *bdiag = b->diag, *bjtmp;
1344: PetscInt *ai = a->i, *aj = a->j, *ajtmp;
1345: PetscInt k, jmin, jmax, *c2r, *il, col, nexti, ili, nz;
1346: MatScalar *rtmp, *ba = b->a, *bval, *aa = a->a, dk, uikdi;
1347: FactorShiftCtx sctx;
1348: PetscReal rs;
1349: MatScalar d, *v;
1350: const PetscInt *adiag;
1352: PetscFunctionBegin;
1353: PetscCall(MatGetDiagonalMarkers_SeqSBAIJ(A, &adiag, NULL));
1354: PetscCall(PetscMalloc3(mbs, &rtmp, mbs, &il, mbs, &c2r));
1356: /* MatPivotSetUp(): initialize shift context sctx */
1357: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
1359: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
1360: sctx.shift_top = info->zeropivot;
1362: PetscCall(PetscArrayzero(rtmp, mbs));
1364: for (i = 0; i < mbs; i++) {
1365: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
1366: d = aa[adiag[i]];
1367: rtmp[i] += -PetscRealPart(d); /* diagonal entry */
1368: ajtmp = aj + ai[i] + 1; /* exclude diagonal */
1369: v = aa + ai[i] + 1;
1370: nz = ai[i + 1] - ai[i] - 1;
1371: for (j = 0; j < nz; j++) {
1372: rtmp[i] += PetscAbsScalar(v[j]);
1373: rtmp[ajtmp[j]] += PetscAbsScalar(v[j]);
1374: }
1375: if (PetscRealPart(rtmp[i]) > sctx.shift_top) sctx.shift_top = PetscRealPart(rtmp[i]);
1376: }
1377: sctx.shift_top *= 1.1;
1378: sctx.nshift_max = 5;
1379: sctx.shift_lo = 0.;
1380: sctx.shift_hi = 1.;
1381: }
1383: /* allocate working arrays
1384: c2r: linked list, keep track of pivot rows for a given column. c2r[col]: head of the list for a given col
1385: 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
1386: */
1387: do {
1388: sctx.newshift = PETSC_FALSE;
1390: for (i = 0; i < mbs; i++) c2r[i] = mbs;
1391: if (mbs) il[0] = 0;
1393: for (k = 0; k < mbs; k++) {
1394: /* zero rtmp */
1395: nz = bi[k + 1] - bi[k];
1396: bjtmp = bj + bi[k];
1397: for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;
1399: /* load in initial unfactored row */
1400: bval = ba + bi[k];
1401: jmin = ai[k];
1402: jmax = ai[k + 1];
1403: for (j = jmin; j < jmax; j++) {
1404: col = aj[j];
1405: rtmp[col] = aa[j];
1406: *bval++ = 0.0; /* for in-place factorization */
1407: }
1408: /* shift the diagonal of the matrix: ZeropivotApply() */
1409: rtmp[k] += sctx.shift_amount; /* shift the diagonal of the matrix */
1411: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1412: dk = rtmp[k];
1413: i = c2r[k]; /* first row to be added to k_th row */
1415: while (i < k) {
1416: nexti = c2r[i]; /* next row to be added to k_th row */
1418: /* compute multiplier, update diag(k) and U(i,k) */
1419: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1420: uikdi = -ba[ili] * ba[bdiag[i]]; /* diagonal(k) */
1421: dk += uikdi * ba[ili]; /* update diag[k] */
1422: ba[ili] = uikdi; /* -U(i,k) */
1424: /* add multiple of row i to k-th row */
1425: jmin = ili + 1;
1426: jmax = bi[i + 1];
1427: if (jmin < jmax) {
1428: for (j = jmin; j < jmax; j++) rtmp[bj[j]] += uikdi * ba[j];
1429: /* update il and c2r for row i */
1430: il[i] = jmin;
1431: j = bj[jmin];
1432: c2r[i] = c2r[j];
1433: c2r[j] = i;
1434: }
1435: i = nexti;
1436: }
1438: /* copy data into U(k,:) */
1439: rs = 0.0;
1440: jmin = bi[k];
1441: jmax = bi[k + 1] - 1;
1442: if (jmin < jmax) {
1443: for (j = jmin; j < jmax; j++) {
1444: col = bj[j];
1445: ba[j] = rtmp[col];
1446: rs += PetscAbsScalar(ba[j]);
1447: }
1448: /* add the k-th row into il and c2r */
1449: il[k] = jmin;
1450: i = bj[jmin];
1451: c2r[k] = c2r[i];
1452: c2r[i] = k;
1453: }
1455: sctx.rs = rs;
1456: sctx.pv = dk;
1457: PetscCall(MatPivotCheck(B, A, info, &sctx, k));
1458: if (sctx.newshift) break;
1459: dk = sctx.pv;
1461: ba[bdiag[k]] = 1.0 / dk; /* U(k,k) */
1462: }
1463: } while (sctx.newshift);
1465: PetscCall(PetscFree3(rtmp, il, c2r));
1467: B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1468: B->ops->solves = MatSolves_SeqSBAIJ_1;
1469: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1470: B->ops->matsolve = MatMatSolve_SeqSBAIJ_1_NaturalOrdering;
1471: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
1472: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;
1474: B->assembled = PETSC_TRUE;
1475: B->preallocated = PETSC_TRUE;
1477: PetscCall(PetscLogFlops(B->rmap->n));
1479: /* MatPivotView() */
1480: if (sctx.nshift) {
1481: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1482: 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));
1483: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1484: PetscCall(PetscInfo(A, "number of shift_nz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1485: } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
1486: PetscCall(PetscInfo(A, "number of shift_inblocks applied %" PetscInt_FMT ", each shift_amount %g\n", sctx.nshift, (double)info->shiftamount));
1487: }
1488: }
1489: PetscFunctionReturn(PETSC_SUCCESS);
1490: }
1492: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace(Mat C, Mat A, const MatFactorInfo *info)
1493: {
1494: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
1495: PetscInt i, j, mbs = a->mbs;
1496: PetscInt *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
1497: PetscInt k, jmin, *jl, *il, nexti, ili, *acol, *bcol, nz;
1498: MatScalar *rtmp, *ba = b->a, *aa = a->a, dk, uikdi, *aval, *bval;
1499: PetscReal rs;
1500: FactorShiftCtx sctx;
1502: PetscFunctionBegin;
1503: /* MatPivotSetUp(): initialize shift context sctx */
1504: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
1506: /* initialization */
1507: /* il and jl record the first nonzero element in each row of the accessing
1508: window U(0:k, k:mbs-1).
1509: jl: list of rows to be added to uneliminated rows
1510: i>= k: jl(i) is the first row to be added to row i
1511: i< k: jl(i) is the row following row i in some list of rows
1512: jl(i) = mbs indicates the end of a list
1513: il(i): points to the first nonzero element in U(i,k:mbs-1)
1514: */
1515: PetscCall(PetscMalloc1(mbs, &rtmp));
1516: PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
1518: do {
1519: sctx.newshift = PETSC_FALSE;
1520: il[0] = 0;
1521: for (i = 0; i < mbs; i++) {
1522: rtmp[i] = 0.0;
1523: jl[i] = mbs;
1524: }
1526: for (k = 0; k < mbs; k++) {
1527: /*initialize k-th row with elements nonzero in row perm(k) of A */
1528: nz = ai[k + 1] - ai[k];
1529: acol = aj + ai[k];
1530: aval = aa + ai[k];
1531: bval = ba + bi[k];
1532: while (nz--) {
1533: rtmp[*acol++] = *aval++;
1534: *bval++ = 0.0; /* for in-place factorization */
1535: }
1537: /* shift the diagonal of the matrix */
1538: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
1540: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1541: dk = rtmp[k];
1542: i = jl[k]; /* first row to be added to k_th row */
1544: while (i < k) {
1545: nexti = jl[i]; /* next row to be added to k_th row */
1546: /* compute multiplier, update D(k) and U(i,k) */
1547: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1548: uikdi = -ba[ili] * ba[bi[i]];
1549: dk += uikdi * ba[ili];
1550: ba[ili] = uikdi; /* -U(i,k) */
1552: /* add multiple of row i to k-th row ... */
1553: jmin = ili + 1;
1554: nz = bi[i + 1] - jmin;
1555: if (nz > 0) {
1556: bcol = bj + jmin;
1557: bval = ba + jmin;
1558: PetscCall(PetscLogFlops(2.0 * nz));
1559: while (nz--) rtmp[*bcol++] += uikdi * (*bval++);
1561: /* update il and jl for i-th row */
1562: il[i] = jmin;
1563: j = bj[jmin];
1564: jl[i] = jl[j];
1565: jl[j] = i;
1566: }
1567: i = nexti;
1568: }
1570: /* shift the diagonals when zero pivot is detected */
1571: /* compute rs=sum of abs(off-diagonal) */
1572: rs = 0.0;
1573: jmin = bi[k] + 1;
1574: nz = bi[k + 1] - jmin;
1575: if (nz) {
1576: bcol = bj + jmin;
1577: while (nz--) {
1578: rs += PetscAbsScalar(rtmp[*bcol]);
1579: bcol++;
1580: }
1581: }
1583: sctx.rs = rs;
1584: sctx.pv = dk;
1585: PetscCall(MatPivotCheck(C, A, info, &sctx, k));
1586: if (sctx.newshift) break; /* sctx.shift_amount is updated */
1587: dk = sctx.pv;
1589: /* copy data into U(k,:) */
1590: ba[bi[k]] = 1.0 / dk;
1591: jmin = bi[k] + 1;
1592: nz = bi[k + 1] - jmin;
1593: if (nz) {
1594: bcol = bj + jmin;
1595: bval = ba + jmin;
1596: while (nz--) {
1597: *bval++ = rtmp[*bcol];
1598: rtmp[*bcol++] = 0.0;
1599: }
1600: /* add k-th row into il and jl */
1601: il[k] = jmin;
1602: i = bj[jmin];
1603: jl[k] = jl[i];
1604: jl[i] = k;
1605: }
1606: } /* end of for (k = 0; k<mbs; k++) */
1607: } while (sctx.newshift);
1608: PetscCall(PetscFree(rtmp));
1609: PetscCall(PetscFree2(il, jl));
1611: C->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1612: C->ops->solves = MatSolves_SeqSBAIJ_1_inplace;
1613: C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1614: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1615: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1617: C->assembled = PETSC_TRUE;
1618: C->preallocated = PETSC_TRUE;
1620: PetscCall(PetscLogFlops(C->rmap->N));
1621: if (sctx.nshift) {
1622: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1623: PetscCall(PetscInfo(A, "number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1624: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1625: PetscCall(PetscInfo(A, "number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1626: }
1627: }
1628: PetscFunctionReturn(PETSC_SUCCESS);
1629: }
1631: PetscErrorCode MatCholeskyFactor_SeqSBAIJ(Mat A, IS perm, const MatFactorInfo *info)
1632: {
1633: Mat C;
1635: PetscFunctionBegin;
1636: PetscCall(MatGetFactor(A, "petsc", MAT_FACTOR_CHOLESKY, &C));
1637: PetscCall(MatCholeskyFactorSymbolic(C, A, perm, info));
1638: PetscCall(MatCholeskyFactorNumeric(C, A, info));
1640: A->ops->solve = C->ops->solve;
1641: A->ops->solvetranspose = C->ops->solvetranspose;
1643: PetscCall(MatHeaderMerge(A, &C));
1644: PetscFunctionReturn(PETSC_SUCCESS);
1645: }