Actual source code: sbaijfact.c

petsc-3.11.2 2019-05-18
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:       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,&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      = 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,&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 = 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: }