Actual source code: sbaijfact.c

petsc-master 2020-11-28
Report Typos and Errors

  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:       PetscArraycpy(jutmp,ju,iu[k]);
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 = NULL;
177:   b->ilen = NULL;
178:   b->imax = NULL;
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,&current_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      = NULL;
340:   b->imax      = NULL;
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,&current_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 = NULL;
516:   b->ilen = NULL;
517:   b->imax = NULL;
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;

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:     PetscArraycpy(aa,a->a,bs2*ai[mbs]);
586:     PetscMalloc1(ai[mbs],&a2anew);
587:     PetscArraycpy(a2anew,a->a2anew,ai[mbs]);

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:           ap = aa + j*bs2;                     /* ptr to the beginning of j-th block of aa */
603:           for (k=0; k<bs2; k++) dk[k] = ap[k]; /* dk <- j-th block of aa */
604:           for (k=0; k<bs; k++) {               /* j-th block of aa <- dk^T */
605:             for (k1=0; k1<bs; k1++) *ap++ = dk[k + bs*k1];
606:           }
607:         }
608:       }
609:     }
610:     PetscFree(a2anew);
611:   }

613:   /* for each row k */
614:   for (k = 0; k<mbs; k++) {

616:     /*initialize k-th row with elements nonzero in row perm(k) of A */
617:     jmin = ai[perm_ptr[k]]; jmax = ai[perm_ptr[k]+1];

619:     ap = aa + jmin*bs2;
620:     for (j = jmin; j < jmax; j++) {
621:       vj       = perm_ptr[aj[j]];   /* block col. index */
622:       rtmp_ptr = rtmp + vj*bs2;
623:       for (i=0; i<bs2; i++) *rtmp_ptr++ = *ap++;
624:     }

626:     /* modify k-th row by adding in those rows i with U(i,k) != 0 */
627:     PetscArraycpy(dk,rtmp+k*bs2,bs2);
628:     i    = jl[k]; /* first row to be added to k_th row  */

630:     while (i < k) {
631:       nexti = jl[i]; /* next row to be added to k_th row */

633:       /* compute multiplier */
634:       ili = il[i];  /* index of first nonzero element in U(i,k:bms-1) */

636:       /* uik = -inv(Di)*U_bar(i,k) */
637:       diag = ba + i*bs2;
638:       u    = ba + ili*bs2;
639:       PetscArrayzero(uik,bs2);
640:       PetscKernel_A_gets_A_minus_B_times_C(bs,uik,diag,u);

642:       /* update D(k) += -U(i,k)^T * U_bar(i,k) */
643:       PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,dk,uik,u);
644:       PetscLogFlops(4.0*bs*bs2);

646:       /* update -U(i,k) */
647:       PetscArraycpy(ba+ili*bs2,uik,bs2);

649:       /* add multiple of row i to k-th row ... */
650:       jmin = ili + 1; jmax = bi[i+1];
651:       if (jmin < jmax) {
652:         for (j=jmin; j<jmax; j++) {
653:           /* rtmp += -U(i,k)^T * U_bar(i,j) */
654:           rtmp_ptr = rtmp + bj[j]*bs2;
655:           u        = ba + j*bs2;
656:           PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,rtmp_ptr,uik,u);
657:         }
658:         PetscLogFlops(2.0*bs*bs2*(jmax-jmin));

660:         /* ... add i to row list for next nonzero entry */
661:         il[i] = jmin;             /* update il(i) in column k+1, ... mbs-1 */
662:         j     = bj[jmin];
663:         jl[i] = jl[j]; jl[j] = i; /* update jl */
664:       }
665:       i = nexti;
666:     }

668:     /* save nonzero entries in k-th row of U ... */

670:     /* invert diagonal block */
671:     diag = ba+k*bs2;
672:     PetscArraycpy(diag,dk,bs2);

674:     PetscKernel_A_gets_inverse_A(bs,diag,pivots,work,allowzeropivot,&zeropivotdetected);
675:     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

677:     jmin = bi[k]; jmax = bi[k+1];
678:     if (jmin < jmax) {
679:       for (j=jmin; j<jmax; j++) {
680:         vj       = bj[j];      /* block col. index of U */
681:         u        = ba + j*bs2;
682:         rtmp_ptr = rtmp + vj*bs2;
683:         for (k1=0; k1<bs2; k1++) {
684:           *u++        = *rtmp_ptr;
685:           *rtmp_ptr++ = 0.0;
686:         }
687:       }

689:       /* ... add k to row list for first nonzero entry in k-th row */
690:       il[k] = jmin;
691:       i     = bj[jmin];
692:       jl[k] = jl[i]; jl[i] = k;
693:     }
694:   }

696:   PetscFree(rtmp);
697:   PetscFree2(il,jl);
698:   PetscFree3(dk,uik,work);
699:   PetscFree(pivots);
700:   if (a->permute) {
701:     PetscFree(aa);
702:   }

704:   ISRestoreIndices(perm,&perm_ptr);

706:   C->ops->solve          = MatSolve_SeqSBAIJ_N_inplace;
707:   C->ops->solvetranspose = MatSolve_SeqSBAIJ_N_inplace;
708:   C->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_N_inplace;
709:   C->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_N_inplace;

711:   C->assembled    = PETSC_TRUE;
712:   C->preallocated = PETSC_TRUE;

714:   PetscLogFlops(1.3333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
715:   return(0);
716: }

