Actual source code: sbaijfact.c
petsc-3.11.4 2019-09-28
2: #include <../src/mat/impls/baij/seq/baij.h>
3: #include <../src/mat/impls/sbaij/seq/sbaij.h>
4: #include <petsc/private/kernels/blockinvert.h>
5: #include <petscis.h>
7: PetscErrorCode MatGetInertia_SeqSBAIJ(Mat F,PetscInt *nneg,PetscInt *nzero,PetscInt *npos)
8: {
9: Mat_SeqSBAIJ *fact=(Mat_SeqSBAIJ*)F->data;
10: MatScalar *dd=fact->a;
11: PetscInt mbs=fact->mbs,bs=F->rmap->bs,i,nneg_tmp,npos_tmp,*fi=fact->diag;
14: if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for bs: %D >1 yet",bs);
15: if (F->factorerrortype==MAT_FACTOR_NUMERIC_ZEROPIVOT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatFactor fails with numeric zeropivot");
17: nneg_tmp = 0; npos_tmp = 0;
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: if (nneg) *nneg = nneg_tmp;
24: if (npos) *npos = npos_tmp;
25: if (nzero) *nzero = mbs - nneg_tmp - npos_tmp;
26: return(0);
27: }
29: /*
30: Symbolic U^T*D*U factorization for SBAIJ format. Modified from SSF of YSMP.
31: Use Modified Sparse Row (MSR) storage for u and ju. See page 85, "Iterative Methods ..." by Saad.
32: */
33: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ_MSR(Mat F,Mat A,IS perm,const MatFactorInfo *info)
34: {
35: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data,*b;
37: const PetscInt *rip,*ai,*aj;
38: PetscInt i,mbs = a->mbs,*jutmp,bs = A->rmap->bs,bs2=a->bs2;
39: PetscInt m,reallocs = 0,prow;
40: PetscInt *jl,*q,jmin,jmax,juidx,nzk,qm,*iu,*ju,k,j,vj,umax,maxadd;
41: PetscReal f = info->fill;
42: PetscBool perm_identity;
45: /* check whether perm is the identity mapping */
46: ISIdentity(perm,&perm_identity);
47: ISGetIndices(perm,&rip);
49: if (perm_identity) { /* without permutation */
50: a->permute = PETSC_FALSE;
52: ai = a->i; aj = a->j;
53: } else { /* non-trivial permutation */
54: a->permute = PETSC_TRUE;
56: MatReorderingSeqSBAIJ(A,perm);
58: ai = a->inew; aj = a->jnew;
59: }
61: /* initialization */
62: PetscMalloc1(mbs+1,&iu);
63: umax = (PetscInt)(f*ai[mbs] + 1); umax += mbs + 1;
64: PetscMalloc1(umax,&ju);
65: iu[0] = mbs+1;
66: juidx = mbs + 1; /* index for ju */
67: /* jl linked list for pivot row -- linked list for col index */
68: PetscMalloc2(mbs,&jl,mbs,&q);
69: for (i=0; i<mbs; i++) {
70: jl[i] = mbs;
71: q[i] = 0;
72: }
74: /* for each row k */
75: for (k=0; k<mbs; k++) {
76: for (i=0; i<mbs; i++) q[i] = 0; /* to be removed! */
77: nzk = 0; /* num. of nz blocks in k-th block row with diagonal block excluded */
78: q[k] = mbs;
79: /* initialize nonzero structure of k-th row to row rip[k] of A */
80: jmin = ai[rip[k]] +1; /* exclude diag[k] */
81: jmax = ai[rip[k]+1];
82: for (j=jmin; j<jmax; j++) {
83: vj = rip[aj[j]]; /* col. value */
84: if (vj > k) {
85: qm = k;
86: do {
87: m = qm; qm = q[m];
88: } while (qm < vj);
89: if (qm == vj) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Duplicate entry in A\n");
90: nzk++;
91: q[m] = vj;
92: q[vj] = qm;
93: } /* if (vj > k) */
94: } /* for (j=jmin; j<jmax; j++) */
96: /* modify nonzero structure of k-th row by computing fill-in
97: for each row i to be merged in */
98: prow = k;
99: prow = jl[prow]; /* next pivot row (== mbs for symbolic factorization) */
101: while (prow < k) {
102: /* merge row prow into k-th row */
103: jmin = iu[prow] + 1; jmax = iu[prow+1];
104: qm = k;
105: for (j=jmin; j<jmax; j++) {
106: vj = ju[j];
107: do {
108: m = qm; qm = q[m];
109: } while (qm < vj);
110: if (qm != vj) {
111: nzk++; q[m] = vj; q[vj] = qm; qm = vj;
112: }
113: }
114: prow = jl[prow]; /* next pivot row */
115: }
117: /* add k to row list for first nonzero element in k-th row */
118: if (nzk > 0) {
119: i = q[k]; /* col value of first nonzero element in U(k, k+1:mbs-1) */
120: jl[k] = jl[i]; jl[i] = k;
121: }
122: iu[k+1] = iu[k] + nzk;
124: /* allocate more space to ju if needed */
125: if (iu[k+1] > umax) {
126: /* estimate how much additional space we will need */
127: /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
128: /* just double the memory each time */
129: maxadd = umax;
130: if (maxadd < nzk) maxadd = (mbs-k)*(nzk+1)/2;
131: umax += maxadd;
133: /* allocate a longer ju */
134: PetscMalloc1(umax,&jutmp);
135: PetscMemcpy(jutmp,ju,iu[k]*sizeof(PetscInt));
136: PetscFree(ju);
137: ju = jutmp;
138: reallocs++; /* count how many times we realloc */
139: }
141: /* save nonzero structure of k-th row in ju */
142: i=k;
143: while (nzk--) {
144: i = q[i];
145: ju[juidx++] = i;
146: }
147: }
149: #if defined(PETSC_USE_INFO)
150: if (ai[mbs] != 0) {
151: PetscReal af = ((PetscReal)iu[mbs])/((PetscReal)ai[mbs]);
152: PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);
153: PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
154: PetscInfo1(A,"PCFactorSetFill(pc,%g);\n",(double)af);
155: PetscInfo(A,"for best performance.\n");
156: } else {
157: PetscInfo(A,"Empty matrix.\n");
158: }
159: #endif
161: ISRestoreIndices(perm,&rip);
162: PetscFree2(jl,q);
164: /* put together the new matrix */
165: MatSeqSBAIJSetPreallocation(F,bs,MAT_SKIP_ALLOCATION,NULL);
167: /* PetscLogObjectParent((PetscObject)B,(PetscObject)iperm); */
168: b = (Mat_SeqSBAIJ*)(F)->data;
169: b->singlemalloc = PETSC_FALSE;
170: b->free_a = PETSC_TRUE;
171: b->free_ij = PETSC_TRUE;
173: PetscMalloc1((iu[mbs]+1)*bs2,&b->a);
174: b->j = ju;
175: b->i = iu;
176: b->diag = 0;
177: b->ilen = 0;
178: b->imax = 0;
179: b->row = perm;
181: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
183: PetscObjectReference((PetscObject)perm);
185: b->icol = perm;
186: PetscObjectReference((PetscObject)perm);
187: PetscMalloc1(bs*mbs+bs,&b->solve_work);
188: /* In b structure: Free imax, ilen, old a, old j.
189: Allocate idnew, solve_work, new a, new j */
190: PetscLogObjectMemory((PetscObject)F,(iu[mbs]-mbs)*(sizeof(PetscInt)+sizeof(MatScalar)));
191: b->maxnz = b->nz = iu[mbs];
193: (F)->info.factor_mallocs = reallocs;
194: (F)->info.fill_ratio_given = f;
195: if (ai[mbs] != 0) {
196: (F)->info.fill_ratio_needed = ((PetscReal)iu[mbs])/((PetscReal)ai[mbs]);
197: } else {
198: (F)->info.fill_ratio_needed = 0.0;
199: }
200: MatSeqSBAIJSetNumericFactorization_inplace(F,perm_identity);
201: return(0);
202: }
203: /*
204: Symbolic U^T*D*U factorization for SBAIJ format.
205: See MatICCFactorSymbolic_SeqAIJ() for description of its data structure.
206: */
207: #include <petscbt.h>
208: #include <../src/mat/utils/freespace.h>
209: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
210: {
211: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
212: Mat_SeqSBAIJ *b;
213: PetscErrorCode ierr;
214: PetscBool perm_identity,missing;
215: PetscReal fill = info->fill;
216: const PetscInt *rip,*ai=a->i,*aj=a->j;
217: PetscInt i,mbs=a->mbs,bs=A->rmap->bs,reallocs=0,prow;
218: PetscInt *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
219: PetscInt nlnk,*lnk,ncols,*cols,*uj,**ui_ptr,*uj_ptr,*udiag;
220: PetscFreeSpaceList free_space=NULL,current_space=NULL;
221: PetscBT lnkbt;
224: if (A->rmap->n != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %D columns %D",A->rmap->n,A->cmap->n);
225: MatMissingDiagonal(A,&missing,&i);
226: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);
227: if (bs > 1) {
228: MatCholeskyFactorSymbolic_SeqSBAIJ_inplace(fact,A,perm,info);
229: return(0);
230: }
232: /* check whether perm is the identity mapping */
233: ISIdentity(perm,&perm_identity);
234: if (!perm_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported for sbaij matrix. Use aij format");
235: a->permute = PETSC_FALSE;
236: ISGetIndices(perm,&rip);
238: /* initialization */
239: PetscMalloc1(mbs+1,&ui);
240: PetscMalloc1(mbs+1,&udiag);
241: ui[0] = 0;
243: /* jl: linked list for storing indices of the pivot rows
244: il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
245: PetscMalloc4(mbs,&ui_ptr,mbs,&il,mbs,&jl,mbs,&cols);
246: for (i=0; i<mbs; i++) {
247: jl[i] = mbs; il[i] = 0;
248: }
250: /* create and initialize a linked list for storing column indices of the active row k */
251: nlnk = mbs + 1;
252: PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);
254: /* initial FreeSpace size is fill*(ai[mbs]+1) */
255: PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,ai[mbs]+1),&free_space);
256: current_space = free_space;
258: for (k=0; k<mbs; k++) { /* for each active row k */
259: /* initialize lnk by the column indices of row rip[k] of A */
260: nzk = 0;
261: ncols = ai[k+1] - ai[k];
262: if (!ncols) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Empty row %D in matrix ",k);
263: for (j=0; j<ncols; j++) {
264: i = *(aj + ai[k] + j);
265: cols[j] = i;
266: }
267: PetscLLAdd(ncols,cols,mbs,nlnk,lnk,lnkbt);
268: nzk += nlnk;
270: /* update lnk by computing fill-in for each pivot row to be merged in */
271: prow = jl[k]; /* 1st pivot row */
273: while (prow < k) {
274: nextprow = jl[prow];
275: /* merge prow into k-th row */
276: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
277: jmax = ui[prow+1];
278: ncols = jmax-jmin;
279: uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
280: PetscLLAddSorted(ncols,uj_ptr,mbs,nlnk,lnk,lnkbt);
281: nzk += nlnk;
283: /* update il and jl for prow */
284: if (jmin < jmax) {
285: il[prow] = jmin;
286: j = *uj_ptr; jl[prow] = jl[j]; jl[j] = prow;
287: }
288: prow = nextprow;
289: }
291: /* if free space is not available, make more free space */
292: if (current_space->local_remaining<nzk) {
293: i = mbs - k + 1; /* num of unfactored rows */
294: i = PetscIntMultTruncate(i,PetscMin(nzk, i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
295: PetscFreeSpaceGet(i,¤t_space);
296: reallocs++;
297: }
299: /* copy data into free space, then initialize lnk */
300: PetscLLClean(mbs,mbs,nzk,lnk,current_space->array,lnkbt);
302: /* add the k-th row into il and jl */
303: if (nzk > 1) {
304: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
305: jl[k] = jl[i]; jl[i] = k;
306: il[k] = ui[k] + 1;
307: }
308: ui_ptr[k] = current_space->array;
310: current_space->array += nzk;
311: current_space->local_used += nzk;
312: current_space->local_remaining -= nzk;
314: ui[k+1] = ui[k] + nzk;
315: }
317: ISRestoreIndices(perm,&rip);
318: PetscFree4(ui_ptr,il,jl,cols);
320: /* destroy list of free space and other temporary array(s) */
321: PetscMalloc1(ui[mbs]+1,&uj);
322: PetscFreeSpaceContiguous_Cholesky(&free_space,uj,mbs,ui,udiag); /* store matrix factor */
323: PetscLLDestroy(lnk,lnkbt);
325: /* put together the new matrix in MATSEQSBAIJ format */
326: MatSeqSBAIJSetPreallocation(fact,bs,MAT_SKIP_ALLOCATION,NULL);
328: b = (Mat_SeqSBAIJ*)fact->data;
329: b->singlemalloc = PETSC_FALSE;
330: b->free_a = PETSC_TRUE;
331: b->free_ij = PETSC_TRUE;
333: PetscMalloc1(ui[mbs]+1,&b->a);
335: b->j = uj;
336: b->i = ui;
337: b->diag = udiag;
338: b->free_diag = PETSC_TRUE;
339: b->ilen = 0;
340: b->imax = 0;
341: b->row = perm;
342: b->icol = perm;
344: PetscObjectReference((PetscObject)perm);
345: PetscObjectReference((PetscObject)perm);
347: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
349: PetscMalloc1(mbs+1,&b->solve_work);
350: PetscLogObjectMemory((PetscObject)fact,ui[mbs]*(sizeof(PetscInt)+sizeof(MatScalar)));
352: b->maxnz = b->nz = ui[mbs];
354: fact->info.factor_mallocs = reallocs;
355: fact->info.fill_ratio_given = fill;
356: if (ai[mbs] != 0) {
357: fact->info.fill_ratio_needed = ((PetscReal)ui[mbs])/ai[mbs];
358: } else {
359: fact->info.fill_ratio_needed = 0.0;
360: }
361: #if defined(PETSC_USE_INFO)
362: if (ai[mbs] != 0) {
363: PetscReal af = fact->info.fill_ratio_needed;
364: PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
365: PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
366: PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
367: } else {
368: PetscInfo(A,"Empty matrix.\n");
369: }
370: #endif
371: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering;
372: return(0);
373: }
375: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ_inplace(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
376: {
377: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
378: Mat_SeqSBAIJ *b;
379: PetscErrorCode ierr;
380: PetscBool perm_identity,missing;
381: PetscReal fill = info->fill;
382: const PetscInt *rip,*ai,*aj;
383: PetscInt i,mbs=a->mbs,bs=A->rmap->bs,reallocs=0,prow,d;
384: PetscInt *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
385: PetscInt nlnk,*lnk,ncols,*cols,*uj,**ui_ptr,*uj_ptr;
386: PetscFreeSpaceList free_space=NULL,current_space=NULL;
387: PetscBT lnkbt;
390: MatMissingDiagonal(A,&missing,&d);
391: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
393: /*
394: This code originally uses Modified Sparse Row (MSR) storage
395: (see page 85, "Iterative Methods ..." by Saad) for the output matrix B - bad choise!
396: Then it is rewritten so the factor B takes seqsbaij format. However the associated
397: MatCholeskyFactorNumeric_() have not been modified for the cases of bs>1 or !perm_identity,
398: thus the original code in MSR format is still used for these cases.
399: The code below should replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR() whenever
400: MatCholeskyFactorNumeric_() is modified for using sbaij symbolic factor.
401: */
402: if (bs > 1) {
403: MatCholeskyFactorSymbolic_SeqSBAIJ_MSR(fact,A,perm,info);
404: return(0);
405: }
407: /* check whether perm is the identity mapping */
408: ISIdentity(perm,&perm_identity);
410: if (perm_identity) {
411: a->permute = PETSC_FALSE;
413: ai = a->i; aj = a->j;
414: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported for sbaij matrix. Use aij format");
415: ISGetIndices(perm,&rip);
417: /* initialization */
418: PetscMalloc1(mbs+1,&ui);
419: ui[0] = 0;
421: /* jl: linked list for storing indices of the pivot rows
422: il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
423: PetscMalloc4(mbs,&ui_ptr,mbs,&il,mbs,&jl,mbs,&cols);
424: for (i=0; i<mbs; i++) {
425: jl[i] = mbs; il[i] = 0;
426: }
428: /* create and initialize a linked list for storing column indices of the active row k */
429: nlnk = mbs + 1;
430: PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);
432: /* initial FreeSpace size is fill*(ai[mbs]+1) */
433: PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,ai[mbs]+1),&free_space);
434: current_space = free_space;
436: for (k=0; k<mbs; k++) { /* for each active row k */
437: /* initialize lnk by the column indices of row rip[k] of A */
438: nzk = 0;
439: ncols = ai[rip[k]+1] - ai[rip[k]];
440: for (j=0; j<ncols; j++) {
441: i = *(aj + ai[rip[k]] + j);
442: cols[j] = rip[i];
443: }
444: PetscLLAdd(ncols,cols,mbs,nlnk,lnk,lnkbt);
445: nzk += nlnk;
447: /* update lnk by computing fill-in for each pivot row to be merged in */
448: prow = jl[k]; /* 1st pivot row */
450: while (prow < k) {
451: nextprow = jl[prow];
452: /* merge prow into k-th row */
453: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
454: jmax = ui[prow+1];
455: ncols = jmax-jmin;
456: uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
457: PetscLLAddSorted(ncols,uj_ptr,mbs,nlnk,lnk,lnkbt);
458: nzk += nlnk;
460: /* update il and jl for prow */
461: if (jmin < jmax) {
462: il[prow] = jmin;
464: j = *uj_ptr; jl[prow] = jl[j]; jl[j] = prow;
465: }
466: prow = nextprow;
467: }
469: /* if free space is not available, make more free space */
470: if (current_space->local_remaining<nzk) {
471: i = mbs - k + 1; /* num of unfactored rows */
472: i = PetscMin(PetscIntMultTruncate(i,nzk), PetscIntMultTruncate(i,i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
473: PetscFreeSpaceGet(i,¤t_space);
474: reallocs++;
475: }
477: /* copy data into free space, then initialize lnk */
478: PetscLLClean(mbs,mbs,nzk,lnk,current_space->array,lnkbt);
480: /* add the k-th row into il and jl */
481: if (nzk-1 > 0) {
482: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
483: jl[k] = jl[i]; jl[i] = k;
484: il[k] = ui[k] + 1;
485: }
486: ui_ptr[k] = current_space->array;
488: current_space->array += nzk;
489: current_space->local_used += nzk;
490: current_space->local_remaining -= nzk;
492: ui[k+1] = ui[k] + nzk;
493: }
495: ISRestoreIndices(perm,&rip);
496: PetscFree4(ui_ptr,il,jl,cols);
498: /* destroy list of free space and other temporary array(s) */
499: PetscMalloc1(ui[mbs]+1,&uj);
500: PetscFreeSpaceContiguous(&free_space,uj);
501: PetscLLDestroy(lnk,lnkbt);
503: /* put together the new matrix in MATSEQSBAIJ format */
504: MatSeqSBAIJSetPreallocation(fact,bs,MAT_SKIP_ALLOCATION,NULL);
506: b = (Mat_SeqSBAIJ*)fact->data;
507: b->singlemalloc = PETSC_FALSE;
508: b->free_a = PETSC_TRUE;
509: b->free_ij = PETSC_TRUE;
511: PetscMalloc1(ui[mbs]+1,&b->a);
513: b->j = uj;
514: b->i = ui;
515: b->diag = 0;
516: b->ilen = 0;
517: b->imax = 0;
518: b->row = perm;
520: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
522: PetscObjectReference((PetscObject)perm);
523: b->icol = perm;
524: PetscObjectReference((PetscObject)perm);
525: PetscMalloc1(mbs+1,&b->solve_work);
526: PetscLogObjectMemory((PetscObject)fact,(ui[mbs]-mbs)*(sizeof(PetscInt)+sizeof(MatScalar)));
527: b->maxnz = b->nz = ui[mbs];
529: fact->info.factor_mallocs = reallocs;
530: fact->info.fill_ratio_given = fill;
531: if (ai[mbs] != 0) {
532: fact->info.fill_ratio_needed = ((PetscReal)ui[mbs])/ai[mbs];
533: } else {
534: fact->info.fill_ratio_needed = 0.0;
535: }
536: #if defined(PETSC_USE_INFO)
537: if (ai[mbs] != 0) {
538: PetscReal af = fact->info.fill_ratio_needed;
539: PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
540: PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
541: PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
542: } else {
543: PetscInfo(A,"Empty matrix.\n");
544: }
545: #endif
546: MatSeqSBAIJSetNumericFactorization_inplace(fact,perm_identity);
547: return(0);
548: }
550: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat C,Mat A,const MatFactorInfo *info)
551: {
552: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ*)C->data;
553: 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;
566: /* initialization */
567: PetscCalloc1(bs2*mbs,&rtmp);
568: PetscMalloc2(mbs,&il,mbs,&jl);
569: allowzeropivot = PetscNot(A->erroriffailure);
571: il[0] = 0;
572: for (i=0; i<mbs; i++) jl[i] = mbs;
573:
574: PetscMalloc3(bs2,&dk,bs2,&uik,bs,&work);
575: PetscMalloc1(bs,&pivots);
577: ISGetIndices(perm,&perm_ptr);
579: /* check permutation */
580: if (!a->permute) {
581: ai = a->i; aj = a->j; aa = a->a;
582: } else {
583: ai = a->inew; aj = a->jnew;
584: PetscMalloc1(bs2*ai[mbs],&aa);
585: PetscMemcpy(aa,a->a,bs2*ai[mbs]*sizeof(MatScalar));
586: PetscMalloc1(ai[mbs],&a2anew);
587: PetscMemcpy(a2anew,a->a2anew,(ai[mbs])*sizeof(PetscInt));
589: for (i=0; i<mbs; i++) {
590: jmin = ai[i]; jmax = ai[i+1];
591: for (j=jmin; j<jmax; j++) {
592: while (a2anew[j] != j) {
593: k = a2anew[j]; a2anew[j] = a2anew[k]; a2anew[k] = k;
594: for (k1=0; k1<bs2; k1++) {
595: dk[k1] = aa[k*bs2+k1];
596: aa[k*bs2+k1] = aa[j*bs2+k1];
597: aa[j*bs2+k1] = dk[k1];
598: }
599: }
600: /* transform columnoriented blocks that lie in the lower triangle to roworiented blocks */
601: if (i > aj[j]) {
602: /* printf("change orientation, row: %d, col: %d\n",i,aj[j]); */
603: ap = aa + j*bs2; /* ptr to the beginning of j-th block of aa */
604: for (k=0; k<bs2; k++) dk[k] = ap[k]; /* dk <- j-th block of aa */
605: for (k=0; k<bs; k++) { /* j-th block of aa <- dk^T */
606: for (k1=0; k1<bs; k1++) *ap++ = dk[k + bs*k1];
607: }
608: }
609: }
610: }
611: PetscFree(a2anew);
612: }
614: /* for each row k */
615: for (k = 0; k<mbs; k++) {
617: /*initialize k-th row with elements nonzero in row perm(k) of A */
618: jmin = ai[perm_ptr[k]]; jmax = ai[perm_ptr[k]+1];
620: ap = aa + jmin*bs2;
621: for (j = jmin; j < jmax; j++) {
622: vj = perm_ptr[aj[j]]; /* block col. index */
623: rtmp_ptr = rtmp + vj*bs2;
624: for (i=0; i<bs2; i++) *rtmp_ptr++ = *ap++;
625: }
627: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
628: PetscMemcpy(dk,rtmp+k*bs2,bs2*sizeof(MatScalar));
629: i = jl[k]; /* first row to be added to k_th row */
631: while (i < k) {
632: nexti = jl[i]; /* next row to be added to k_th row */
634: /* compute multiplier */
635: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
637: /* uik = -inv(Di)*U_bar(i,k) */
638: diag = ba + i*bs2;
639: u = ba + ili*bs2;
640: PetscMemzero(uik,bs2*sizeof(MatScalar));
641: PetscKernel_A_gets_A_minus_B_times_C(bs,uik,diag,u);
643: /* update D(k) += -U(i,k)^T * U_bar(i,k) */
644: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,dk,uik,u);
645: PetscLogFlops(4.0*bs*bs2);
647: /* update -U(i,k) */
648: PetscMemcpy(ba+ili*bs2,uik,bs2*sizeof(MatScalar));
650: /* add multiple of row i to k-th row ... */
651: jmin = ili + 1; jmax = bi[i+1];
652: if (jmin < jmax) {
653: for (j=jmin; j<jmax; j++) {
654: /* rtmp += -U(i,k)^T * U_bar(i,j) */
655: rtmp_ptr = rtmp + bj[j]*bs2;
656: u = ba + j*bs2;
657: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,rtmp_ptr,uik,u);
658: }
659: PetscLogFlops(2.0*bs*bs2*(jmax-jmin));
661: /* ... add i to row list for next nonzero entry */
662: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
663: j = bj[jmin];
664: jl[i] = jl[j]; jl[j] = i; /* update jl */
665: }
666: i = nexti;
667: }
669: /* save nonzero entries in k-th row of U ... */
671: /* invert diagonal block */
672: diag = ba+k*bs2;
673: PetscMemcpy(diag,dk,bs2*sizeof(MatScalar));
675: PetscKernel_A_gets_inverse_A(bs,diag,pivots,work,allowzeropivot,&zeropivotdetected);
676: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
678: jmin = bi[k]; jmax = bi[k+1];
679: if (jmin < jmax) {
680: for (j=jmin; j<jmax; j++) {
681: vj = bj[j]; /* block col. index of U */
682: u = ba + j*bs2;
683: rtmp_ptr = rtmp + vj*bs2;
684: for (k1=0; k1<bs2; k1++) {
685: *u++ = *rtmp_ptr;
686: *rtmp_ptr++ = 0.0;
687: }
688: }
690: /* ... add k to row list for first nonzero entry in k-th row */
691: il[k] = jmin;
692: i = bj[jmin];
693: jl[k] = jl[i]; jl[i] = k;
694: }
695: }
697: PetscFree(rtmp);
698: PetscFree2(il,jl);
699: PetscFree3(dk,uik,work);
700: PetscFree(pivots);
701: if (a->permute) {
702: PetscFree(aa);
703: }
705: ISRestoreIndices(perm,&perm_ptr);
707: C->ops->solve = MatSolve_SeqSBAIJ_N_inplace;
708: C->ops->solvetranspose = MatSolve_SeqSBAIJ_N_inplace;
709: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_N_inplace;
710: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_N_inplace;
712: C->assembled = PETSC_TRUE;
713: C->preallocated = PETSC_TRUE;
715: PetscLogFlops(1.3333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
716: return(0);
717: }
719: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info)
720: {
721: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ*)C->data;
723: PetscInt i,j,mbs=a->mbs,*bi=b->i,*bj=b->j;
724: PetscInt *ai,*aj,k,k1,jmin,jmax,*jl,*il,vj,nexti,ili;
725: PetscInt bs =A->rmap->bs,bs2 = a->bs2;
726: MatScalar *ba = b->a,*aa,*ap,*dk,*uik;
727: MatScalar *u,*diag,*rtmp,*rtmp_ptr;
728: MatScalar *work;
729: PetscInt *pivots;
730: PetscBool allowzeropivot,zeropivotdetected;
733: PetscCalloc1(bs2*mbs,&rtmp);
734: PetscMalloc2(mbs,&il,mbs,&jl);
735: il[0] = 0;
736: for (i=0; i<mbs; i++) jl[i] = mbs;
737:
738: PetscMalloc3(bs2,&dk,bs2,&uik,bs,&work);
739: PetscMalloc1(bs,&pivots);
740: allowzeropivot = PetscNot(A->erroriffailure);
742: ai = a->i; aj = a->j; aa = a->a;
744: /* for each row k */
745: for (k = 0; k<mbs; k++) {
747: /*initialize k-th row with elements nonzero in row k of A */
748: jmin = ai[k]; jmax = ai[k+1];
749: ap = aa + jmin*bs2;
750: for (j = jmin; j < jmax; j++) {
751: vj = aj[j]; /* block col. index */
752: rtmp_ptr = rtmp + vj*bs2;
753: for (i=0; i<bs2; i++) *rtmp_ptr++ = *ap++;
754: }
756: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
757: PetscMemcpy(dk,rtmp+k*bs2,bs2*sizeof(MatScalar));
758: i = jl[k]; /* first row to be added to k_th row */
760: while (i < k) {
761: nexti = jl[i]; /* next row to be added to k_th row */
763: /* compute multiplier */
764: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
766: /* uik = -inv(Di)*U_bar(i,k) */
767: diag = ba + i*bs2;
768: u = ba + ili*bs2;
769: PetscMemzero(uik,bs2*sizeof(MatScalar));
770: PetscKernel_A_gets_A_minus_B_times_C(bs,uik,diag,u);
772: /* update D(k) += -U(i,k)^T * U_bar(i,k) */
773: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,dk,uik,u);
774: PetscLogFlops(2.0*bs*bs2);
776: /* update -U(i,k) */
777: PetscMemcpy(ba+ili*bs2,uik,bs2*sizeof(MatScalar));
779: /* add multiple of row i to k-th row ... */
780: jmin = ili + 1; jmax = bi[i+1];
781: if (jmin < jmax) {
782: for (j=jmin; j<jmax; j++) {
783: /* rtmp += -U(i,k)^T * U_bar(i,j) */
784: rtmp_ptr = rtmp + bj[j]*bs2;
785: u = ba + j*bs2;
786: PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,rtmp_ptr,uik,u);
787: }
788: PetscLogFlops(2.0*bs*bs2*(jmax-jmin));
790: /* ... add i to row list for next nonzero entry */
791: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
792: j = bj[jmin];
793: jl[i] = jl[j]; jl[j] = i; /* update jl */
794: }
795: i = nexti;
796: }
798: /* save nonzero entries in k-th row of U ... */
800: /* invert diagonal block */
801: diag = ba+k*bs2;
802: PetscMemcpy(diag,dk,bs2*sizeof(MatScalar));
804: PetscKernel_A_gets_inverse_A(bs,diag,pivots,work,allowzeropivot,&zeropivotdetected);
805: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
807: jmin = bi[k]; jmax = bi[k+1];
808: if (jmin < jmax) {
809: for (j=jmin; j<jmax; j++) {
810: vj = bj[j]; /* block col. index of U */
811: u = ba + j*bs2;
812: rtmp_ptr = rtmp + vj*bs2;
813: for (k1=0; k1<bs2; k1++) {
814: *u++ = *rtmp_ptr;
815: *rtmp_ptr++ = 0.0;
816: }
817: }
819: /* ... add k to row list for first nonzero entry in k-th row */
820: il[k] = jmin;
821: i = bj[jmin];
822: jl[k] = jl[i]; jl[i] = k;
823: }
824: }
826: PetscFree(rtmp);
827: PetscFree2(il,jl);
828: PetscFree3(dk,uik,work);
829: PetscFree(pivots);
831: C->ops->solve = MatSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
832: C->ops->solvetranspose = MatSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
833: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
834: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
835: C->assembled = PETSC_TRUE;
836: C->preallocated = PETSC_TRUE;
838: PetscLogFlops(1.3333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
839: return(0);
840: }
842: /*
843: Numeric U^T*D*U factorization for SBAIJ format. Modified from SNF of YSMP.
844: Version for blocks 2 by 2.
845: */
846: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat C,Mat A,const MatFactorInfo *info)
847: {
848: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ*)C->data;
849: IS perm = b->row;
851: const PetscInt *ai,*aj,*perm_ptr;
852: PetscInt i,j,mbs=a->mbs,*bi=b->i,*bj=b->j;
853: PetscInt *a2anew,k,k1,jmin,jmax,*jl,*il,vj,nexti,ili;
854: MatScalar *ba = b->a,*aa,*ap;
855: MatScalar *u,*diag,*rtmp,*rtmp_ptr,dk[4],uik[4];
856: PetscReal shift = info->shiftamount;
857: PetscBool allowzeropivot,zeropivotdetected;
860: allowzeropivot = PetscNot(A->erroriffailure);
862: /* initialization */
863: /* il and jl record the first nonzero element in each row of the accessing
864: window U(0:k, k:mbs-1).
865: jl: list of rows to be added to uneliminated rows
866: i>= k: jl(i) is the first row to be added to row i
867: i< k: jl(i) is the row following row i in some list of rows
868: jl(i) = mbs indicates the end of a list
869: il(i): points to the first nonzero element in columns k,...,mbs-1 of
870: row i of U */
871: PetscCalloc1(4*mbs,&rtmp);
872: PetscMalloc2(mbs,&il,mbs,&jl);
873: il[0] = 0;
874: for (i=0; i<mbs; i++) jl[i] = mbs;
875:
876: ISGetIndices(perm,&perm_ptr);
878: /* check permutation */
879: if (!a->permute) {
880: ai = a->i; aj = a->j; aa = a->a;
881: } else {
882: ai = a->inew; aj = a->jnew;
883: PetscMalloc1(4*ai[mbs],&aa);
884: PetscMemcpy(aa,a->a,4*ai[mbs]*sizeof(MatScalar));
885: PetscMalloc1(ai[mbs],&a2anew);
886: PetscMemcpy(a2anew,a->a2anew,(ai[mbs])*sizeof(PetscInt));
888: for (i=0; i<mbs; i++) {
889: jmin = ai[i]; jmax = ai[i+1];
890: for (j=jmin; j<jmax; j++) {
891: while (a2anew[j] != j) {
892: k = a2anew[j]; a2anew[j] = a2anew[k]; a2anew[k] = k;
893: for (k1=0; k1<4; k1++) {
894: dk[k1] = aa[k*4+k1];
895: aa[k*4+k1] = aa[j*4+k1];
896: aa[j*4+k1] = dk[k1];
897: }
898: }
899: /* transform columnoriented blocks that lie in the lower triangle to roworiented blocks */
900: if (i > aj[j]) {
901: /* printf("change orientation, row: %d, col: %d\n",i,aj[j]); */
902: ap = aa + j*4; /* ptr to the beginning of the block */
903: dk[1] = ap[1]; /* swap ap[1] and ap[2] */
904: ap[1] = ap[2];
905: ap[2] = dk[1];
906: }
907: }
908: }
909: PetscFree(a2anew);
910: }
912: /* for each row k */
913: for (k = 0; k<mbs; k++) {
915: /*initialize k-th row with elements nonzero in row perm(k) of A */
916: jmin = ai[perm_ptr[k]]; jmax = ai[perm_ptr[k]+1];
917: ap = aa + jmin*4;
918: for (j = jmin; j < jmax; j++) {
919: vj = perm_ptr[aj[j]]; /* block col. index */
920: rtmp_ptr = rtmp + vj*4;
921: for (i=0; i<4; i++) *rtmp_ptr++ = *ap++;
922: }
924: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
925: PetscMemcpy(dk,rtmp+k*4,4*sizeof(MatScalar));
926: i = jl[k]; /* first row to be added to k_th row */
928: while (i < k) {
929: nexti = jl[i]; /* next row to be added to k_th row */
931: /* compute multiplier */
932: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
934: /* uik = -inv(Di)*U_bar(i,k): - ba[ili]*ba[i] */
935: diag = ba + i*4;
936: u = ba + ili*4;
937: uik[0] = -(diag[0]*u[0] + diag[2]*u[1]);
938: uik[1] = -(diag[1]*u[0] + diag[3]*u[1]);
939: uik[2] = -(diag[0]*u[2] + diag[2]*u[3]);
940: uik[3] = -(diag[1]*u[2] + diag[3]*u[3]);
942: /* update D(k) += -U(i,k)^T * U_bar(i,k): dk += uik*ba[ili] */
943: dk[0] += uik[0]*u[0] + uik[1]*u[1];
944: dk[1] += uik[2]*u[0] + uik[3]*u[1];
945: dk[2] += uik[0]*u[2] + uik[1]*u[3];
946: dk[3] += uik[2]*u[2] + uik[3]*u[3];
948: PetscLogFlops(16.0*2.0);
950: /* update -U(i,k): ba[ili] = uik */
951: PetscMemcpy(ba+ili*4,uik,4*sizeof(MatScalar));
953: /* add multiple of row i to k-th row ... */
954: jmin = ili + 1; jmax = bi[i+1];
955: if (jmin < jmax) {
956: for (j=jmin; j<jmax; j++) {
957: /* rtmp += -U(i,k)^T * U_bar(i,j): rtmp[bj[j]] += uik*ba[j]; */
958: rtmp_ptr = rtmp + bj[j]*4;
959: u = ba + j*4;
960: rtmp_ptr[0] += uik[0]*u[0] + uik[1]*u[1];
961: rtmp_ptr[1] += uik[2]*u[0] + uik[3]*u[1];
962: rtmp_ptr[2] += uik[0]*u[2] + uik[1]*u[3];
963: rtmp_ptr[3] += uik[2]*u[2] + uik[3]*u[3];
964: }
965: PetscLogFlops(16.0*(jmax-jmin));
967: /* ... add i to row list for next nonzero entry */
968: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
969: j = bj[jmin];
970: jl[i] = jl[j]; jl[j] = i; /* update jl */
971: }
972: i = nexti;
973: }
975: /* save nonzero entries in k-th row of U ... */
977: /* invert diagonal block */
978: diag = ba+k*4;
979: PetscMemcpy(diag,dk,4*sizeof(MatScalar));
980: PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);
981: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
983: jmin = bi[k]; jmax = bi[k+1];
984: if (jmin < jmax) {
985: for (j=jmin; j<jmax; j++) {
986: vj = bj[j]; /* block col. index of U */
987: u = ba + j*4;
988: rtmp_ptr = rtmp + vj*4;
989: for (k1=0; k1<4; k1++) {
990: *u++ = *rtmp_ptr;
991: *rtmp_ptr++ = 0.0;
992: }
993: }
995: /* ... add k to row list for first nonzero entry in k-th row */
996: il[k] = jmin;
997: i = bj[jmin];
998: jl[k] = jl[i]; jl[i] = k;
999: }
1000: }
1002: PetscFree(rtmp);
1003: PetscFree2(il,jl);
1004: if (a->permute) {
1005: PetscFree(aa);
1006: }
1007: ISRestoreIndices(perm,&perm_ptr);
1009: C->ops->solve = MatSolve_SeqSBAIJ_2_inplace;
1010: C->ops->solvetranspose = MatSolve_SeqSBAIJ_2_inplace;
1011: C->assembled = PETSC_TRUE;
1012: C->preallocated = PETSC_TRUE;
1014: PetscLogFlops(1.3333*8*b->mbs); /* from inverting diagonal blocks */
1015: return(0);
1016: }
1018: /*
1019: Version for when blocks are 2 by 2 Using natural ordering
1020: */
1021: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info)
1022: {
1023: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ*)C->data;
1025: PetscInt i,j,mbs=a->mbs,*bi=b->i,*bj=b->j;
1026: PetscInt *ai,*aj,k,k1,jmin,jmax,*jl,*il,vj,nexti,ili;
1027: MatScalar *ba = b->a,*aa,*ap,dk[8],uik[8];
1028: MatScalar *u,*diag,*rtmp,*rtmp_ptr;
1029: PetscReal shift = info->shiftamount;
1030: PetscBool allowzeropivot,zeropivotdetected;
1033: allowzeropivot = PetscNot(A->erroriffailure);
1035: /* initialization */
1036: /* il and jl record the first nonzero element in each row of the accessing
1037: window U(0:k, k:mbs-1).
1038: jl: list of rows to be added to uneliminated rows
1039: i>= k: jl(i) is the first row to be added to row i
1040: i< k: jl(i) is the row following row i in some list of rows
1041: jl(i) = mbs indicates the end of a list
1042: il(i): points to the first nonzero element in columns k,...,mbs-1 of
1043: row i of U */
1044: PetscCalloc1(4*mbs,&rtmp);
1045: PetscMalloc2(mbs,&il,mbs,&jl);
1046: il[0] = 0;
1047: for (i=0; i<mbs; i++) jl[i] = mbs;
1048:
1049: ai = a->i; aj = a->j; aa = a->a;
1051: /* for each row k */
1052: for (k = 0; k<mbs; k++) {
1054: /*initialize k-th row with elements nonzero in row k of A */
1055: jmin = ai[k]; jmax = ai[k+1];
1056: ap = aa + jmin*4;
1057: for (j = jmin; j < jmax; j++) {
1058: vj = aj[j]; /* block col. index */
1059: rtmp_ptr = rtmp + vj*4;
1060: for (i=0; i<4; i++) *rtmp_ptr++ = *ap++;
1061: }
1063: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
1064: PetscMemcpy(dk,rtmp+k*4,4*sizeof(MatScalar));
1065: i = jl[k]; /* first row to be added to k_th row */
1067: while (i < k) {
1068: nexti = jl[i]; /* next row to be added to k_th row */
1070: /* compute multiplier */
1071: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1073: /* uik = -inv(Di)*U_bar(i,k): - ba[ili]*ba[i] */
1074: diag = ba + i*4;
1075: u = ba + ili*4;
1076: uik[0] = -(diag[0]*u[0] + diag[2]*u[1]);
1077: uik[1] = -(diag[1]*u[0] + diag[3]*u[1]);
1078: uik[2] = -(diag[0]*u[2] + diag[2]*u[3]);
1079: uik[3] = -(diag[1]*u[2] + diag[3]*u[3]);
1081: /* update D(k) += -U(i,k)^T * U_bar(i,k): dk += uik*ba[ili] */
1082: dk[0] += uik[0]*u[0] + uik[1]*u[1];
1083: dk[1] += uik[2]*u[0] + uik[3]*u[1];
1084: dk[2] += uik[0]*u[2] + uik[1]*u[3];
1085: dk[3] += uik[2]*u[2] + uik[3]*u[3];
1087: PetscLogFlops(16.0*2.0);
1089: /* update -U(i,k): ba[ili] = uik */
1090: PetscMemcpy(ba+ili*4,uik,4*sizeof(MatScalar));
1092: /* add multiple of row i to k-th row ... */
1093: jmin = ili + 1; jmax = bi[i+1];
1094: if (jmin < jmax) {
1095: for (j=jmin; j<jmax; j++) {
1096: /* rtmp += -U(i,k)^T * U_bar(i,j): rtmp[bj[j]] += uik*ba[j]; */
1097: rtmp_ptr = rtmp + bj[j]*4;
1098: u = ba + j*4;
1099: rtmp_ptr[0] += uik[0]*u[0] + uik[1]*u[1];
1100: rtmp_ptr[1] += uik[2]*u[0] + uik[3]*u[1];
1101: rtmp_ptr[2] += uik[0]*u[2] + uik[1]*u[3];
1102: rtmp_ptr[3] += uik[2]*u[2] + uik[3]*u[3];
1103: }
1104: PetscLogFlops(16.0*(jmax-jmin));
1106: /* ... add i to row list for next nonzero entry */
1107: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
1108: j = bj[jmin];
1109: jl[i] = jl[j]; jl[j] = i; /* update jl */
1110: }
1111: i = nexti;
1112: }
1114: /* save nonzero entries in k-th row of U ... */
1116: /* invert diagonal block */
1117: diag = ba+k*4;
1118: PetscMemcpy(diag,dk,4*sizeof(MatScalar));
1119: PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);
1120: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1122: jmin = bi[k]; jmax = bi[k+1];
1123: if (jmin < jmax) {
1124: for (j=jmin; j<jmax; j++) {
1125: vj = bj[j]; /* block col. index of U */
1126: u = ba + j*4;
1127: rtmp_ptr = rtmp + vj*4;
1128: for (k1=0; k1<4; k1++) {
1129: *u++ = *rtmp_ptr;
1130: *rtmp_ptr++ = 0.0;
1131: }
1132: }
1134: /* ... add k to row list for first nonzero entry in k-th row */
1135: il[k] = jmin;
1136: i = bj[jmin];
1137: jl[k] = jl[i]; jl[i] = k;
1138: }
1139: }
1141: PetscFree(rtmp);
1142: PetscFree2(il,jl);
1144: C->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1145: C->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1146: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1147: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1148: C->assembled = PETSC_TRUE;
1149: C->preallocated = PETSC_TRUE;
1151: PetscLogFlops(1.3333*8*b->mbs); /* from inverting diagonal blocks */
1152: return(0);
1153: }
1155: /*
1156: Numeric U^T*D*U factorization for SBAIJ format.
1157: Version for blocks are 1 by 1.
1158: */
1159: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace(Mat C,Mat A,const MatFactorInfo *info)
1160: {
1161: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data,*b=(Mat_SeqSBAIJ*)C->data;
1162: IS ip=b->row;
1164: const PetscInt *ai,*aj,*rip;
1165: PetscInt *a2anew,i,j,mbs=a->mbs,*bi=b->i,*bj=b->j,*bcol;
1166: PetscInt k,jmin,jmax,*jl,*il,col,nexti,ili,nz;
1167: MatScalar *rtmp,*ba=b->a,*bval,*aa,dk,uikdi;
1168: PetscReal rs;
1169: FactorShiftCtx sctx;
1172: /* MatPivotSetUp(): initialize shift context sctx */
1173: PetscMemzero(&sctx,sizeof(FactorShiftCtx));
1175: ISGetIndices(ip,&rip);
1176: if (!a->permute) {
1177: ai = a->i; aj = a->j; aa = a->a;
1178: } else {
1179: ai = a->inew; aj = a->jnew;
1180: nz = ai[mbs];
1181: PetscMalloc1(nz,&aa);
1182: a2anew = a->a2anew;
1183: bval = a->a;
1184: for (j=0; j<nz; j++) {
1185: aa[a2anew[j]] = *(bval++);
1186: }
1187: }
1189: /* initialization */
1190: /* il and jl record the first nonzero element in each row of the accessing
1191: window U(0:k, k:mbs-1).
1192: jl: list of rows to be added to uneliminated rows
1193: i>= k: jl(i) is the first row to be added to row i
1194: i< k: jl(i) is the row following row i in some list of rows
1195: jl(i) = mbs indicates the end of a list
1196: il(i): points to the first nonzero element in columns k,...,mbs-1 of
1197: row i of U */
1198: PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&jl);
1200: do {
1201: sctx.newshift = PETSC_FALSE;
1202: il[0] = 0;
1203: for (i=0; i<mbs; i++) {
1204: rtmp[i] = 0.0; jl[i] = mbs;
1205: }
1207: for (k = 0; k<mbs; k++) {
1208: /*initialize k-th row by the perm[k]-th row of A */
1209: jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
1210: bval = ba + bi[k];
1211: for (j = jmin; j < jmax; j++) {
1212: col = rip[aj[j]];
1213: rtmp[col] = aa[j];
1214: *bval++ = 0.0; /* for in-place factorization */
1215: }
1217: /* shift the diagonal of the matrix */
1218: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
1220: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1221: dk = rtmp[k];
1222: i = jl[k]; /* first row to be added to k_th row */
1224: while (i < k) {
1225: nexti = jl[i]; /* next row to be added to k_th row */
1227: /* compute multiplier, update diag(k) and U(i,k) */
1228: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1229: uikdi = -ba[ili]*ba[bi[i]]; /* diagonal(k) */
1230: dk += uikdi*ba[ili];
1231: ba[ili] = uikdi; /* -U(i,k) */
1233: /* add multiple of row i to k-th row */
1234: jmin = ili + 1; jmax = bi[i+1];
1235: if (jmin < jmax) {
1236: for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
1237: PetscLogFlops(2.0*(jmax-jmin));
1239: /* update il and jl for row i */
1240: il[i] = jmin;
1241: j = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
1242: }
1243: i = nexti;
1244: }
1246: /* shift the diagonals when zero pivot is detected */
1247: /* compute rs=sum of abs(off-diagonal) */
1248: rs = 0.0;
1249: jmin = bi[k]+1;
1250: nz = bi[k+1] - jmin;
1251: if (nz) {
1252: bcol = bj + jmin;
1253: while (nz--) {
1254: rs += PetscAbsScalar(rtmp[*bcol]);
1255: bcol++;
1256: }
1257: }
1259: sctx.rs = rs;
1260: sctx.pv = dk;
1261: MatPivotCheck(C,A,info,&sctx,k);
1262: if (sctx.newshift) break; /* sctx.shift_amount is updated */
1263: dk = sctx.pv;
1265: /* copy data into U(k,:) */
1266: ba[bi[k]] = 1.0/dk; /* U(k,k) */
1267: jmin = bi[k]+1; jmax = bi[k+1];
1268: if (jmin < jmax) {
1269: for (j=jmin; j<jmax; j++) {
1270: col = bj[j]; ba[j] = rtmp[col]; rtmp[col] = 0.0;
1271: }
1272: /* add the k-th row into il and jl */
1273: il[k] = jmin;
1274: i = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
1275: }
1276: }
1277: } while (sctx.newshift);
1278: PetscFree3(rtmp,il,jl);
1279: if (a->permute) {PetscFree(aa);}
1281: ISRestoreIndices(ip,&rip);
1283: C->ops->solve = MatSolve_SeqSBAIJ_1_inplace;
1284: C->ops->solves = MatSolves_SeqSBAIJ_1_inplace;
1285: C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
1286: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_inplace;
1287: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_inplace;
1288: C->assembled = PETSC_TRUE;
1289: C->preallocated = PETSC_TRUE;
1291: PetscLogFlops(C->rmap->N);
1292: if (sctx.nshift) {
1293: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1294: PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
1295: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1296: PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
1297: }
1298: }
1299: return(0);
1300: }
1302: /*
1303: Version for when blocks are 1 by 1 Using natural ordering under new datastructure
1304: Modified from MatCholeskyFactorNumeric_SeqAIJ()
1305: */
1306: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info)
1307: {
1308: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data;
1309: Mat_SeqSBAIJ *b=(Mat_SeqSBAIJ*)B->data;
1311: PetscInt i,j,mbs=A->rmap->n,*bi=b->i,*bj=b->j,*bdiag=b->diag,*bjtmp;
1312: PetscInt *ai=a->i,*aj=a->j,*ajtmp;
1313: PetscInt k,jmin,jmax,*c2r,*il,col,nexti,ili,nz;
1314: MatScalar *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
1315: FactorShiftCtx sctx;
1316: PetscReal rs;
1317: MatScalar d,*v;
1320: PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&c2r);
1322: /* MatPivotSetUp(): initialize shift context sctx */
1323: PetscMemzero(&sctx,sizeof(FactorShiftCtx));
1325: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
1326: sctx.shift_top = info->zeropivot;
1328: PetscMemzero(rtmp,mbs*sizeof(MatScalar));
1330: for (i=0; i<mbs; i++) {
1331: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
1332: d = (aa)[a->diag[i]];
1333: rtmp[i] += -PetscRealPart(d); /* diagonal entry */
1334: ajtmp = aj + ai[i] + 1; /* exclude diagonal */
1335: v = aa + ai[i] + 1;
1336: nz = ai[i+1] - ai[i] - 1;
1337: for (j=0; j<nz; j++) {
1338: rtmp[i] += PetscAbsScalar(v[j]);
1339: rtmp[ajtmp[j]] += PetscAbsScalar(v[j]);
1340: }
1341: if (PetscRealPart(rtmp[i]) > sctx.shift_top) sctx.shift_top = PetscRealPart(rtmp[i]);
1342: }
1343: sctx.shift_top *= 1.1;
1344: sctx.nshift_max = 5;
1345: sctx.shift_lo = 0.;
1346: sctx.shift_hi = 1.;
1347: }
1349: /* allocate working arrays
1350: c2r: linked list, keep track of pivot rows for a given column. c2r[col]: head of the list for a given col
1351: 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
1352: */
1353: do {
1354: sctx.newshift = PETSC_FALSE;
1356: for (i=0; i<mbs; i++) c2r[i] = mbs;
1357: if (mbs) il[0] = 0;
1359: for (k = 0; k<mbs; k++) {
1360: /* zero rtmp */
1361: nz = bi[k+1] - bi[k];
1362: bjtmp = bj + bi[k];
1363: for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;
1365: /* load in initial unfactored row */
1366: bval = ba + bi[k];
1367: jmin = ai[k]; jmax = ai[k+1];
1368: for (j = jmin; j < jmax; j++) {
1369: col = aj[j];
1370: rtmp[col] = aa[j];
1371: *bval++ = 0.0; /* for in-place factorization */
1372: }
1373: /* shift the diagonal of the matrix: ZeropivotApply() */
1374: rtmp[k] += sctx.shift_amount; /* shift the diagonal of the matrix */
1376: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1377: dk = rtmp[k];
1378: i = c2r[k]; /* first row to be added to k_th row */
1380: while (i < k) {
1381: nexti = c2r[i]; /* next row to be added to k_th row */
1383: /* compute multiplier, update diag(k) and U(i,k) */
1384: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1385: uikdi = -ba[ili]*ba[bdiag[i]]; /* diagonal(k) */
1386: dk += uikdi*ba[ili]; /* update diag[k] */
1387: ba[ili] = uikdi; /* -U(i,k) */
1389: /* add multiple of row i to k-th row */
1390: jmin = ili + 1; jmax = bi[i+1];
1391: if (jmin < jmax) {
1392: for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
1393: /* update il and c2r for row i */
1394: il[i] = jmin;
1395: j = bj[jmin]; c2r[i] = c2r[j]; c2r[j] = i;
1396: }
1397: i = nexti;
1398: }
1400: /* copy data into U(k,:) */
1401: rs = 0.0;
1402: jmin = bi[k]; jmax = bi[k+1]-1;
1403: if (jmin < jmax) {
1404: for (j=jmin; j<jmax; j++) {
1405: col = bj[j]; ba[j] = rtmp[col]; rs += PetscAbsScalar(ba[j]);
1406: }
1407: /* add the k-th row into il and c2r */
1408: il[k] = jmin;
1409: i = bj[jmin]; c2r[k] = c2r[i]; c2r[i] = k;
1410: }
1412: sctx.rs = rs;
1413: sctx.pv = dk;
1414: MatPivotCheck(B,A,info,&sctx,k);
1415: if (sctx.newshift) break;
1416: dk = sctx.pv;
1418: ba[bdiag[k]] = 1.0/dk; /* U(k,k) */
1419: }
1420: } while (sctx.newshift);
1422: PetscFree3(rtmp,il,c2r);
1424: B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1425: B->ops->solves = MatSolves_SeqSBAIJ_1;
1426: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1427: B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
1428: B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;
1430: B->assembled = PETSC_TRUE;
1431: B->preallocated = PETSC_TRUE;
1433: PetscLogFlops(B->rmap->n);
1435: /* MatPivotView() */
1436: if (sctx.nshift) {
1437: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1438: PetscInfo4(A,"number of shift_pd tries %D, 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);
1439: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1440: PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
1441: } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
1442: PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %g\n",sctx.nshift,(double)info->shiftamount);
1443: }
1444: }
1445: return(0);
1446: }
1448: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo *info)
1449: {
1450: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data,*b=(Mat_SeqSBAIJ*)C->data;
1452: PetscInt i,j,mbs = a->mbs;
1453: PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
1454: PetscInt k,jmin,*jl,*il,nexti,ili,*acol,*bcol,nz;
1455: MatScalar *rtmp,*ba=b->a,*aa=a->a,dk,uikdi,*aval,*bval;
1456: PetscReal rs;
1457: FactorShiftCtx sctx;
1460: /* MatPivotSetUp(): initialize shift context sctx */
1461: PetscMemzero(&sctx,sizeof(FactorShiftCtx));
1463: /* initialization */
1464: /* il and jl record the first nonzero element in each row of the accessing
1465: window U(0:k, k:mbs-1).
1466: jl: list of rows to be added to uneliminated rows
1467: i>= k: jl(i) is the first row to be added to row i
1468: i< k: jl(i) is the row following row i in some list of rows
1469: jl(i) = mbs indicates the end of a list
1470: il(i): points to the first nonzero element in U(i,k:mbs-1)
1471: */
1472: PetscMalloc1(mbs,&rtmp);
1473: PetscMalloc2(mbs,&il,mbs,&jl);
1475: do {
1476: sctx.newshift = PETSC_FALSE;
1477: il[0] = 0;
1478: for (i=0; i<mbs; i++) {
1479: rtmp[i] = 0.0; jl[i] = mbs;
1480: }
1482: for (k = 0; k<mbs; k++) {
1483: /*initialize k-th row with elements nonzero in row perm(k) of A */
1484: nz = ai[k+1] - ai[k];
1485: acol = aj + ai[k];
1486: aval = aa + ai[k];
1487: bval = ba + bi[k];
1488: while (nz--) {
1489: rtmp[*acol++] = *aval++;
1490: *bval++ = 0.0; /* for in-place factorization */
1491: }
1493: /* shift the diagonal of the matrix */
1494: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
1496: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1497: dk = rtmp[k];
1498: i = jl[k]; /* first row to be added to k_th row */
1500: while (i < k) {
1501: nexti = jl[i]; /* next row to be added to k_th row */
1502: /* compute multiplier, update D(k) and U(i,k) */
1503: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1504: uikdi = -ba[ili]*ba[bi[i]];
1505: dk += uikdi*ba[ili];
1506: ba[ili] = uikdi; /* -U(i,k) */
1508: /* add multiple of row i to k-th row ... */
1509: jmin = ili + 1;
1510: nz = bi[i+1] - jmin;
1511: if (nz > 0) {
1512: bcol = bj + jmin;
1513: bval = ba + jmin;
1514: PetscLogFlops(2.0*nz);
1515: while (nz--) rtmp[*bcol++] += uikdi*(*bval++);
1517: /* update il and jl for i-th row */
1518: il[i] = jmin;
1519: j = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
1520: }
1521: i = nexti;
1522: }
1524: /* shift the diagonals when zero pivot is detected */
1525: /* compute rs=sum of abs(off-diagonal) */
1526: rs = 0.0;
1527: jmin = bi[k]+1;
1528: nz = bi[k+1] - jmin;
1529: if (nz) {
1530: bcol = bj + jmin;
1531: while (nz--) {
1532: rs += PetscAbsScalar(rtmp[*bcol]);
1533: bcol++;
1534: }
1535: }
1537: sctx.rs = rs;
1538: sctx.pv = dk;
1539: MatPivotCheck(C,A,info,&sctx,k);
1540: if (sctx.newshift) break; /* sctx.shift_amount is updated */
1541: dk = sctx.pv;
1543: /* copy data into U(k,:) */
1544: ba[bi[k]] = 1.0/dk;
1545: jmin = bi[k]+1;
1546: nz = bi[k+1] - jmin;
1547: if (nz) {
1548: bcol = bj + jmin;
1549: bval = ba + jmin;
1550: while (nz--) {
1551: *bval++ = rtmp[*bcol];
1552: rtmp[*bcol++] = 0.0;
1553: }
1554: /* add k-th row into il and jl */
1555: il[k] = jmin;
1556: i = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
1557: }
1558: } /* end of for (k = 0; k<mbs; k++) */
1559: } while (sctx.newshift);
1560: PetscFree(rtmp);
1561: PetscFree2(il,jl);
1563: C->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1564: C->ops->solves = MatSolves_SeqSBAIJ_1_inplace;
1565: C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1566: C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1567: C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1569: C->assembled = PETSC_TRUE;
1570: C->preallocated = PETSC_TRUE;
1572: PetscLogFlops(C->rmap->N);
1573: if (sctx.nshift) {
1574: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1575: PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
1576: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1577: PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
1578: }
1579: }
1580: return(0);
1581: }
1583: PetscErrorCode MatCholeskyFactor_SeqSBAIJ(Mat A,IS perm,const MatFactorInfo *info)
1584: {
1586: Mat C;
1589: MatGetFactor(A,"petsc",MAT_FACTOR_CHOLESKY,&C);
1590: MatCholeskyFactorSymbolic(C,A,perm,info);
1591: MatCholeskyFactorNumeric(C,A,info);
1593: A->ops->solve = C->ops->solve;
1594: A->ops->solvetranspose = C->ops->solvetranspose;
1596: MatHeaderMerge(A,&C);
1597: return(0);
1598: }