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 <petscbt.h>
217: #include <../src/mat/utils/freespace.h>
218: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
219: {
220: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
221: Mat_SeqSBAIJ *b;
222: PetscBool perm_identity;
223: PetscReal fill = info->fill;
224: const PetscInt *rip, *ai = a->i, *aj = a->j;
225: PetscInt i, mbs = a->mbs, bs = A->rmap->bs, reallocs = 0, prow;
226: PetscInt *jl, jmin, jmax, nzk, *ui, k, j, *il, nextprow;
227: PetscInt nlnk, *lnk, ncols, *cols, *uj, **ui_ptr, *uj_ptr, *udiag;
228: PetscFreeSpaceList free_space = NULL, current_space = NULL;
229: PetscBT lnkbt;
230: PetscBool diagDense;
232: PetscFunctionBegin;
233: 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);
234: PetscCall(MatGetDiagonalMarkers_SeqSBAIJ(A, NULL, &diagDense));
235: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry");
236: if (bs > 1) {
237: PetscCall(MatCholeskyFactorSymbolic_SeqSBAIJ_inplace(fact, A, perm, info));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: /* check whether perm is the identity mapping */
242: PetscCall(ISIdentity(perm, &perm_identity));
243: PetscCheck(perm_identity, PETSC_COMM_SELF, PETSC_ERR_SUP, "Matrix reordering is not supported for sbaij matrix. Use aij format");
244: a->permute = PETSC_FALSE;
245: PetscCall(ISGetIndices(perm, &rip));
247: /* initialization */
248: PetscCall(PetscShmgetAllocateArray(mbs + 1, sizeof(PetscInt), (void **)&ui));
249: PetscCall(PetscMalloc1(mbs + 1, &udiag));
250: ui[0] = 0;
252: /* jl: linked list for storing indices of the pivot rows
253: il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
254: PetscCall(PetscMalloc4(mbs, &ui_ptr, mbs, &il, mbs, &jl, mbs, &cols));
255: for (i = 0; i < mbs; i++) {
256: jl[i] = mbs;
257: il[i] = 0;
258: }
260: /* create and initialize a linked list for storing column indices of the active row k */
261: nlnk = mbs + 1;
262: PetscCall(PetscLLCreate(mbs, mbs, nlnk, lnk, lnkbt));
264: /* initial FreeSpace size is fill*(ai[mbs]+1) */
265: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, ai[mbs] + 1), &free_space));
266: current_space = free_space;
268: for (k = 0; k < mbs; k++) { /* for each active row k */
269: /* initialize lnk by the column indices of row rip[k] of A */
270: nzk = 0;
271: ncols = ai[k + 1] - ai[k];
272: PetscCheck(ncols, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Empty row %" PetscInt_FMT " in matrix ", k);
273: for (j = 0; j < ncols; j++) {
274: i = *(aj + ai[k] + j);
275: cols[j] = i;
276: }
277: PetscCall(PetscLLAdd(ncols, cols, mbs, &nlnk, lnk, lnkbt));
278: nzk += nlnk;
280: /* update lnk by computing fill-in for each pivot row to be merged in */
281: prow = jl[k]; /* 1st pivot row */
283: while (prow < k) {
284: nextprow = jl[prow];
285: /* merge prow into k-th row */
286: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
287: jmax = ui[prow + 1];
288: ncols = jmax - jmin;
289: uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
290: PetscCall(PetscLLAddSorted(ncols, uj_ptr, mbs, &nlnk, lnk, lnkbt));
291: nzk += nlnk;
293: /* update il and jl for prow */
294: if (jmin < jmax) {
295: il[prow] = jmin;
296: j = *uj_ptr;
297: jl[prow] = jl[j];
298: jl[j] = prow;
299: }
300: prow = nextprow;
301: }
303: /* if free space is not available, make more free space */
304: if (current_space->local_remaining < nzk) {
305: i = mbs - k + 1; /* num of unfactored rows */
306: i = PetscIntMultTruncate(i, PetscMin(nzk, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
307: PetscCall(PetscFreeSpaceGet(i, ¤t_space));
308: reallocs++;
309: }
311: /* copy data into free space, then initialize lnk */
312: PetscCall(PetscLLClean(mbs, mbs, nzk, lnk, current_space->array, lnkbt));
314: /* add the k-th row into il and jl */
315: if (nzk > 1) {
316: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
317: jl[k] = jl[i];
318: jl[i] = k;
319: il[k] = ui[k] + 1;
320: }
321: ui_ptr[k] = current_space->array;
323: current_space->array += nzk;
324: current_space->local_used += nzk;
325: current_space->local_remaining -= nzk;
327: ui[k + 1] = ui[k] + nzk;
328: }
330: PetscCall(ISRestoreIndices(perm, &rip));
331: PetscCall(PetscFree4(ui_ptr, il, jl, cols));
333: /* destroy list of free space and other temporary array(s) */
334: PetscCall(PetscShmgetAllocateArray(ui[mbs], sizeof(PetscInt), (void **)&uj));
335: PetscCall(PetscFreeSpaceContiguous_Cholesky(&free_space, uj, mbs, ui, udiag)); /* store matrix factor */
336: PetscCall(PetscLLDestroy(lnk, lnkbt));
338: /* put together the new matrix in MATSEQSBAIJ format */
339: PetscCall(MatSeqSBAIJSetPreallocation(fact, bs, MAT_SKIP_ALLOCATION, NULL));
341: b = (Mat_SeqSBAIJ *)fact->data;
342: b->free_ij = PETSC_TRUE;
343: PetscCall(PetscShmgetAllocateArray(ui[mbs], sizeof(PetscScalar), (void **)&b->a));
344: b->free_a = PETSC_TRUE;
345: b->j = uj;
346: b->i = ui;
347: b->diag = udiag;
348: b->ilen = NULL;
349: b->imax = NULL;
350: b->row = perm;
351: b->icol = perm;
353: PetscCall(PetscObjectReference((PetscObject)perm));
354: PetscCall(PetscObjectReference((PetscObject)perm));
356: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
358: PetscCall(PetscMalloc1(mbs + 1, &b->solve_work));
360: b->maxnz = b->nz = ui[mbs];
362: fact->info.factor_mallocs = reallocs;
363: fact->info.fill_ratio_given = fill;
364: if (ai[mbs] != 0) {
365: fact->info.fill_ratio_needed = ((PetscReal)ui[mbs]) / ai[mbs];
366: } else {
367: fact->info.fill_ratio_needed = 0.0;
368: }
369: #if PetscDefined(USE_INFO)
370: if (ai[mbs] != 0) {
371: PetscReal af = fact->info.fill_ratio_needed;
372: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
373: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
374: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
375: } else PetscCall(PetscInfo(A, "Empty matrix\n"));
376: #endif
377: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering;
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ_inplace(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
382: {
383: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
384: Mat_SeqSBAIJ *b;
385: PetscBool perm_identity;
386: PetscReal fill = info->fill;
387: const PetscInt *rip, *ai, *aj;
388: PetscInt i, mbs = a->mbs, bs = A->rmap->bs, reallocs = 0, prow;
389: PetscInt *jl, jmin, jmax, nzk, *ui, k, j, *il, nextprow;
390: PetscInt nlnk, *lnk, ncols, *cols, *uj, **ui_ptr, *uj_ptr;
391: PetscFreeSpaceList free_space = NULL, current_space = NULL;
392: PetscBT lnkbt;
393: PetscBool diagDense;
395: PetscFunctionBegin;
396: PetscCall(MatGetDiagonalMarkers_SeqSBAIJ(A, NULL, &diagDense));
397: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry");
399: /*
400: This code originally uses Modified Sparse Row (MSR) storage
401: (see page 85, "Iterative Methods ..." by Saad) for the output matrix B - bad choice!
402: Then it is rewritten so the factor B takes seqsbaij format. However the associated
403: MatCholeskyFactorNumeric_() have not been modified for the cases of bs>1 or !perm_identity,
404: thus the original code in MSR format is still used for these cases.
405: The code below should replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR() whenever
406: MatCholeskyFactorNumeric_() is modified for using sbaij symbolic factor.
407: */
408: if (bs > 1) {
409: PetscCall(MatCholeskyFactorSymbolic_SeqSBAIJ_MSR(fact, A, perm, info));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /* check whether perm is the identity mapping */
414: PetscCall(ISIdentity(perm, &perm_identity));
415: PetscCheck(perm_identity, PETSC_COMM_SELF, PETSC_ERR_SUP, "Matrix reordering is not supported for sbaij matrix. Use aij format");
416: a->permute = PETSC_FALSE;
417: ai = a->i;
418: aj = a->j;
419: PetscCall(ISGetIndices(perm, &rip));
421: /* initialization */
422: PetscCall(PetscShmgetAllocateArray(mbs + 1, sizeof(PetscInt), (void **)&ui));
423: ui[0] = 0;
425: /* jl: linked list for storing indices of the pivot rows
426: il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
427: PetscCall(PetscMalloc4(mbs, &ui_ptr, mbs, &il, mbs, &jl, mbs, &cols));
428: for (i = 0; i < mbs; i++) {
429: jl[i] = mbs;
430: il[i] = 0;
431: }
433: /* create and initialize a linked list for storing column indices of the active row k */
434: nlnk = mbs + 1;
435: PetscCall(PetscLLCreate(mbs, mbs, nlnk, lnk, lnkbt));
437: /* initial FreeSpace size is fill*(ai[mbs]+1) */
438: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, ai[mbs] + 1), &free_space));
439: current_space = free_space;
441: for (k = 0; k < mbs; k++) { /* for each active row k */
442: /* initialize lnk by the column indices of row rip[k] of A */
443: nzk = 0;
444: ncols = ai[rip[k] + 1] - ai[rip[k]];
445: for (j = 0; j < ncols; j++) {
446: i = *(aj + ai[rip[k]] + j);
447: cols[j] = rip[i];
448: }
449: PetscCall(PetscLLAdd(ncols, cols, mbs, &nlnk, lnk, lnkbt));
450: nzk += nlnk;
452: /* update lnk by computing fill-in for each pivot row to be merged in */
453: prow = jl[k]; /* 1st pivot row */
455: while (prow < k) {
456: nextprow = jl[prow];
457: /* merge prow into k-th row */
458: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
459: jmax = ui[prow + 1];
460: ncols = jmax - jmin;
461: uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
462: PetscCall(PetscLLAddSorted(ncols, uj_ptr, mbs, &nlnk, lnk, lnkbt));
463: nzk += nlnk;
465: /* update il and jl for prow */
466: if (jmin < jmax) {
467: il[prow] = jmin;
469: j = *uj_ptr;
470: jl[prow] = jl[j];
471: jl[j] = prow;
472: }
473: prow = nextprow;
474: }
476: /* if free space is not available, make more free space */
477: if (current_space->local_remaining < nzk) {
478: i = mbs - k + 1; /* num of unfactored rows */
479: i = PetscMin(PetscIntMultTruncate(i, nzk), PetscIntMultTruncate(i, i - 1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
480: PetscCall(PetscFreeSpaceGet(i, ¤t_space));
481: reallocs++;
482: }
484: /* copy data into free space, then initialize lnk */
485: PetscCall(PetscLLClean(mbs, mbs, nzk, lnk, current_space->array, lnkbt));
487: /* add the k-th row into il and jl */
488: if (nzk - 1 > 0) {
489: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
490: jl[k] = jl[i];
491: jl[i] = k;
492: il[k] = ui[k] + 1;
493: }
494: ui_ptr[k] = current_space->array;
496: current_space->array += nzk;
497: current_space->local_used += nzk;
498: current_space->local_remaining -= nzk;
500: ui[k + 1] = ui[k] + nzk;
501: }
503: PetscCall(ISRestoreIndices(perm, &rip));
504: PetscCall(PetscFree4(ui_ptr, il, jl, cols));
506: /* destroy list of free space and other temporary array(s) */
507: PetscCall(PetscShmgetAllocateArray(ui[mbs] + 1, sizeof(PetscInt), (void **)&uj));
508: PetscCall(PetscFreeSpaceContiguous(&free_space, uj));
509: PetscCall(PetscLLDestroy(lnk, lnkbt));
511: /* put together the new matrix in MATSEQSBAIJ format */
512: PetscCall(MatSeqSBAIJSetPreallocation(fact, bs, MAT_SKIP_ALLOCATION, NULL));
514: b = (Mat_SeqSBAIJ *)fact->data;
515: PetscCall(PetscShmgetAllocateArray(ui[mbs] + 1, sizeof(PetscScalar), (void **)&b->a));
516: b->free_a = PETSC_TRUE;
517: b->free_ij = PETSC_TRUE;
518: b->j = uj;
519: b->i = ui;
520: b->diag = NULL;
521: b->ilen = NULL;
522: b->imax = NULL;
523: b->row = perm;
525: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
527: PetscCall(PetscObjectReference((PetscObject)perm));
528: b->icol = perm;
529: PetscCall(PetscObjectReference((PetscObject)perm));
530: PetscCall(PetscMalloc1(mbs + 1, &b->solve_work));
531: b->maxnz = b->nz = ui[mbs];
533: fact->info.factor_mallocs = reallocs;
534: fact->info.fill_ratio_given = fill;
535: if (ai[mbs] != 0) {
536: fact->info.fill_ratio_needed = ((PetscReal)ui[mbs]) / ai[mbs];
537: } else {
538: fact->info.fill_ratio_needed = 0.0;
539: }
540: #if PetscDefined(USE_INFO)
541: if (ai[mbs] != 0) {
542: PetscReal af = fact->info.fill_ratio_needed;
543: PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)fill, (double)af));
544: PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
545: PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g) for best performance.\n", (double)af));
546: } else PetscCall(PetscInfo(A, "Empty matrix\n"));
547: #endif
548: PetscCall(MatSeqSBAIJSetNumericFactorization_inplace(fact, perm_identity));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat C, Mat A, const MatFactorInfo *info)
553: {
554: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
555: IS perm = b->row;
556: const PetscInt *ai, *aj, *perm_ptr, mbs = a->mbs, *bi = b->i, *bj = b->j;
557: PetscInt i, j;
558: PetscInt *a2anew, k, k1, jmin, jmax, *jl, *il, vj, nexti, ili;
559: PetscInt bs = A->rmap->bs, bs2 = a->bs2;
560: MatScalar *ba = b->a, *aa, *ap, *dk, *uik;
561: MatScalar *u, *diag, *rtmp, *rtmp_ptr;
562: MatScalar *work;
563: PetscInt *pivots;
564: PetscBool allowzeropivot, zeropivotdetected;
566: PetscFunctionBegin;
567: /* initialization */
568: PetscCall(PetscCalloc1(bs2 * mbs, &rtmp));
569: PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
570: allowzeropivot = PetscNot(A->erroriffailure);
572: il[0] = 0;
573: for (i = 0; i < mbs; i++) jl[i] = mbs;
575: PetscCall(PetscMalloc3(bs2, &dk, bs2, &uik, bs, &work));
576: PetscCall(PetscMalloc1(bs, &pivots));
578: PetscCall(ISGetIndices(perm, &perm_ptr));
580: /* check permutation */
581: if (!a->permute) {
582: ai = a->i;
583: aj = a->j;
584: aa = a->a;
585: } else {
586: ai = a->inew;
587: aj = a->jnew;
588: PetscCall(PetscMalloc1(bs2 * ai[mbs], &aa));
589: PetscCall(PetscArraycpy(aa, a->a, bs2 * ai[mbs]));
590: PetscCall(PetscMalloc1(ai[mbs], &a2anew));
591: PetscCall(PetscArraycpy(a2anew, a->a2anew, ai[mbs]));
593: for (i = 0; i < mbs; i++) {
594: jmin = ai[i];
595: jmax = ai[i + 1];
596: for (j = jmin; j < jmax; j++) {
597: while (a2anew[j] != j) {
598: k = a2anew[j];
599: a2anew[j] = a2anew[k];
600: a2anew[k] = k;
601: for (k1 = 0; k1 < bs2; k1++) {
602: dk[k1] = aa[k * bs2 + k1];
603: aa[k * bs2 + k1] = aa[j * bs2 + k1];
604: aa[j * bs2 + k1] = dk[k1];
605: }
606: }
607: /* transform column-oriented blocks that lie in the lower triangle to row-oriented blocks */
608: if (i > aj[j]) {
609: ap = aa + j * bs2; /* ptr to the beginning of j-th block of aa */
610: for (k = 0; k < bs2; k++) dk[k] = ap[k]; /* dk <- j-th block of aa */
611: for (k = 0; k < bs; k++) { /* j-th block of aa <- dk^T */
612: for (k1 = 0; k1 < bs; k1++) *ap++ = dk[k + bs * k1];
613: }
614: }
615: }
616: }
617: PetscCall(PetscFree(a2anew));
618: }
620: /* for each row k */
621: for (k = 0; k < mbs; k++) {
622: /*initialize k-th row with elements nonzero in row perm(k) of A */
623: jmin = ai[perm_ptr[k]];
624: jmax = ai[perm_ptr[k] + 1];
626: ap = aa + jmin * bs2;
627: for (j = jmin; j < jmax; j++) {
628: vj = perm_ptr[aj[j]]; /* block col. index */
629: rtmp_ptr = rtmp + vj * bs2;
630: for (i = 0; i < bs2; i++) *rtmp_ptr++ = *ap++;
631: }
633: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
634: PetscCall(PetscArraycpy(dk, rtmp + k * bs2, bs2));
635: i = jl[k]; /* first row to be added to k_th row */
637: while (i < k) {
638: nexti = jl[i]; /* next row to be added to k_th row */
640: /* compute multiplier */
641: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
643: /* uik = -inv(Di)*U_bar(i,k) */
644: diag = ba + i * bs2;
645: u = ba + ili * bs2;
646: PetscCall(PetscArrayzero(uik, bs2));
647: PetscKernel_A_gets_A_minus_B_times_C(bs, uik, diag, u);
649: /* update D(k) += -U(i,k)^T * U_bar(i,k) */
650: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs, dk, uik, u);
651: PetscCall(PetscLogFlops(4.0 * bs * bs2));
653: /* update -U(i,k) */
654: PetscCall(PetscArraycpy(ba + ili * bs2, uik, bs2));
656: /* add multiple of row i to k-th row ... */
657: jmin = ili + 1;
658: jmax = bi[i + 1];
659: if (jmin < jmax) {
660: for (j = jmin; j < jmax; j++) {
661: /* rtmp += -U(i,k)^T * U_bar(i,j) */
662: rtmp_ptr = rtmp + bj[j] * bs2;
663: u = ba + j * bs2;
664: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs, rtmp_ptr, uik, u);
665: }
666: PetscCall(PetscLogFlops(2.0 * bs * bs2 * (jmax - jmin)));
668: /* ... add i to row list for next nonzero entry */
669: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
670: j = bj[jmin];
671: jl[i] = jl[j];
672: jl[j] = i; /* update jl */
673: }
674: i = nexti;
675: }
677: /* save nonzero entries in k-th row of U ... */
679: /* invert diagonal block */
680: diag = ba + k * bs2;
681: PetscCall(PetscArraycpy(diag, dk, bs2));
683: PetscCall(PetscKernel_A_gets_inverse_A(bs, diag, pivots, work, allowzeropivot, &zeropivotdetected));
684: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
686: jmin = bi[k];
687: jmax = bi[k + 1];
688: if (jmin < jmax) {
689: for (j = jmin; j < jmax; j++) {
690: vj = bj[j]; /* block col. index of U */
691: u = ba + j * bs2;
692: rtmp_ptr = rtmp + vj * bs2;
693: for (k1 = 0; k1 < bs2; k1++) {
694: *u++ = *rtmp_ptr;
695: *rtmp_ptr++ = 0.0;
696: }
697: }
699: /* ... add k to row list for first nonzero entry in k-th row */
700: il[k] = jmin;
701: i = bj[jmin];
702: jl[k] = jl[i];
703: jl[i] = k;
704: }
705: }
707: PetscCall(PetscFree(rtmp));
708: PetscCall(PetscFree2(il, jl));
709: PetscCall(PetscFree3(dk, uik, work));
710: PetscCall(PetscFree(pivots));
711: if (a->permute) PetscCall(PetscFree(aa));
713: PetscCall(ISRestoreIndices(perm, &perm_ptr));
715: C->ops->solve = MatSolve_SeqSBAIJ_N_inplace;
716: C->ops->solvetranspose = MatSolve_SeqSBAIJ_N_inplace;
717: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_N_inplace;
718: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_N_inplace;
720: C->assembled = PETSC_TRUE;
721: C->preallocated = PETSC_TRUE;
723: PetscCall(PetscLogFlops(1.3333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering(Mat C, Mat A, const MatFactorInfo *info)
728: {
729: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
730: PetscInt i, j, mbs = a->mbs, *bi = b->i, *bj = b->j;
731: PetscInt *ai, *aj, k, k1, jmin, jmax, *jl, *il, vj, nexti, ili;
732: PetscInt bs = A->rmap->bs, bs2 = a->bs2;
733: MatScalar *ba = b->a, *aa, *ap, *dk, *uik;
734: MatScalar *u, *diag, *rtmp, *rtmp_ptr;
735: MatScalar *work;
736: PetscInt *pivots;
737: PetscBool allowzeropivot, zeropivotdetected;
739: PetscFunctionBegin;
740: PetscCall(PetscCalloc1(bs2 * mbs, &rtmp));
741: PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
742: il[0] = 0;
743: for (i = 0; i < mbs; i++) jl[i] = mbs;
745: PetscCall(PetscMalloc3(bs2, &dk, bs2, &uik, bs, &work));
746: PetscCall(PetscMalloc1(bs, &pivots));
747: allowzeropivot = PetscNot(A->erroriffailure);
749: ai = a->i;
750: aj = a->j;
751: aa = a->a;
753: /* for each row k */
754: for (k = 0; k < mbs; k++) {
755: /*initialize k-th row with elements nonzero in row k of A */
756: jmin = ai[k];
757: jmax = ai[k + 1];
758: ap = aa + jmin * bs2;
759: for (j = jmin; j < jmax; j++) {
760: vj = aj[j]; /* block col. index */
761: rtmp_ptr = rtmp + vj * bs2;
762: for (i = 0; i < bs2; i++) *rtmp_ptr++ = *ap++;
763: }
765: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
766: PetscCall(PetscArraycpy(dk, rtmp + k * bs2, bs2));
767: i = jl[k]; /* first row to be added to k_th row */
769: while (i < k) {
770: nexti = jl[i]; /* next row to be added to k_th row */
772: /* compute multiplier */
773: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
775: /* uik = -inv(Di)*U_bar(i,k) */
776: diag = ba + i * bs2;
777: u = ba + ili * bs2;
778: PetscCall(PetscArrayzero(uik, bs2));
779: PetscKernel_A_gets_A_minus_B_times_C(bs, uik, diag, u);
781: /* update D(k) += -U(i,k)^T * U_bar(i,k) */
782: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs, dk, uik, u);
783: PetscCall(PetscLogFlops(2.0 * bs * bs2));
785: /* update -U(i,k) */
786: PetscCall(PetscArraycpy(ba + ili * bs2, uik, bs2));
788: /* add multiple of row i to k-th row ... */
789: jmin = ili + 1;
790: jmax = bi[i + 1];
791: if (jmin < jmax) {
792: for (j = jmin; j < jmax; j++) {
793: /* rtmp += -U(i,k)^T * U_bar(i,j) */
794: rtmp_ptr = rtmp + bj[j] * bs2;
795: u = ba + j * bs2;
796: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs, rtmp_ptr, uik, u);
797: }
798: PetscCall(PetscLogFlops(2.0 * bs * bs2 * (jmax - jmin)));
800: /* ... add i to row list for next nonzero entry */
801: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
802: j = bj[jmin];
803: jl[i] = jl[j];
804: jl[j] = i; /* update jl */
805: }
806: i = nexti;
807: }
809: /* save nonzero entries in k-th row of U ... */
811: /* invert diagonal block */
812: diag = ba + k * bs2;
813: PetscCall(PetscArraycpy(diag, dk, bs2));
815: PetscCall(PetscKernel_A_gets_inverse_A(bs, diag, pivots, work, allowzeropivot, &zeropivotdetected));
816: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
818: jmin = bi[k];
819: jmax = bi[k + 1];
820: if (jmin < jmax) {
821: for (j = jmin; j < jmax; j++) {
822: vj = bj[j]; /* block col. index of U */
823: u = ba + j * bs2;
824: rtmp_ptr = rtmp + vj * bs2;
825: for (k1 = 0; k1 < bs2; k1++) {
826: *u++ = *rtmp_ptr;
827: *rtmp_ptr++ = 0.0;
828: }
829: }
831: /* ... add k to row list for first nonzero entry in k-th row */
832: il[k] = jmin;
833: i = bj[jmin];
834: jl[k] = jl[i];
835: jl[i] = k;
836: }
837: }
839: PetscCall(PetscFree(rtmp));
840: PetscCall(PetscFree2(il, jl));
841: PetscCall(PetscFree3(dk, uik, work));
842: PetscCall(PetscFree(pivots));
844: C->ops->solve = MatSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
845: C->ops->solvetranspose = MatSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
846: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
847: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
848: C->assembled = PETSC_TRUE;
849: C->preallocated = PETSC_TRUE;
851: PetscCall(PetscLogFlops(1.3333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: /*
856: Numeric U^T*D*U factorization for SBAIJ format. Modified from SNF of YSMP.
857: Version for blocks 2 by 2.
858: */
859: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat C, Mat A, const MatFactorInfo *info)
860: {
861: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
862: IS perm = b->row;
863: const PetscInt *ai, *aj, *perm_ptr;
864: PetscInt i, j, mbs = a->mbs, *bi = b->i, *bj = b->j;
865: PetscInt *a2anew, k, k1, jmin, jmax, *jl, *il, vj, nexti, ili;
866: MatScalar *ba = b->a, *aa, *ap;
867: MatScalar *u, *diag, *rtmp, *rtmp_ptr, dk[4], uik[4];
868: PetscReal shift = info->shiftamount;
869: PetscBool allowzeropivot, zeropivotdetected;
871: PetscFunctionBegin;
872: allowzeropivot = PetscNot(A->erroriffailure);
874: /* initialization */
875: /* il and jl record the first nonzero element in each row of the accessing
876: window U(0:k, k:mbs-1).
877: jl: list of rows to be added to uneliminated rows
878: i>= k: jl(i) is the first row to be added to row i
879: i< k: jl(i) is the row following row i in some list of rows
880: jl(i) = mbs indicates the end of a list
881: il(i): points to the first nonzero element in columns k,...,mbs-1 of
882: row i of U */
883: PetscCall(PetscCalloc1(4 * mbs, &rtmp));
884: PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
885: il[0] = 0;
886: for (i = 0; i < mbs; i++) jl[i] = mbs;
888: PetscCall(ISGetIndices(perm, &perm_ptr));
890: /* check permutation */
891: if (!a->permute) {
892: ai = a->i;
893: aj = a->j;
894: aa = a->a;
895: } else {
896: ai = a->inew;
897: aj = a->jnew;
898: PetscCall(PetscMalloc1(4 * ai[mbs], &aa));
899: PetscCall(PetscArraycpy(aa, a->a, 4 * ai[mbs]));
900: PetscCall(PetscMalloc1(ai[mbs], &a2anew));
901: PetscCall(PetscArraycpy(a2anew, a->a2anew, ai[mbs]));
903: for (i = 0; i < mbs; i++) {
904: jmin = ai[i];
905: jmax = ai[i + 1];
906: for (j = jmin; j < jmax; j++) {
907: while (a2anew[j] != j) {
908: k = a2anew[j];
909: a2anew[j] = a2anew[k];
910: a2anew[k] = k;
911: for (k1 = 0; k1 < 4; k1++) {
912: dk[k1] = aa[k * 4 + k1];
913: aa[k * 4 + k1] = aa[j * 4 + k1];
914: aa[j * 4 + k1] = dk[k1];
915: }
916: }
917: /* transform column-oriented blocks that lie in the lower triangle to row-oriented blocks */
918: if (i > aj[j]) {
919: ap = aa + j * 4; /* ptr to the beginning of the block */
920: dk[1] = ap[1]; /* swap ap[1] and ap[2] */
921: ap[1] = ap[2];
922: ap[2] = dk[1];
923: }
924: }
925: }
926: PetscCall(PetscFree(a2anew));
927: }
929: /* for each row k */
930: for (k = 0; k < mbs; k++) {
931: /*initialize k-th row with elements nonzero in row perm(k) of A */
932: jmin = ai[perm_ptr[k]];
933: jmax = ai[perm_ptr[k] + 1];
934: ap = aa + jmin * 4;
935: for (j = jmin; j < jmax; j++) {
936: vj = perm_ptr[aj[j]]; /* block col. index */
937: rtmp_ptr = rtmp + vj * 4;
938: for (i = 0; i < 4; i++) *rtmp_ptr++ = *ap++;
939: }
941: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
942: PetscCall(PetscArraycpy(dk, rtmp + k * 4, 4));
943: i = jl[k]; /* first row to be added to k_th row */
945: while (i < k) {
946: nexti = jl[i]; /* next row to be added to k_th row */
948: /* compute multiplier */
949: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
951: /* uik = -inv(Di)*U_bar(i,k): - ba[ili]*ba[i] */
952: diag = ba + i * 4;
953: u = ba + ili * 4;
954: uik[0] = -(diag[0] * u[0] + diag[2] * u[1]);
955: uik[1] = -(diag[1] * u[0] + diag[3] * u[1]);
956: uik[2] = -(diag[0] * u[2] + diag[2] * u[3]);
957: uik[3] = -(diag[1] * u[2] + diag[3] * u[3]);
959: /* update D(k) += -U(i,k)^T * U_bar(i,k): dk += uik*ba[ili] */
960: dk[0] += uik[0] * u[0] + uik[1] * u[1];
961: dk[1] += uik[2] * u[0] + uik[3] * u[1];
962: dk[2] += uik[0] * u[2] + uik[1] * u[3];
963: dk[3] += uik[2] * u[2] + uik[3] * u[3];
965: PetscCall(PetscLogFlops(16.0 * 2.0));
967: /* update -U(i,k): ba[ili] = uik */
968: PetscCall(PetscArraycpy(ba + ili * 4, uik, 4));
970: /* add multiple of row i to k-th row ... */
971: jmin = ili + 1;
972: jmax = bi[i + 1];
973: if (jmin < jmax) {
974: for (j = jmin; j < jmax; j++) {
975: /* rtmp += -U(i,k)^T * U_bar(i,j): rtmp[bj[j]] += uik*ba[j]; */
976: rtmp_ptr = rtmp + bj[j] * 4;
977: u = ba + j * 4;
978: rtmp_ptr[0] += uik[0] * u[0] + uik[1] * u[1];
979: rtmp_ptr[1] += uik[2] * u[0] + uik[3] * u[1];
980: rtmp_ptr[2] += uik[0] * u[2] + uik[1] * u[3];
981: rtmp_ptr[3] += uik[2] * u[2] + uik[3] * u[3];
982: }
983: PetscCall(PetscLogFlops(16.0 * (jmax - jmin)));
985: /* ... add i to row list for next nonzero entry */
986: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
987: j = bj[jmin];
988: jl[i] = jl[j];
989: jl[j] = i; /* update jl */
990: }
991: i = nexti;
992: }
994: /* save nonzero entries in k-th row of U ... */
996: /* invert diagonal block */
997: diag = ba + k * 4;
998: PetscCall(PetscArraycpy(diag, dk, 4));
999: PetscCall(PetscKernel_A_gets_inverse_A_2(diag, shift, allowzeropivot, &zeropivotdetected));
1000: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1002: jmin = bi[k];
1003: jmax = bi[k + 1];
1004: if (jmin < jmax) {
1005: for (j = jmin; j < jmax; j++) {
1006: vj = bj[j]; /* block col. index of U */
1007: u = ba + j * 4;
1008: rtmp_ptr = rtmp + vj * 4;
1009: for (k1 = 0; k1 < 4; k1++) {
1010: *u++ = *rtmp_ptr;
1011: *rtmp_ptr++ = 0.0;
1012: }
1013: }
1015: /* ... add k to row list for first nonzero entry in k-th row */
1016: il[k] = jmin;
1017: i = bj[jmin];
1018: jl[k] = jl[i];
1019: jl[i] = k;
1020: }
1021: }
1023: PetscCall(PetscFree(rtmp));
1024: PetscCall(PetscFree2(il, jl));
1025: if (a->permute) PetscCall(PetscFree(aa));
1026: PetscCall(ISRestoreIndices(perm, &perm_ptr));
1028: C->ops->solve = MatSolve_SeqSBAIJ_2_inplace;
1029: C->ops->solvetranspose = MatSolve_SeqSBAIJ_2_inplace;
1030: C->assembled = PETSC_TRUE;
1031: C->preallocated = PETSC_TRUE;
1033: PetscCall(PetscLogFlops(1.3333 * 8 * b->mbs)); /* from inverting diagonal blocks */
1034: PetscFunctionReturn(PETSC_SUCCESS);
1035: }
1037: /*
1038: Version for when blocks are 2 by 2 Using natural ordering
1039: */
1040: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering(Mat C, Mat A, const MatFactorInfo *info)
1041: {
1042: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
1043: PetscInt i, j, mbs = a->mbs, *bi = b->i, *bj = b->j;
1044: PetscInt *ai, *aj, k, k1, jmin, jmax, *jl, *il, vj, nexti, ili;
1045: MatScalar *ba = b->a, *aa, *ap, dk[8], uik[8];
1046: MatScalar *u, *diag, *rtmp, *rtmp_ptr;
1047: PetscReal shift = info->shiftamount;
1048: PetscBool allowzeropivot, zeropivotdetected;
1050: PetscFunctionBegin;
1051: allowzeropivot = PetscNot(A->erroriffailure);
1053: /* initialization */
1054: /* il and jl record the first nonzero element in each row of the accessing
1055: window U(0:k, k:mbs-1).
1056: jl: list of rows to be added to uneliminated rows
1057: i>= k: jl(i) is the first row to be added to row i
1058: i< k: jl(i) is the row following row i in some list of rows
1059: jl(i) = mbs indicates the end of a list
1060: il(i): points to the first nonzero element in columns k,...,mbs-1 of
1061: row i of U */
1062: PetscCall(PetscCalloc1(4 * mbs, &rtmp));
1063: PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
1064: il[0] = 0;
1065: for (i = 0; i < mbs; i++) jl[i] = mbs;
1067: ai = a->i;
1068: aj = a->j;
1069: aa = a->a;
1071: /* for each row k */
1072: for (k = 0; k < mbs; k++) {
1073: /*initialize k-th row with elements nonzero in row k of A */
1074: jmin = ai[k];
1075: jmax = ai[k + 1];
1076: ap = aa + jmin * 4;
1077: for (j = jmin; j < jmax; j++) {
1078: vj = aj[j]; /* block col. index */
1079: rtmp_ptr = rtmp + vj * 4;
1080: for (i = 0; i < 4; i++) *rtmp_ptr++ = *ap++;
1081: }
1083: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
1084: PetscCall(PetscArraycpy(dk, rtmp + k * 4, 4));
1085: i = jl[k]; /* first row to be added to k_th row */
1087: while (i < k) {
1088: nexti = jl[i]; /* next row to be added to k_th row */
1090: /* compute multiplier */
1091: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1093: /* uik = -inv(Di)*U_bar(i,k): - ba[ili]*ba[i] */
1094: diag = ba + i * 4;
1095: u = ba + ili * 4;
1096: uik[0] = -(diag[0] * u[0] + diag[2] * u[1]);
1097: uik[1] = -(diag[1] * u[0] + diag[3] * u[1]);
1098: uik[2] = -(diag[0] * u[2] + diag[2] * u[3]);
1099: uik[3] = -(diag[1] * u[2] + diag[3] * u[3]);
1101: /* update D(k) += -U(i,k)^T * U_bar(i,k): dk += uik*ba[ili] */
1102: dk[0] += uik[0] * u[0] + uik[1] * u[1];
1103: dk[1] += uik[2] * u[0] + uik[3] * u[1];
1104: dk[2] += uik[0] * u[2] + uik[1] * u[3];
1105: dk[3] += uik[2] * u[2] + uik[3] * u[3];
1107: PetscCall(PetscLogFlops(16.0 * 2.0));
1109: /* update -U(i,k): ba[ili] = uik */
1110: PetscCall(PetscArraycpy(ba + ili * 4, uik, 4));
1112: /* add multiple of row i to k-th row ... */
1113: jmin = ili + 1;
1114: jmax = bi[i + 1];
1115: if (jmin < jmax) {
1116: for (j = jmin; j < jmax; j++) {
1117: /* rtmp += -U(i,k)^T * U_bar(i,j): rtmp[bj[j]] += uik*ba[j]; */
1118: rtmp_ptr = rtmp + bj[j] * 4;
1119: u = ba + j * 4;
1120: rtmp_ptr[0] += uik[0] * u[0] + uik[1] * u[1];
1121: rtmp_ptr[1] += uik[2] * u[0] + uik[3] * u[1];
1122: rtmp_ptr[2] += uik[0] * u[2] + uik[1] * u[3];
1123: rtmp_ptr[3] += uik[2] * u[2] + uik[3] * u[3];
1124: }
1125: PetscCall(PetscLogFlops(16.0 * (jmax - jmin)));
1127: /* ... add i to row list for next nonzero entry */
1128: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
1129: j = bj[jmin];
1130: jl[i] = jl[j];
1131: jl[j] = i; /* update jl */
1132: }
1133: i = nexti;
1134: }
1136: /* save nonzero entries in k-th row of U ... */
1138: /* invert diagonal block */
1139: diag = ba + k * 4;
1140: PetscCall(PetscArraycpy(diag, dk, 4));
1141: PetscCall(PetscKernel_A_gets_inverse_A_2(diag, shift, allowzeropivot, &zeropivotdetected));
1142: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1144: jmin = bi[k];
1145: jmax = bi[k + 1];
1146: if (jmin < jmax) {
1147: for (j = jmin; j < jmax; j++) {
1148: vj = bj[j]; /* block col. index of U */
1149: u = ba + j * 4;
1150: rtmp_ptr = rtmp + vj * 4;
1151: for (k1 = 0; k1 < 4; k1++) {
1152: *u++ = *rtmp_ptr;
1153: *rtmp_ptr++ = 0.0;
1154: }
1155: }
1157: /* ... add k to row list for first nonzero entry in k-th row */
1158: il[k] = jmin;
1159: i = bj[jmin];
1160: jl[k] = jl[i];
1161: jl[i] = k;
1162: }
1163: }
1165: PetscCall(PetscFree(rtmp));
1166: PetscCall(PetscFree2(il, jl));
1168: C->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1169: C->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1170: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1171: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1172: C->assembled = PETSC_TRUE;
1173: C->preallocated = PETSC_TRUE;
1175: PetscCall(PetscLogFlops(1.3333 * 8 * b->mbs)); /* from inverting diagonal blocks */
1176: PetscFunctionReturn(PETSC_SUCCESS);
1177: }
1179: /*
1180: Numeric U^T*D*U factorization for SBAIJ format.
1181: Version for blocks are 1 by 1.
1182: */
1183: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace(Mat C, Mat A, const MatFactorInfo *info)
1184: {
1185: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
1186: IS ip = b->row;
1187: const PetscInt *ai, *aj, *rip;
1188: PetscInt *a2anew, i, j, mbs = a->mbs, *bi = b->i, *bj = b->j, *bcol;
1189: PetscInt k, jmin, jmax, *jl, *il, col, nexti, ili, nz;
1190: MatScalar *rtmp, *ba = b->a, *bval, *aa, dk, uikdi;
1191: PetscReal rs;
1192: FactorShiftCtx sctx;
1194: PetscFunctionBegin;
1195: /* MatPivotSetUp(): initialize shift context sctx */
1196: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
1198: PetscCall(ISGetIndices(ip, &rip));
1199: if (!a->permute) {
1200: ai = a->i;
1201: aj = a->j;
1202: aa = a->a;
1203: } else {
1204: ai = a->inew;
1205: aj = a->jnew;
1206: nz = ai[mbs];
1207: PetscCall(PetscMalloc1(nz, &aa));
1208: a2anew = a->a2anew;
1209: bval = a->a;
1210: for (j = 0; j < nz; j++) aa[a2anew[j]] = *(bval++);
1211: }
1213: /* initialization */
1214: /* il and jl record the first nonzero element in each row of the accessing
1215: window U(0:k, k:mbs-1).
1216: jl: list of rows to be added to uneliminated rows
1217: i>= k: jl(i) is the first row to be added to row i
1218: i< k: jl(i) is the row following row i in some list of rows
1219: jl(i) = mbs indicates the end of a list
1220: il(i): points to the first nonzero element in columns k,...,mbs-1 of
1221: row i of U */
1222: PetscCall(PetscMalloc3(mbs, &rtmp, mbs, &il, mbs, &jl));
1224: do {
1225: sctx.newshift = PETSC_FALSE;
1226: il[0] = 0;
1227: for (i = 0; i < mbs; i++) {
1228: rtmp[i] = 0.0;
1229: jl[i] = mbs;
1230: }
1232: for (k = 0; k < mbs; k++) {
1233: /*initialize k-th row by the perm[k]-th row of A */
1234: jmin = ai[rip[k]];
1235: jmax = ai[rip[k] + 1];
1236: bval = ba + bi[k];
1237: for (j = jmin; j < jmax; j++) {
1238: col = rip[aj[j]];
1239: rtmp[col] = aa[j];
1240: *bval++ = 0.0; /* for in-place factorization */
1241: }
1243: /* shift the diagonal of the matrix */
1244: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
1246: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1247: dk = rtmp[k];
1248: i = jl[k]; /* first row to be added to k_th row */
1250: while (i < k) {
1251: nexti = jl[i]; /* next row to be added to k_th row */
1253: /* compute multiplier, update diag(k) and U(i,k) */
1254: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1255: uikdi = -ba[ili] * ba[bi[i]]; /* diagonal(k) */
1256: dk += uikdi * ba[ili];
1257: ba[ili] = uikdi; /* -U(i,k) */
1259: /* add multiple of row i to k-th row */
1260: jmin = ili + 1;
1261: jmax = bi[i + 1];
1262: if (jmin < jmax) {
1263: for (j = jmin; j < jmax; j++) rtmp[bj[j]] += uikdi * ba[j];
1264: PetscCall(PetscLogFlops(2.0 * (jmax - jmin)));
1266: /* update il and jl for row i */
1267: il[i] = jmin;
1268: j = bj[jmin];
1269: jl[i] = jl[j];
1270: jl[j] = i;
1271: }
1272: i = nexti;
1273: }
1275: /* shift the diagonals when zero pivot is detected */
1276: /* compute rs=sum of abs(off-diagonal) */
1277: rs = 0.0;
1278: jmin = bi[k] + 1;
1279: nz = bi[k + 1] - jmin;
1280: if (nz) {
1281: bcol = bj + jmin;
1282: while (nz--) {
1283: rs += PetscAbsScalar(rtmp[*bcol]);
1284: bcol++;
1285: }
1286: }
1288: sctx.rs = rs;
1289: sctx.pv = dk;
1290: PetscCall(MatPivotCheck(C, A, info, &sctx, k));
1291: if (sctx.newshift) break; /* sctx.shift_amount is updated */
1292: dk = sctx.pv;
1294: /* copy data into U(k,:) */
1295: ba[bi[k]] = 1.0 / dk; /* U(k,k) */
1296: jmin = bi[k] + 1;
1297: jmax = bi[k + 1];
1298: if (jmin < jmax) {
1299: for (j = jmin; j < jmax; j++) {
1300: col = bj[j];
1301: ba[j] = rtmp[col];
1302: rtmp[col] = 0.0;
1303: }
1304: /* add the k-th row into il and jl */
1305: il[k] = jmin;
1306: i = bj[jmin];
1307: jl[k] = jl[i];
1308: jl[i] = k;
1309: }
1310: }
1311: } while (sctx.newshift);
1312: PetscCall(PetscFree3(rtmp, il, jl));
1313: if (a->permute) PetscCall(PetscFree(aa));
1315: PetscCall(ISRestoreIndices(ip, &rip));
1317: C->ops->solve = MatSolve_SeqSBAIJ_1_inplace;
1318: C->ops->solves = MatSolves_SeqSBAIJ_1_inplace;
1319: C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
1320: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_inplace;
1321: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_inplace;
1322: C->assembled = PETSC_TRUE;
1323: C->preallocated = PETSC_TRUE;
1325: PetscCall(PetscLogFlops(C->rmap->N));
1326: if (sctx.nshift) {
1327: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1328: PetscCall(PetscInfo(A, "number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1329: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1330: PetscCall(PetscInfo(A, "number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1331: }
1332: }
1333: PetscFunctionReturn(PETSC_SUCCESS);
1334: }
1336: /*
1337: Version for when blocks are 1 by 1 Using natural ordering under new datastructure
1338: Modified from MatCholeskyFactorNumeric_SeqAIJ()
1339: */
1340: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info)
1341: {
1342: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1343: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data;
1344: PetscInt i, j, mbs = A->rmap->n, *bi = b->i, *bj = b->j, *bdiag = b->diag, *bjtmp;
1345: PetscInt *ai = a->i, *aj = a->j, *ajtmp;
1346: PetscInt k, jmin, jmax, *c2r, *il, col, nexti, ili, nz;
1347: MatScalar *rtmp, *ba = b->a, *bval, *aa = a->a, dk, uikdi;
1348: FactorShiftCtx sctx;
1349: PetscReal rs;
1350: MatScalar d, *v;
1351: const PetscInt *adiag;
1353: PetscFunctionBegin;
1354: PetscCall(MatGetDiagonalMarkers_SeqSBAIJ(A, &adiag, NULL));
1355: PetscCall(PetscMalloc3(mbs, &rtmp, mbs, &il, mbs, &c2r));
1357: /* MatPivotSetUp(): initialize shift context sctx */
1358: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
1360: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
1361: sctx.shift_top = info->zeropivot;
1363: PetscCall(PetscArrayzero(rtmp, mbs));
1365: for (i = 0; i < mbs; i++) {
1366: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
1367: d = aa[adiag[i]];
1368: rtmp[i] += -PetscRealPart(d); /* diagonal entry */
1369: ajtmp = aj + ai[i] + 1; /* exclude diagonal */
1370: v = aa + ai[i] + 1;
1371: nz = ai[i + 1] - ai[i] - 1;
1372: for (j = 0; j < nz; j++) {
1373: rtmp[i] += PetscAbsScalar(v[j]);
1374: rtmp[ajtmp[j]] += PetscAbsScalar(v[j]);
1375: }
1376: if (PetscRealPart(rtmp[i]) > sctx.shift_top) sctx.shift_top = PetscRealPart(rtmp[i]);
1377: }
1378: sctx.shift_top *= 1.1;
1379: sctx.nshift_max = 5;
1380: sctx.shift_lo = 0.;
1381: sctx.shift_hi = 1.;
1382: }
1384: /* allocate working arrays
1385: c2r: linked list, keep track of pivot rows for a given column. c2r[col]: head of the list for a given col
1386: 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
1387: */
1388: do {
1389: sctx.newshift = PETSC_FALSE;
1391: for (i = 0; i < mbs; i++) c2r[i] = mbs;
1392: if (mbs) il[0] = 0;
1394: for (k = 0; k < mbs; k++) {
1395: /* zero rtmp */
1396: nz = bi[k + 1] - bi[k];
1397: bjtmp = bj + bi[k];
1398: for (j = 0; j < nz; j++) rtmp[bjtmp[j]] = 0.0;
1400: /* load in initial unfactored row */
1401: bval = ba + bi[k];
1402: jmin = ai[k];
1403: jmax = ai[k + 1];
1404: for (j = jmin; j < jmax; j++) {
1405: col = aj[j];
1406: rtmp[col] = aa[j];
1407: *bval++ = 0.0; /* for in-place factorization */
1408: }
1409: /* shift the diagonal of the matrix: ZeropivotApply() */
1410: rtmp[k] += sctx.shift_amount; /* shift the diagonal of the matrix */
1412: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1413: dk = rtmp[k];
1414: i = c2r[k]; /* first row to be added to k_th row */
1416: while (i < k) {
1417: nexti = c2r[i]; /* next row to be added to k_th row */
1419: /* compute multiplier, update diag(k) and U(i,k) */
1420: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1421: uikdi = -ba[ili] * ba[bdiag[i]]; /* diagonal(k) */
1422: dk += uikdi * ba[ili]; /* update diag[k] */
1423: ba[ili] = uikdi; /* -U(i,k) */
1425: /* add multiple of row i to k-th row */
1426: jmin = ili + 1;
1427: jmax = bi[i + 1];
1428: if (jmin < jmax) {
1429: for (j = jmin; j < jmax; j++) rtmp[bj[j]] += uikdi * ba[j];
1430: /* update il and c2r for row i */
1431: il[i] = jmin;
1432: j = bj[jmin];
1433: c2r[i] = c2r[j];
1434: c2r[j] = i;
1435: }
1436: i = nexti;
1437: }
1439: /* copy data into U(k,:) */
1440: rs = 0.0;
1441: jmin = bi[k];
1442: jmax = bi[k + 1] - 1;
1443: if (jmin < jmax) {
1444: for (j = jmin; j < jmax; j++) {
1445: col = bj[j];
1446: ba[j] = rtmp[col];
1447: rs += PetscAbsScalar(ba[j]);
1448: }
1449: /* add the k-th row into il and c2r */
1450: il[k] = jmin;
1451: i = bj[jmin];
1452: c2r[k] = c2r[i];
1453: c2r[i] = k;
1454: }
1456: sctx.rs = rs;
1457: sctx.pv = dk;
1458: PetscCall(MatPivotCheck(B, A, info, &sctx, k));
1459: if (sctx.newshift) break;
1460: dk = sctx.pv;
1462: ba[bdiag[k]] = 1.0 / dk; /* U(k,k) */
1463: }
1464: } while (sctx.newshift);
1466: PetscCall(PetscFree3(rtmp, il, c2r));
1468: B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1469: B->ops->solves = MatSolves_SeqSBAIJ_1;
1470: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1471: B->ops->matsolve = MatMatSolve_SeqSBAIJ_1_NaturalOrdering;
1472: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
1473: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;
1475: B->assembled = PETSC_TRUE;
1476: B->preallocated = PETSC_TRUE;
1478: PetscCall(PetscLogFlops(B->rmap->n));
1480: /* MatPivotView() */
1481: if (sctx.nshift) {
1482: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1483: 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));
1484: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1485: PetscCall(PetscInfo(A, "number of shift_nz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1486: } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
1487: PetscCall(PetscInfo(A, "number of shift_inblocks applied %" PetscInt_FMT ", each shift_amount %g\n", sctx.nshift, (double)info->shiftamount));
1488: }
1489: }
1490: PetscFunctionReturn(PETSC_SUCCESS);
1491: }
1493: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace(Mat C, Mat A, const MatFactorInfo *info)
1494: {
1495: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
1496: PetscInt i, j, mbs = a->mbs;
1497: PetscInt *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
1498: PetscInt k, jmin, *jl, *il, nexti, ili, *acol, *bcol, nz;
1499: MatScalar *rtmp, *ba = b->a, *aa = a->a, dk, uikdi, *aval, *bval;
1500: PetscReal rs;
1501: FactorShiftCtx sctx;
1503: PetscFunctionBegin;
1504: /* MatPivotSetUp(): initialize shift context sctx */
1505: PetscCall(PetscMemzero(&sctx, sizeof(FactorShiftCtx)));
1507: /* initialization */
1508: /* il and jl record the first nonzero element in each row of the accessing
1509: window U(0:k, k:mbs-1).
1510: jl: list of rows to be added to uneliminated rows
1511: i>= k: jl(i) is the first row to be added to row i
1512: i< k: jl(i) is the row following row i in some list of rows
1513: jl(i) = mbs indicates the end of a list
1514: il(i): points to the first nonzero element in U(i,k:mbs-1)
1515: */
1516: PetscCall(PetscMalloc1(mbs, &rtmp));
1517: PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
1519: do {
1520: sctx.newshift = PETSC_FALSE;
1521: il[0] = 0;
1522: for (i = 0; i < mbs; i++) {
1523: rtmp[i] = 0.0;
1524: jl[i] = mbs;
1525: }
1527: for (k = 0; k < mbs; k++) {
1528: /*initialize k-th row with elements nonzero in row perm(k) of A */
1529: nz = ai[k + 1] - ai[k];
1530: acol = aj + ai[k];
1531: aval = aa + ai[k];
1532: bval = ba + bi[k];
1533: while (nz--) {
1534: rtmp[*acol++] = *aval++;
1535: *bval++ = 0.0; /* for in-place factorization */
1536: }
1538: /* shift the diagonal of the matrix */
1539: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
1541: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1542: dk = rtmp[k];
1543: i = jl[k]; /* first row to be added to k_th row */
1545: while (i < k) {
1546: nexti = jl[i]; /* next row to be added to k_th row */
1547: /* compute multiplier, update D(k) and U(i,k) */
1548: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1549: uikdi = -ba[ili] * ba[bi[i]];
1550: dk += uikdi * ba[ili];
1551: ba[ili] = uikdi; /* -U(i,k) */
1553: /* add multiple of row i to k-th row ... */
1554: jmin = ili + 1;
1555: nz = bi[i + 1] - jmin;
1556: if (nz > 0) {
1557: bcol = bj + jmin;
1558: bval = ba + jmin;
1559: PetscCall(PetscLogFlops(2.0 * nz));
1560: while (nz--) rtmp[*bcol++] += uikdi * (*bval++);
1562: /* update il and jl for i-th row */
1563: il[i] = jmin;
1564: j = bj[jmin];
1565: jl[i] = jl[j];
1566: jl[j] = i;
1567: }
1568: i = nexti;
1569: }
1571: /* shift the diagonals when zero pivot is detected */
1572: /* compute rs=sum of abs(off-diagonal) */
1573: rs = 0.0;
1574: jmin = bi[k] + 1;
1575: nz = bi[k + 1] - jmin;
1576: if (nz) {
1577: bcol = bj + jmin;
1578: while (nz--) {
1579: rs += PetscAbsScalar(rtmp[*bcol]);
1580: bcol++;
1581: }
1582: }
1584: sctx.rs = rs;
1585: sctx.pv = dk;
1586: PetscCall(MatPivotCheck(C, A, info, &sctx, k));
1587: if (sctx.newshift) break; /* sctx.shift_amount is updated */
1588: dk = sctx.pv;
1590: /* copy data into U(k,:) */
1591: ba[bi[k]] = 1.0 / dk;
1592: jmin = bi[k] + 1;
1593: nz = bi[k + 1] - jmin;
1594: if (nz) {
1595: bcol = bj + jmin;
1596: bval = ba + jmin;
1597: while (nz--) {
1598: *bval++ = rtmp[*bcol];
1599: rtmp[*bcol++] = 0.0;
1600: }
1601: /* add k-th row into il and jl */
1602: il[k] = jmin;
1603: i = bj[jmin];
1604: jl[k] = jl[i];
1605: jl[i] = k;
1606: }
1607: } /* end of for (k = 0; k<mbs; k++) */
1608: } while (sctx.newshift);
1609: PetscCall(PetscFree(rtmp));
1610: PetscCall(PetscFree2(il, jl));
1612: C->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1613: C->ops->solves = MatSolves_SeqSBAIJ_1_inplace;
1614: C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1615: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1616: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1618: C->assembled = PETSC_TRUE;
1619: C->preallocated = PETSC_TRUE;
1621: PetscCall(PetscLogFlops(C->rmap->N));
1622: if (sctx.nshift) {
1623: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1624: PetscCall(PetscInfo(A, "number of shiftnz tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1625: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1626: PetscCall(PetscInfo(A, "number of shiftpd tries %" PetscInt_FMT ", shift_amount %g\n", sctx.nshift, (double)sctx.shift_amount));
1627: }
1628: }
1629: PetscFunctionReturn(PETSC_SUCCESS);
1630: }
1632: PetscErrorCode MatCholeskyFactor_SeqSBAIJ(Mat A, IS perm, const MatFactorInfo *info)
1633: {
1634: Mat C;
1636: PetscFunctionBegin;
1637: PetscCall(MatGetFactor(A, "petsc", MAT_FACTOR_CHOLESKY, &C));
1638: PetscCall(MatCholeskyFactorSymbolic(C, A, perm, info));
1639: PetscCall(MatCholeskyFactorNumeric(C, A, info));
1641: A->ops->solve = C->ops->solve;
1642: A->ops->solvetranspose = C->ops->solvetranspose;
1644: PetscCall(MatHeaderMerge(A, &C));
1645: PetscFunctionReturn(PETSC_SUCCESS);
1646: }