718: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info)
719: {
720:   Mat_SeqSBAIJ   *a = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ*)C->data;
722:   PetscInt       i,j,mbs=a->mbs,*bi=b->i,*bj=b->j;
723:   PetscInt       *ai,*aj,k,k1,jmin,jmax,*jl,*il,vj,nexti,ili;
724:   PetscInt       bs  =A->rmap->bs,bs2 = a->bs2;
725:   MatScalar      *ba = b->a,*aa,*ap,*dk,*uik;
726:   MatScalar      *u,*diag,*rtmp,*rtmp_ptr;
727:   MatScalar      *work;
728:   PetscInt       *pivots;
729:   PetscBool      allowzeropivot,zeropivotdetected;

732:   PetscCalloc1(bs2*mbs,&rtmp);
733:   PetscMalloc2(mbs,&il,mbs,&jl);
734:   il[0] = 0;
735:   for (i=0; i<mbs; i++) jl[i] = mbs;

737:   PetscMalloc3(bs2,&dk,bs2,&uik,bs,&work);
738:   PetscMalloc1(bs,&pivots);
739:   allowzeropivot = PetscNot(A->erroriffailure);

741:   ai = a->i; aj = a->j; aa = a->a;

743:   /* for each row k */
744:   for (k = 0; k<mbs; k++) {

746:     /*initialize k-th row with elements nonzero in row k of A */
747:     jmin = ai[k]; jmax = ai[k+1];
748:     ap   = aa + jmin*bs2;
749:     for (j = jmin; j < jmax; j++) {
750:       vj       = aj[j];   /* block col. index */
751:       rtmp_ptr = rtmp + vj*bs2;
752:       for (i=0; i<bs2; i++) *rtmp_ptr++ = *ap++;
753:     }

755:     /* modify k-th row by adding in those rows i with U(i,k) != 0 */
756:     PetscArraycpy(dk,rtmp+k*bs2,bs2);
757:     i    = jl[k]; /* first row to be added to k_th row  */

759:     while (i < k) {
760:       nexti = jl[i]; /* next row to be added to k_th row */

762:       /* compute multiplier */
763:       ili = il[i];  /* index of first nonzero element in U(i,k:bms-1) */

765:       /* uik = -inv(Di)*U_bar(i,k) */
766:       diag = ba + i*bs2;
767:       u    = ba + ili*bs2;
768:       PetscArrayzero(uik,bs2);
769:       PetscKernel_A_gets_A_minus_B_times_C(bs,uik,diag,u);

771:       /* update D(k) += -U(i,k)^T * U_bar(i,k) */
772:       PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,dk,uik,u);
773:       PetscLogFlops(2.0*bs*bs2);

775:       /* update -U(i,k) */
776:       PetscArraycpy(ba+ili*bs2,uik,bs2);

778:       /* add multiple of row i to k-th row ... */
779:       jmin = ili + 1; jmax = bi[i+1];
780:       if (jmin < jmax) {
781:         for (j=jmin; j<jmax; j++) {
782:           /* rtmp += -U(i,k)^T * U_bar(i,j) */
783:           rtmp_ptr = rtmp + bj[j]*bs2;
784:           u        = ba + j*bs2;
785:           PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,rtmp_ptr,uik,u);
786:         }
787:         PetscLogFlops(2.0*bs*bs2*(jmax-jmin));

789:         /* ... add i to row list for next nonzero entry */
790:         il[i] = jmin;             /* update il(i) in column k+1, ... mbs-1 */
791:         j     = bj[jmin];
792:         jl[i] = jl[j]; jl[j] = i; /* update jl */
793:       }
794:       i = nexti;
795:     }

797:     /* save nonzero entries in k-th row of U ... */

799:     /* invert diagonal block */
800:     diag = ba+k*bs2;
801:     PetscArraycpy(diag,dk,bs2);

803:     PetscKernel_A_gets_inverse_A(bs,diag,pivots,work,allowzeropivot,&zeropivotdetected);
804:     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

806:     jmin = bi[k]; jmax = bi[k+1];
807:     if (jmin < jmax) {
808:       for (j=jmin; j<jmax; j++) {
809:         vj       = bj[j];      /* block col. index of U */
810:         u        = ba + j*bs2;
811:         rtmp_ptr = rtmp + vj*bs2;
812:         for (k1=0; k1<bs2; k1++) {
813:           *u++        = *rtmp_ptr;
814:           *rtmp_ptr++ = 0.0;
815:         }
816:       }

818:       /* ... add k to row list for first nonzero entry in k-th row */
819:       il[k] = jmin;
820:       i     = bj[jmin];
821:       jl[k] = jl[i]; jl[i] = k;
822:     }
823:   }

825:   PetscFree(rtmp);
826:   PetscFree2(il,jl);
827:   PetscFree3(dk,uik,work);
828:   PetscFree(pivots);

830:   C->ops->solve          = MatSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
831:   C->ops->solvetranspose = MatSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
832:   C->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
833:   C->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_N_NaturalOrdering_inplace;
834:   C->assembled           = PETSC_TRUE;
835:   C->preallocated        = PETSC_TRUE;

837:   PetscLogFlops(1.3333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
838:   return(0);
839: }

841: /*
842:     Numeric U^T*D*U factorization for SBAIJ format. Modified from SNF of YSMP.
843:     Version for blocks 2 by 2.
844: */
845: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat C,Mat A,const MatFactorInfo *info)
846: {
847:   Mat_SeqSBAIJ   *a   = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ*)C->data;
848:   IS             perm = b->row;
850:   const PetscInt *ai,*aj,*perm_ptr;
851:   PetscInt       i,j,mbs=a->mbs,*bi=b->i,*bj=b->j;
852:   PetscInt       *a2anew,k,k1,jmin,jmax,*jl,*il,vj,nexti,ili;
853:   MatScalar      *ba = b->a,*aa,*ap;
854:   MatScalar      *u,*diag,*rtmp,*rtmp_ptr,dk[4],uik[4];
855:   PetscReal      shift = info->shiftamount;
856:   PetscBool      allowzeropivot,zeropivotdetected;

859:   allowzeropivot = PetscNot(A->erroriffailure);

861:   /* initialization */
862:   /* il and jl record the first nonzero element in each row of the accessing
863:      window U(0:k, k:mbs-1).
864:      jl:    list of rows to be added to uneliminated rows
865:             i>= k: jl(i) is the first row to be added to row i
866:             i<  k: jl(i) is the row following row i in some list of rows
867:             jl(i) = mbs indicates the end of a list
868:      il(i): points to the first nonzero element in columns k,...,mbs-1 of
869:             row i of U */
870:   PetscCalloc1(4*mbs,&rtmp);
871:   PetscMalloc2(mbs,&il,mbs,&jl);
872:   il[0] = 0;
873:   for (i=0; i<mbs; i++) jl[i] = mbs;

875:   ISGetIndices(perm,&perm_ptr);

877:   /* check permutation */
878:   if (!a->permute) {
879:     ai = a->i; aj = a->j; aa = a->a;
880:   } else {
881:     ai   = a->inew; aj = a->jnew;
882:     PetscMalloc1(4*ai[mbs],&aa);
883:     PetscArraycpy(aa,a->a,4*ai[mbs]);
884:     PetscMalloc1(ai[mbs],&a2anew);
885:     PetscArraycpy(a2anew,a->a2anew,ai[mbs]);

887:     for (i=0; i<mbs; i++) {
888:       jmin = ai[i]; jmax = ai[i+1];
889:       for (j=jmin; j<jmax; j++) {
890:         while (a2anew[j] != j) {
891:           k = a2anew[j]; a2anew[j] = a2anew[k]; a2anew[k] = k;
892:           for (k1=0; k1<4; k1++) {
893:             dk[k1]     = aa[k*4+k1];
894:             aa[k*4+k1] = aa[j*4+k1];
895:             aa[j*4+k1] = dk[k1];
896:           }
897:         }
898:         /* transform columnoriented blocks that lie in the lower triangle to roworiented blocks */
899:         if (i > aj[j]) {
900:           ap    = aa + j*4;  /* ptr to the beginning of the block */
901:           dk[1] = ap[1];     /* swap ap[1] and ap[2] */
902:           ap[1] = ap[2];
903:           ap[2] = dk[1];
904:         }
905:       }
906:     }
907:     PetscFree(a2anew);
908:   }

910:   /* for each row k */
911:   for (k = 0; k<mbs; k++) {

913:     /*initialize k-th row with elements nonzero in row perm(k) of A */
914:     jmin = ai[perm_ptr[k]]; jmax = ai[perm_ptr[k]+1];
915:     ap   = aa + jmin*4;
916:     for (j = jmin; j < jmax; j++) {
917:       vj       = perm_ptr[aj[j]];   /* block col. index */
918:       rtmp_ptr = rtmp + vj*4;
919:       for (i=0; i<4; i++) *rtmp_ptr++ = *ap++;
920:     }

922:     /* modify k-th row by adding in those rows i with U(i,k) != 0 */
923:     PetscArraycpy(dk,rtmp+k*4,4);
924:     i    = jl[k]; /* first row to be added to k_th row  */

926:     while (i < k) {
927:       nexti = jl[i]; /* next row to be added to k_th row */

929:       /* compute multiplier */
930:       ili = il[i];  /* index of first nonzero element in U(i,k:bms-1) */

932:       /* uik = -inv(Di)*U_bar(i,k): - ba[ili]*ba[i] */
933:       diag   = ba + i*4;
934:       u      = ba + ili*4;
935:       uik[0] = -(diag[0]*u[0] + diag[2]*u[1]);
936:       uik[1] = -(diag[1]*u[0] + diag[3]*u[1]);
937:       uik[2] = -(diag[0]*u[2] + diag[2]*u[3]);
938:       uik[3] = -(diag[1]*u[2] + diag[3]*u[3]);

940:       /* update D(k) += -U(i,k)^T * U_bar(i,k): dk += uik*ba[ili] */
941:       dk[0] += uik[0]*u[0] + uik[1]*u[1];
942:       dk[1] += uik[2]*u[0] + uik[3]*u[1];
943:       dk[2] += uik[0]*u[2] + uik[1]*u[3];
944:       dk[3] += uik[2]*u[2] + uik[3]*u[3];

946:       PetscLogFlops(16.0*2.0);

948:       /* update -U(i,k): ba[ili] = uik */
949:       PetscArraycpy(ba+ili*4,uik,4);

951:       /* add multiple of row i to k-th row ... */
952:       jmin = ili + 1; jmax = bi[i+1];
953:       if (jmin < jmax) {
954:         for (j=jmin; j<jmax; j++) {
955:           /* rtmp += -U(i,k)^T * U_bar(i,j): rtmp[bj[j]] += uik*ba[j]; */
956:           rtmp_ptr     = rtmp + bj[j]*4;
957:           u            = ba + j*4;
958:           rtmp_ptr[0] += uik[0]*u[0] + uik[1]*u[1];
959:           rtmp_ptr[1] += uik[2]*u[0] + uik[3]*u[1];
960:           rtmp_ptr[2] += uik[0]*u[2] + uik[1]*u[3];
961:           rtmp_ptr[3] += uik[2]*u[2] + uik[3]*u[3];
962:         }
963:         PetscLogFlops(16.0*(jmax-jmin));

965:         /* ... add i to row list for next nonzero entry */
966:         il[i] = jmin;             /* update il(i) in column k+1, ... mbs-1 */
967:         j     = bj[jmin];
968:         jl[i] = jl[j]; jl[j] = i; /* update jl */
969:       }
970:       i = nexti;
971:     }

973:     /* save nonzero entries in k-th row of U ... */

975:     /* invert diagonal block */
976:     diag = ba+k*4;
977:     PetscArraycpy(diag,dk,4);
978:     PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);
979:     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

981:     jmin = bi[k]; jmax = bi[k+1];
982:     if (jmin < jmax) {
983:       for (j=jmin; j<jmax; j++) {
984:         vj       = bj[j];      /* block col. index of U */
985:         u        = ba + j*4;
986:         rtmp_ptr = rtmp + vj*4;
987:         for (k1=0; k1<4; k1++) {
988:           *u++        = *rtmp_ptr;
989:           *rtmp_ptr++ = 0.0;
990:         }
991:       }

993:       /* ... add k to row list for first nonzero entry in k-th row */
994:       il[k] = jmin;
995:       i     = bj[jmin];
996:       jl[k] = jl[i]; jl[i] = k;
997:     }
998:   }

1000:   PetscFree(rtmp);
1001:   PetscFree2(il,jl);
1002:   if (a->permute) {
1003:     PetscFree(aa);
1004:   }
1005:   ISRestoreIndices(perm,&perm_ptr);

1007:   C->ops->solve          = MatSolve_SeqSBAIJ_2_inplace;
1008:   C->ops->solvetranspose = MatSolve_SeqSBAIJ_2_inplace;
1009:   C->assembled           = PETSC_TRUE;
1010:   C->preallocated        = PETSC_TRUE;

1012:   PetscLogFlops(1.3333*8*b->mbs); /* from inverting diagonal blocks */
1013:   return(0);
1014: }

1016: /*
1017:       Version for when blocks are 2 by 2 Using natural ordering
1018: */
1019: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info)
1020: {
1021:   Mat_SeqSBAIJ   *a = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ*)C->data;
1023:   PetscInt       i,j,mbs=a->mbs,*bi=b->i,*bj=b->j;
1024:   PetscInt       *ai,*aj,k,k1,jmin,jmax,*jl,*il,vj,nexti,ili;
1025:   MatScalar      *ba = b->a,*aa,*ap,dk[8],uik[8];
1026:   MatScalar      *u,*diag,*rtmp,*rtmp_ptr;
1027:   PetscReal      shift = info->shiftamount;
1028:   PetscBool      allowzeropivot,zeropivotdetected;

1031:   allowzeropivot = PetscNot(A->erroriffailure);

1033:   /* initialization */
1034:   /* il and jl record the first nonzero element in each row of the accessing
1035:      window U(0:k, k:mbs-1).
1036:      jl:    list of rows to be added to uneliminated rows
1037:             i>= k: jl(i) is the first row to be added to row i
1038:             i<  k: jl(i) is the row following row i in some list of rows
1039:             jl(i) = mbs indicates the end of a list
1040:      il(i): points to the first nonzero element in columns k,...,mbs-1 of
1041:             row i of U */
1042:   PetscCalloc1(4*mbs,&rtmp);
1043:   PetscMalloc2(mbs,&il,mbs,&jl);
1044:   il[0] = 0;
1045:   for (i=0; i<mbs; i++) jl[i] = mbs;

1047:   ai = a->i; aj = a->j; aa = a->a;

1049:   /* for each row k */
1050:   for (k = 0; k<mbs; k++) {

1052:     /*initialize k-th row with elements nonzero in row k of A */
1053:     jmin = ai[k]; jmax = ai[k+1];
1054:     ap   = aa + jmin*4;
1055:     for (j = jmin; j < jmax; j++) {
1056:       vj       = aj[j];   /* block col. index */
1057:       rtmp_ptr = rtmp + vj*4;
1058:       for (i=0; i<4; i++) *rtmp_ptr++ = *ap++;
1059:     }

1061:     /* modify k-th row by adding in those rows i with U(i,k) != 0 */
1062:     PetscArraycpy(dk,rtmp+k*4,4);
1063:     i    = jl[k]; /* first row to be added to k_th row  */

1065:     while (i < k) {
1066:       nexti = jl[i]; /* next row to be added to k_th row */

1068:       /* compute multiplier */
1069:       ili = il[i];  /* index of first nonzero element in U(i,k:bms-1) */

1071:       /* uik = -inv(Di)*U_bar(i,k): - ba[ili]*ba[i] */
1072:       diag   = ba + i*4;
1073:       u      = ba + ili*4;
1074:       uik[0] = -(diag[0]*u[0] + diag[2]*u[1]);
1075:       uik[1] = -(diag[1]*u[0] + diag[3]*u[1]);
1076:       uik[2] = -(diag[0]*u[2] + diag[2]*u[3]);
1077:       uik[3] = -(diag[1]*u[2] + diag[3]*u[3]);

1079:       /* update D(k) += -U(i,k)^T * U_bar(i,k): dk += uik*ba[ili] */
1080:       dk[0] += uik[0]*u[0] + uik[1]*u[1];
1081:       dk[1] += uik[2]*u[0] + uik[3]*u[1];
1082:       dk[2] += uik[0]*u[2] + uik[1]*u[3];
1083:       dk[3] += uik[2]*u[2] + uik[3]*u[3];

1085:       PetscLogFlops(16.0*2.0);

1087:       /* update -U(i,k): ba[ili] = uik */
1088:       PetscArraycpy(ba+ili*4,uik,4);

1090:       /* add multiple of row i to k-th row ... */
1091:       jmin = ili + 1; jmax = bi[i+1];
1092:       if (jmin < jmax) {
1093:         for (j=jmin; j<jmax; j++) {
1094:           /* rtmp += -U(i,k)^T * U_bar(i,j): rtmp[bj[j]] += uik*ba[j]; */
1095:           rtmp_ptr     = rtmp + bj[j]*4;
1096:           u            = ba + j*4;
1097:           rtmp_ptr[0] += uik[0]*u[0] + uik[1]*u[1];
1098:           rtmp_ptr[1] += uik[2]*u[0] + uik[3]*u[1];
1099:           rtmp_ptr[2] += uik[0]*u[2] + uik[1]*u[3];
1100:           rtmp_ptr[3] += uik[2]*u[2] + uik[3]*u[3];
1101:         }
1102:         PetscLogFlops(16.0*(jmax-jmin));

1104:         /* ... add i to row list for next nonzero entry */
1105:         il[i] = jmin;             /* update il(i) in column k+1, ... mbs-1 */
1106:         j     = bj[jmin];
1107:         jl[i] = jl[j]; jl[j] = i; /* update jl */
1108:       }
1109:       i = nexti;
1110:     }

1112:     /* save nonzero entries in k-th row of U ... */

1114:     /* invert diagonal block */
1115:     diag = ba+k*4;
1116:     PetscArraycpy(diag,dk,4);
1117:     PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);
1118:     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

1120:     jmin = bi[k]; jmax = bi[k+1];
1121:     if (jmin < jmax) {
1122:       for (j=jmin; j<jmax; j++) {
1123:         vj       = bj[j];      /* block col. index of U */
1124:         u        = ba + j*4;
1125:         rtmp_ptr = rtmp + vj*4;
1126:         for (k1=0; k1<4; k1++) {
1127:           *u++        = *rtmp_ptr;
1128:           *rtmp_ptr++ = 0.0;
1129:         }
1130:       }

1132:       /* ... add k to row list for first nonzero entry in k-th row */
1133:       il[k] = jmin;
1134:       i     = bj[jmin];
1135:       jl[k] = jl[i]; jl[i] = k;
1136:     }
1137:   }

1139:   PetscFree(rtmp);
1140:   PetscFree2(il,jl);

1142:   C->ops->solve          = MatSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1143:   C->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1144:   C->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1145:   C->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_2_NaturalOrdering_inplace;
1146:   C->assembled           = PETSC_TRUE;
1147:   C->preallocated        = PETSC_TRUE;

1149:   PetscLogFlops(1.3333*8*b->mbs); /* from inverting diagonal blocks */
1150:   return(0);
1151: }

1153: /*
1154:     Numeric U^T*D*U factorization for SBAIJ format.
1155:     Version for blocks are 1 by 1.
1156: */
1157: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace(Mat C,Mat A,const MatFactorInfo *info)
1158: {
1159:   Mat_SeqSBAIJ   *a=(Mat_SeqSBAIJ*)A->data,*b=(Mat_SeqSBAIJ*)C->data;
1160:   IS             ip=b->row;
1162:   const PetscInt *ai,*aj,*rip;
1163:   PetscInt       *a2anew,i,j,mbs=a->mbs,*bi=b->i,*bj=b->j,*bcol;
1164:   PetscInt       k,jmin,jmax,*jl,*il,col,nexti,ili,nz;
1165:   MatScalar      *rtmp,*ba=b->a,*bval,*aa,dk,uikdi;
1166:   PetscReal      rs;
1167:   FactorShiftCtx sctx;

1170:   /* MatPivotSetUp(): initialize shift context sctx */
1171:   PetscMemzero(&sctx,sizeof(FactorShiftCtx));

1173:   ISGetIndices(ip,&rip);
1174:   if (!a->permute) {
1175:     ai = a->i; aj = a->j; aa = a->a;
1176:   } else {
1177:     ai     = a->inew; aj = a->jnew;
1178:     nz     = ai[mbs];
1179:     PetscMalloc1(nz,&aa);
1180:     a2anew = a->a2anew;
1181:     bval   = a->a;
1182:     for (j=0; j<nz; j++) {
1183:       aa[a2anew[j]] = *(bval++);
1184:     }
1185:   }

1187:   /* initialization */
1188:   /* il and jl record the first nonzero element in each row of the accessing
1189:      window U(0:k, k:mbs-1).
1190:      jl:    list of rows to be added to uneliminated rows
1191:             i>= k: jl(i) is the first row to be added to row i
1192:             i<  k: jl(i) is the row following row i in some list of rows
1193:             jl(i) = mbs indicates the end of a list
1194:      il(i): points to the first nonzero element in columns k,...,mbs-1 of
1195:             row i of U */
1196:   PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&jl);

1198:   do {
1199:     sctx.newshift = PETSC_FALSE;
1200:     il[0] = 0;
1201:     for (i=0; i<mbs; i++) {
1202:       rtmp[i] = 0.0; jl[i] = mbs;
1203:     }

1205:     for (k = 0; k<mbs; k++) {
1206:       /*initialize k-th row by the perm[k]-th row of A */
1207:       jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
1208:       bval = ba + bi[k];
1209:       for (j = jmin; j < jmax; j++) {
1210:         col       = rip[aj[j]];
1211:         rtmp[col] = aa[j];
1212:         *bval++   = 0.0; /* for in-place factorization */
1213:       }

1215:       /* shift the diagonal of the matrix */
1216:       if (sctx.nshift) rtmp[k] += sctx.shift_amount;

1218:       /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1219:       dk = rtmp[k];
1220:       i  = jl[k]; /* first row to be added to k_th row  */

1222:       while (i < k) {
1223:         nexti = jl[i]; /* next row to be added to k_th row */

1225:         /* compute multiplier, update diag(k) and U(i,k) */
1226:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1227:         uikdi   = -ba[ili]*ba[bi[i]]; /* diagonal(k) */
1228:         dk     += uikdi*ba[ili];
1229:         ba[ili] = uikdi; /* -U(i,k) */

1231:         /* add multiple of row i to k-th row */
1232:         jmin = ili + 1; jmax = bi[i+1];
1233:         if (jmin < jmax) {
1234:           for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
1235:           PetscLogFlops(2.0*(jmax-jmin));

1237:           /* update il and jl for row i */
1238:           il[i] = jmin;
1239:           j     = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
1240:         }
1241:         i = nexti;
1242:       }

1244:       /* shift the diagonals when zero pivot is detected */
1245:       /* compute rs=sum of abs(off-diagonal) */
1246:       rs   = 0.0;
1247:       jmin = bi[k]+1;
1248:       nz   = bi[k+1] - jmin;
1249:       if (nz) {
1250:         bcol = bj + jmin;
1251:         while (nz--) {
1252:           rs += PetscAbsScalar(rtmp[*bcol]);
1253:           bcol++;
1254:         }
1255:       }

1257:       sctx.rs = rs;
1258:       sctx.pv = dk;
1259:       MatPivotCheck(C,A,info,&sctx,k);
1260:       if (sctx.newshift) break;    /* sctx.shift_amount is updated */
1261:       dk = sctx.pv;

1263:       /* copy data into U(k,:) */
1264:       ba[bi[k]] = 1.0/dk; /* U(k,k) */
1265:       jmin      = bi[k]+1; jmax = bi[k+1];
1266:       if (jmin < jmax) {
1267:         for (j=jmin; j<jmax; j++) {
1268:           col = bj[j]; ba[j] = rtmp[col]; rtmp[col] = 0.0;
1269:         }
1270:         /* add the k-th row into il and jl */
1271:         il[k] = jmin;
1272:         i     = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
1273:       }
1274:     }
1275:   } while (sctx.newshift);
1276:   PetscFree3(rtmp,il,jl);
1277:   if (a->permute) {PetscFree(aa);}

1279:   ISRestoreIndices(ip,&rip);

1281:   C->ops->solve          = MatSolve_SeqSBAIJ_1_inplace;
1282:   C->ops->solves         = MatSolves_SeqSBAIJ_1_inplace;
1283:   C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
1284:   C->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_inplace;
1285:   C->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_inplace;
1286:   C->assembled           = PETSC_TRUE;
1287:   C->preallocated        = PETSC_TRUE;

1289:   PetscLogFlops(C->rmap->N);
1290:   if (sctx.nshift) {
1291:     if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1292:       PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
1293:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1294:       PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
1295:     }
1296:   }
1297:   return(0);
1298: }

1300: /*
1301:   Version for when blocks are 1 by 1 Using natural ordering under new datastructure
1302:   Modified from MatCholeskyFactorNumeric_SeqAIJ()
1303: */
1304: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info)
1305: {
1306:   Mat_SeqSBAIJ   *a=(Mat_SeqSBAIJ*)A->data;
1307:   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)B->data;
1309:   PetscInt       i,j,mbs=A->rmap->n,*bi=b->i,*bj=b->j,*bdiag=b->diag,*bjtmp;
1310:   PetscInt       *ai=a->i,*aj=a->j,*ajtmp;
1311:   PetscInt       k,jmin,jmax,*c2r,*il,col,nexti,ili,nz;
1312:   MatScalar      *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
1313:   FactorShiftCtx sctx;
1314:   PetscReal      rs;
1315:   MatScalar      d,*v;

1318:   PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&c2r);

1320:   /* MatPivotSetUp(): initialize shift context sctx */
1321:   PetscMemzero(&sctx,sizeof(FactorShiftCtx));

1323:   if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
1324:     sctx.shift_top = info->zeropivot;

1326:     PetscArrayzero(rtmp,mbs);

1328:     for (i=0; i<mbs; i++) {
1329:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
1330:       d        = (aa)[a->diag[i]];
1331:       rtmp[i] += -PetscRealPart(d);  /* diagonal entry */
1332:       ajtmp    = aj + ai[i] + 1;     /* exclude diagonal */
1333:       v        = aa + ai[i] + 1;
1334:       nz       = ai[i+1] - ai[i] - 1;
1335:       for (j=0; j<nz; j++) {
1336:         rtmp[i]        += PetscAbsScalar(v[j]);
1337:         rtmp[ajtmp[j]] += PetscAbsScalar(v[j]);
1338:       }
1339:       if (PetscRealPart(rtmp[i]) > sctx.shift_top) sctx.shift_top = PetscRealPart(rtmp[i]);
1340:     }
1341:     sctx.shift_top *= 1.1;
1342:     sctx.nshift_max = 5;
1343:     sctx.shift_lo   = 0.;
1344:     sctx.shift_hi   = 1.;
1345:   }

1347:   /* allocate working arrays
1348:      c2r: linked list, keep track of pivot rows for a given column. c2r[col]: head of the list for a given col
1349:      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
1350:   */
1351:   do {
1352:     sctx.newshift = PETSC_FALSE;

1354:     for (i=0; i<mbs; i++) c2r[i] = mbs;
1355:     if (mbs) il[0] = 0;

1357:     for (k = 0; k<mbs; k++) {
1358:       /* zero rtmp */
1359:       nz    = bi[k+1] - bi[k];
1360:       bjtmp = bj + bi[k];
1361:       for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

1363:       /* load in initial unfactored row */
1364:       bval = ba + bi[k];
1365:       jmin = ai[k]; jmax = ai[k+1];
1366:       for (j = jmin; j < jmax; j++) {
1367:         col       = aj[j];
1368:         rtmp[col] = aa[j];
1369:         *bval++   = 0.0; /* for in-place factorization */
1370:       }
1371:       /* shift the diagonal of the matrix: ZeropivotApply() */
1372:       rtmp[k] += sctx.shift_amount;  /* shift the diagonal of the matrix */

1374:       /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1375:       dk = rtmp[k];
1376:       i  = c2r[k]; /* first row to be added to k_th row  */

1378:       while (i < k) {
1379:         nexti = c2r[i]; /* next row to be added to k_th row */

1381:         /* compute multiplier, update diag(k) and U(i,k) */
1382:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1383:         uikdi   = -ba[ili]*ba[bdiag[i]]; /* diagonal(k) */
1384:         dk     += uikdi*ba[ili]; /* update diag[k] */
1385:         ba[ili] = uikdi; /* -U(i,k) */

1387:         /* add multiple of row i to k-th row */
1388:         jmin = ili + 1; jmax = bi[i+1];
1389:         if (jmin < jmax) {
1390:           for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
1391:           /* update il and c2r for row i */
1392:           il[i] = jmin;
1393:           j     = bj[jmin]; c2r[i] = c2r[j]; c2r[j] = i;
1394:         }
1395:         i = nexti;
1396:       }

1398:       /* copy data into U(k,:) */
1399:       rs   = 0.0;
1400:       jmin = bi[k]; jmax = bi[k+1]-1;
1401:       if (jmin < jmax) {
1402:         for (j=jmin; j<jmax; j++) {
1403:           col = bj[j]; ba[j] = rtmp[col]; rs += PetscAbsScalar(ba[j]);
1404:         }
1405:         /* add the k-th row into il and c2r */
1406:         il[k] = jmin;
1407:         i     = bj[jmin]; c2r[k] = c2r[i]; c2r[i] = k;
1408:       }

1410:       sctx.rs = rs;
1411:       sctx.pv = dk;
1412:       MatPivotCheck(B,A,info,&sctx,k);
1413:       if (sctx.newshift) break;
1414:       dk = sctx.pv;

1416:       ba[bdiag[k]] = 1.0/dk; /* U(k,k) */
1417:     }
1418:   } while (sctx.newshift);

1420:   PetscFree3(rtmp,il,c2r);

1422:   B->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1423:   B->ops->solves         = MatSolves_SeqSBAIJ_1;
1424:   B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1425:   B->ops->matsolve       = MatMatSolve_SeqSBAIJ_1_NaturalOrdering;
1426:   B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
1427:   B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;

1429:   B->assembled    = PETSC_TRUE;
1430:   B->preallocated = PETSC_TRUE;

1432:   PetscLogFlops(B->rmap->n);

1434:   /* MatPivotView() */
1435:   if (sctx.nshift) {
1436:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1437:       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);
1438:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1439:       PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
1440:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
1441:       PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %g\n",sctx.nshift,(double)info->shiftamount);
1442:     }
1443:   }
1444:   return(0);
1445: }

1447: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo *info)
1448: {
1449:   Mat_SeqSBAIJ   *a=(Mat_SeqSBAIJ*)A->data,*b=(Mat_SeqSBAIJ*)C->data;
1451:   PetscInt       i,j,mbs = a->mbs;
1452:   PetscInt       *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
1453:   PetscInt       k,jmin,*jl,*il,nexti,ili,*acol,*bcol,nz;
1454:   MatScalar      *rtmp,*ba=b->a,*aa=a->a,dk,uikdi,*aval,*bval;
1455:   PetscReal      rs;
1456:   FactorShiftCtx sctx;

1459:   /* MatPivotSetUp(): initialize shift context sctx */
1460:   PetscMemzero(&sctx,sizeof(FactorShiftCtx));

1462:   /* initialization */
1463:   /* il and jl record the first nonzero element in each row of the accessing
1464:      window U(0:k, k:mbs-1).
1465:      jl:    list of rows to be added to uneliminated rows
1466:             i>= k: jl(i) is the first row to be added to row i
1467:             i<  k: jl(i) is the row following row i in some list of rows
1468:             jl(i) = mbs indicates the end of a list
1469:      il(i): points to the first nonzero element in U(i,k:mbs-1)
1470:   */
1471:   PetscMalloc1(mbs,&rtmp);
1472:   PetscMalloc2(mbs,&il,mbs,&jl);

1474:   do {
1475:     sctx.newshift = PETSC_FALSE;
1476:     il[0] = 0;
1477:     for (i=0; i<mbs; i++) {
1478:       rtmp[i] = 0.0; jl[i] = mbs;
1479:     }

1481:     for (k = 0; k<mbs; k++) {
1482:       /*initialize k-th row with elements nonzero in row perm(k) of A */
1483:       nz   = ai[k+1] - ai[k];
1484:       acol = aj + ai[k];
1485:       aval = aa + ai[k];
1486:       bval = ba + bi[k];
1487:       while (nz--) {
1488:         rtmp[*acol++] = *aval++;
1489:         *bval++       = 0.0; /* for in-place factorization */
1490:       }

1492:       /* shift the diagonal of the matrix */
1493:       if (sctx.nshift) rtmp[k] += sctx.shift_amount;

1495:       /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1496:       dk = rtmp[k];
1497:       i  = jl[k]; /* first row to be added to k_th row  */

1499:       while (i < k) {
1500:         nexti = jl[i]; /* next row to be added to k_th row */
1501:         /* compute multiplier, update D(k) and U(i,k) */
1502:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
1503:         uikdi   = -ba[ili]*ba[bi[i]];
1504:         dk     += uikdi*ba[ili];
1505:         ba[ili] = uikdi; /* -U(i,k) */

1507:         /* add multiple of row i to k-th row ... */
1508:         jmin = ili + 1;
1509:         nz   = bi[i+1] - jmin;
1510:         if (nz > 0) {
1511:           bcol = bj + jmin;
1512:           bval = ba + jmin;
1513:           PetscLogFlops(2.0*nz);
1514:           while (nz--) rtmp[*bcol++] += uikdi*(*bval++);

1516:           /* update il and jl for i-th row */
1517:           il[i] = jmin;
1518:           j     = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
1519:         }
1520:         i = nexti;
1521:       }

1523:       /* shift the diagonals when zero pivot is detected */
1524:       /* compute rs=sum of abs(off-diagonal) */
1525:       rs   = 0.0;
1526:       jmin = bi[k]+1;
1527:       nz   = bi[k+1] - jmin;
1528:       if (nz) {
1529:         bcol = bj + jmin;
1530:         while (nz--) {
1531:           rs += PetscAbsScalar(rtmp[*bcol]);
1532:           bcol++;
1533:         }
1534:       }

1536:       sctx.rs = rs;
1537:       sctx.pv = dk;
1538:       MatPivotCheck(C,A,info,&sctx,k);
1539:       if (sctx.newshift) break;    /* sctx.shift_amount is updated */
1540:       dk = sctx.pv;

1542:       /* copy data into U(k,:) */
1543:       ba[bi[k]] = 1.0/dk;
1544:       jmin      = bi[k]+1;
1545:       nz        = bi[k+1] - jmin;
1546:       if (nz) {
1547:         bcol = bj + jmin;
1548:         bval = ba + jmin;
1549:         while (nz--) {
1550:           *bval++       = rtmp[*bcol];
1551:           rtmp[*bcol++] = 0.0;
1552:         }
1553:         /* add k-th row into il and jl */
1554:         il[k] = jmin;
1555:         i     = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
1556:       }
1557:     } /* end of for (k = 0; k<mbs; k++) */
1558:   } while (sctx.newshift);
1559:   PetscFree(rtmp);
1560:   PetscFree2(il,jl);

1562:   C->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1563:   C->ops->solves         = MatSolves_SeqSBAIJ_1_inplace;
1564:   C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1565:   C->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1566:   C->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;

1568:   C->assembled    = PETSC_TRUE;
1569:   C->preallocated = PETSC_TRUE;

1571:   PetscLogFlops(C->rmap->N);
1572:   if (sctx.nshift) {
1573:     if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1574:       PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
1575:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
1576:       PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
1577:     }
1578:   }
1579:   return(0);
1580: }

1582: PetscErrorCode MatCholeskyFactor_SeqSBAIJ(Mat A,IS perm,const MatFactorInfo *info)
1583: {
1585:   Mat            C;

1588:   MatGetFactor(A,"petsc",MAT_FACTOR_CHOLESKY,&C);
1589:   MatCholeskyFactorSymbolic(C,A,perm,info);
1590:   MatCholeskyFactorNumeric(C,A,info);

1592:   A->ops->solve          = C->ops->solve;
1593:   A->ops->solvetranspose = C->ops->solvetranspose;

1595:   MatHeaderMerge(A,&C);
1596:   return(0);
1597: }