Actual source code: aijfact.c

petsc-master 2020-09-19
Report Typos and Errors

  2: #include <../src/mat/impls/aij/seq/aij.h>
  3: #include <../src/mat/impls/sbaij/seq/sbaij.h>
  4: #include <petscbt.h>
  5: #include <../src/mat/utils/freespace.h>

  7: /*
  8:       Computes an ordering to get most of the large numerical values in the lower triangular part of the matrix

 10:       This code does not work and is not called anywhere. It would be registered with MatOrderingRegisterAll()
 11: */
 12: PetscErrorCode MatGetOrdering_Flow_SeqAIJ(Mat mat,MatOrderingType type,IS *irow,IS *icol)
 13: {
 14:   Mat_SeqAIJ        *a = (Mat_SeqAIJ*)mat->data;
 15:   PetscErrorCode    ierr;
 16:   PetscInt          i,j,jj,k, kk,n = mat->rmap->n, current = 0, newcurrent = 0,*order;
 17:   const PetscInt    *ai = a->i, *aj = a->j;
 18:   const PetscScalar *aa = a->a;
 19:   PetscBool         *done;
 20:   PetscReal         best,past = 0,future;

 23:   /* pick initial row */
 24:   best = -1;
 25:   for (i=0; i<n; i++) {
 26:     future = 0.0;
 27:     for (j=ai[i]; j<ai[i+1]; j++) {
 28:       if (aj[j] != i) future += PetscAbsScalar(aa[j]);
 29:       else              past  = PetscAbsScalar(aa[j]);
 30:     }
 31:     if (!future) future = 1.e-10; /* if there is zero in the upper diagonal part want to rank this row high */
 32:     if (past/future > best) {
 33:       best    = past/future;
 34:       current = i;
 35:     }
 36:   }

 38:   PetscMalloc1(n,&done);
 39:   PetscArrayzero(done,n);
 40:   PetscMalloc1(n,&order);
 41:   order[0] = current;
 42:   for (i=0; i<n-1; i++) {
 43:     done[current] = PETSC_TRUE;
 44:     best          = -1;
 45:     /* loop over all neighbors of current pivot */
 46:     for (j=ai[current]; j<ai[current+1]; j++) {
 47:       jj = aj[j];
 48:       if (done[jj]) continue;
 49:       /* loop over columns of potential next row computing weights for below and above diagonal */
 50:       past = future = 0.0;
 51:       for (k=ai[jj]; k<ai[jj+1]; k++) {
 52:         kk = aj[k];
 53:         if (done[kk]) past += PetscAbsScalar(aa[k]);
 54:         else if (kk != jj) future += PetscAbsScalar(aa[k]);
 55:       }
 56:       if (!future) future = 1.e-10; /* if there is zero in the upper diagonal part want to rank this row high */
 57:       if (past/future > best) {
 58:         best       = past/future;
 59:         newcurrent = jj;
 60:       }
 61:     }
 62:     if (best == -1) { /* no neighbors to select from so select best of all that remain */
 63:       best = -1;
 64:       for (k=0; k<n; k++) {
 65:         if (done[k]) continue;
 66:         future = 0.0;
 67:         past   = 0.0;
 68:         for (j=ai[k]; j<ai[k+1]; j++) {
 69:           kk = aj[j];
 70:           if (done[kk])       past += PetscAbsScalar(aa[j]);
 71:           else if (kk != k) future += PetscAbsScalar(aa[j]);
 72:         }
 73:         if (!future) future = 1.e-10; /* if there is zero in the upper diagonal part want to rank this row high */
 74:         if (past/future > best) {
 75:           best       = past/future;
 76:           newcurrent = k;
 77:         }
 78:       }
 79:     }
 80:     if (current == newcurrent) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"newcurrent cannot be current");
 81:     current    = newcurrent;
 82:     order[i+1] = current;
 83:   }
 84:   ISCreateGeneral(PETSC_COMM_SELF,n,order,PETSC_COPY_VALUES,irow);
 85:   *icol = *irow;
 86:   PetscObjectReference((PetscObject)*irow);
 87:   PetscFree(done);
 88:   PetscFree(order);
 89:   return(0);
 90: }

 92: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_petsc(Mat A,MatFactorType ftype,Mat *B)
 93: {
 94:   PetscInt       n = A->rmap->n;

 98: #if defined(PETSC_USE_COMPLEX)
 99:   if (A->hermitian && !A->symmetric && (ftype == MAT_FACTOR_CHOLESKY||ftype == MAT_FACTOR_ICC)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian CHOLESKY or ICC Factor is not supported");
100: #endif
101:   MatCreate(PetscObjectComm((PetscObject)A),B);
102:   MatSetSizes(*B,n,n,n,n);
103:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU || ftype == MAT_FACTOR_ILUDT) {
104:     MatSetType(*B,MATSEQAIJ);

106:     (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqAIJ;
107:     (*B)->ops->lufactorsymbolic  = MatLUFactorSymbolic_SeqAIJ;

109:     MatSetBlockSizesFromMats(*B,A,A);
110:   } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
111:     MatSetType(*B,MATSEQSBAIJ);
112:     MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,NULL);

114:     (*B)->ops->iccfactorsymbolic      = MatICCFactorSymbolic_SeqAIJ;
115:     (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqAIJ;
116:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
117:   (*B)->factortype = ftype;

119:   PetscFree((*B)->solvertype);
120:   PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);
121:   (*B)->useordering = PETSC_TRUE;
122:   return(0);
123: }

125: PetscErrorCode MatLUFactorSymbolic_SeqAIJ_inplace(Mat B,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
126: {
127:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data,*b;
128:   IS                 isicol;
129:   PetscErrorCode     ierr;
130:   const PetscInt     *r,*ic;
131:   PetscInt           i,n=A->rmap->n,*ai=a->i,*aj=a->j;
132:   PetscInt           *bi,*bj,*ajtmp;
133:   PetscInt           *bdiag,row,nnz,nzi,reallocs=0,nzbd,*im;
134:   PetscReal          f;
135:   PetscInt           nlnk,*lnk,k,**bi_ptr;
136:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
137:   PetscBT            lnkbt;
138:   PetscBool          missing;

141:   if (A->rmap->N != A->cmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"matrix must be square");
142:   MatMissingDiagonal(A,&missing,&i);
143:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);

145:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
146:   ISGetIndices(isrow,&r);
147:   ISGetIndices(isicol,&ic);

149:   /* get new row pointers */
150:   PetscMalloc1(n+1,&bi);
151:   bi[0] = 0;

153:   /* bdiag is location of diagonal in factor */
154:   PetscMalloc1(n+1,&bdiag);
155:   bdiag[0] = 0;

157:   /* linked list for storing column indices of the active row */
158:   nlnk = n + 1;
159:   PetscLLCreate(n,n,nlnk,lnk,lnkbt);

161:   PetscMalloc2(n+1,&bi_ptr,n+1,&im);

163:   /* initial FreeSpace size is f*(ai[n]+1) */
164:   f             = info->fill;
165:   if (n==1)   f = 1; /* prevent failure in corner case of 1x1 matrix with fill < 0.5 */
166:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
167:   current_space = free_space;

169:   for (i=0; i<n; i++) {
170:     /* copy previous fill into linked list */
171:     nzi = 0;
172:     nnz = ai[r[i]+1] - ai[r[i]];
173:     ajtmp = aj + ai[r[i]];
174:     PetscLLAddPerm(nnz,ajtmp,ic,n,nlnk,lnk,lnkbt);
175:     nzi  += nlnk;

177:     /* add pivot rows into linked list */
178:     row = lnk[n];
179:     while (row < i) {
180:       nzbd  = bdiag[row] - bi[row] + 1;   /* num of entries in the row with column index <= row */
181:       ajtmp = bi_ptr[row] + nzbd;   /* points to the entry next to the diagonal */
182:       PetscLLAddSortedLU(ajtmp,row,nlnk,lnk,lnkbt,i,nzbd,im);
183:       nzi  += nlnk;
184:       row   = lnk[row];
185:     }
186:     bi[i+1] = bi[i] + nzi;
187:     im[i]   = nzi;

189:     /* mark bdiag */
190:     nzbd = 0;
191:     nnz  = nzi;
192:     k    = lnk[n];
193:     while (nnz-- && k < i) {
194:       nzbd++;
195:       k = lnk[k];
196:     }
197:     bdiag[i] = bi[i] + nzbd;

199:     /* if free space is not available, make more free space */
200:     if (current_space->local_remaining<nzi) {
201:       nnz  = PetscIntMultTruncate(n - i,nzi); /* estimated and max additional space needed */
202:       PetscFreeSpaceGet(nnz,&current_space);
203:       reallocs++;
204:     }

206:     /* copy data into free space, then initialize lnk */
207:     PetscLLClean(n,n,nzi,lnk,current_space->array,lnkbt);

209:     bi_ptr[i]                       = current_space->array;
210:     current_space->array           += nzi;
211:     current_space->local_used      += nzi;
212:     current_space->local_remaining -= nzi;
213:   }
214: #if defined(PETSC_USE_INFO)
215:   if (ai[n] != 0) {
216:     PetscReal af = ((PetscReal)bi[n])/((PetscReal)ai[n]);
217:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);
218:     PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
219:     PetscInfo1(A,"PCFactorSetFill(pc,%g);\n",(double)af);
220:     PetscInfo(A,"for best performance.\n");
221:   } else {
222:     PetscInfo(A,"Empty matrix\n");
223:   }
224: #endif

226:   ISRestoreIndices(isrow,&r);
227:   ISRestoreIndices(isicol,&ic);

229:   /* destroy list of free space and other temporary array(s) */
230:   PetscMalloc1(bi[n]+1,&bj);
231:   PetscFreeSpaceContiguous(&free_space,bj);
232:   PetscLLDestroy(lnk,lnkbt);
233:   PetscFree2(bi_ptr,im);

235:   /* put together the new matrix */
236:   MatSeqAIJSetPreallocation_SeqAIJ(B,MAT_SKIP_ALLOCATION,NULL);
237:   PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);
238:   b    = (Mat_SeqAIJ*)(B)->data;

240:   b->free_a       = PETSC_TRUE;
241:   b->free_ij      = PETSC_TRUE;
242:   b->singlemalloc = PETSC_FALSE;

244:   PetscMalloc1(bi[n]+1,&b->a);
245:   b->j    = bj;
246:   b->i    = bi;
247:   b->diag = bdiag;
248:   b->ilen = NULL;
249:   b->imax = NULL;
250:   b->row  = isrow;
251:   b->col  = iscol;
252:   PetscObjectReference((PetscObject)isrow);
253:   PetscObjectReference((PetscObject)iscol);
254:   b->icol = isicol;
255:   PetscMalloc1(n+1,&b->solve_work);

257:   /* In b structure:  Free imax, ilen, old a, old j.  Allocate solve_work, new a, new j */
258:   PetscLogObjectMemory((PetscObject)B,(bi[n]-n)*(sizeof(PetscInt)+sizeof(PetscScalar)));
259:   b->maxnz = b->nz = bi[n];

261:   (B)->factortype            = MAT_FACTOR_LU;
262:   (B)->info.factor_mallocs   = reallocs;
263:   (B)->info.fill_ratio_given = f;

265:   if (ai[n]) {
266:     (B)->info.fill_ratio_needed = ((PetscReal)bi[n])/((PetscReal)ai[n]);
267:   } else {
268:     (B)->info.fill_ratio_needed = 0.0;
269:   }
270:   (B)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_inplace;
271:   if (a->inode.size) {
272:     (B)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode_inplace;
273:   }
274:   return(0);
275: }

277: PetscErrorCode MatLUFactorSymbolic_SeqAIJ(Mat B,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
278: {
279:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data,*b;
280:   IS                 isicol;
281:   PetscErrorCode     ierr;
282:   const PetscInt     *r,*ic,*ai=a->i,*aj=a->j,*ajtmp;
283:   PetscInt           i,n=A->rmap->n;
284:   PetscInt           *bi,*bj;
285:   PetscInt           *bdiag,row,nnz,nzi,reallocs=0,nzbd,*im;
286:   PetscReal          f;
287:   PetscInt           nlnk,*lnk,k,**bi_ptr;
288:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
289:   PetscBT            lnkbt;
290:   PetscBool          missing;

293:   if (A->rmap->N != A->cmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"matrix must be square");
294:   MatMissingDiagonal(A,&missing,&i);
295:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);

297:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
298:   ISGetIndices(isrow,&r);
299:   ISGetIndices(isicol,&ic);

301:   /* get new row and diagonal pointers, must be allocated separately because they will be given to the Mat_SeqAIJ and freed separately */
302:   PetscMalloc1(n+1,&bi);
303:   PetscMalloc1(n+1,&bdiag);
304:   bi[0] = bdiag[0] = 0;

306:   /* linked list for storing column indices of the active row */
307:   nlnk = n + 1;
308:   PetscLLCreate(n,n,nlnk,lnk,lnkbt);

310:   PetscMalloc2(n+1,&bi_ptr,n+1,&im);

312:   /* initial FreeSpace size is f*(ai[n]+1) */
313:   f             = info->fill;
314:   if (n==1)   f = 1; /* prevent failure in corner case of 1x1 matrix with fill < 0.5 */
315:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
316:   current_space = free_space;

318:   for (i=0; i<n; i++) {
319:     /* copy previous fill into linked list */
320:     nzi = 0;
321:     nnz = ai[r[i]+1] - ai[r[i]];
322:     ajtmp = aj + ai[r[i]];
323:     PetscLLAddPerm(nnz,ajtmp,ic,n,nlnk,lnk,lnkbt);
324:     nzi  += nlnk;

326:     /* add pivot rows into linked list */
327:     row = lnk[n];
328:     while (row < i) {
329:       nzbd  = bdiag[row] + 1; /* num of entries in the row with column index <= row */
330:       ajtmp = bi_ptr[row] + nzbd; /* points to the entry next to the diagonal */
331:       PetscLLAddSortedLU(ajtmp,row,nlnk,lnk,lnkbt,i,nzbd,im);
332:       nzi  += nlnk;
333:       row   = lnk[row];
334:     }
335:     bi[i+1] = bi[i] + nzi;
336:     im[i]   = nzi;

338:     /* mark bdiag */
339:     nzbd = 0;
340:     nnz  = nzi;
341:     k    = lnk[n];
342:     while (nnz-- && k < i) {
343:       nzbd++;
344:       k = lnk[k];
345:     }
346:     bdiag[i] = nzbd; /* note: bdiag[i] = nnzL as input for PetscFreeSpaceContiguous_LU() */

348:     /* if free space is not available, make more free space */
349:     if (current_space->local_remaining<nzi) {
350:       /* estimated additional space needed */
351:       nnz  = PetscIntMultTruncate(2,PetscIntMultTruncate(n-1,nzi));
352:       PetscFreeSpaceGet(nnz,&current_space);
353:       reallocs++;
354:     }

356:     /* copy data into free space, then initialize lnk */
357:     PetscLLClean(n,n,nzi,lnk,current_space->array,lnkbt);

359:     bi_ptr[i]                       = current_space->array;
360:     current_space->array           += nzi;
361:     current_space->local_used      += nzi;
362:     current_space->local_remaining -= nzi;
363:   }

365:   ISRestoreIndices(isrow,&r);
366:   ISRestoreIndices(isicol,&ic);

368:   /*   copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
369:   PetscMalloc1(bi[n]+1,&bj);
370:   PetscFreeSpaceContiguous_LU(&free_space,bj,n,bi,bdiag);
371:   PetscLLDestroy(lnk,lnkbt);
372:   PetscFree2(bi_ptr,im);

374:   /* put together the new matrix */
375:   MatSeqAIJSetPreallocation_SeqAIJ(B,MAT_SKIP_ALLOCATION,NULL);
376:   PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);
377:   b    = (Mat_SeqAIJ*)(B)->data;

379:   b->free_a       = PETSC_TRUE;
380:   b->free_ij      = PETSC_TRUE;
381:   b->singlemalloc = PETSC_FALSE;

383:   PetscMalloc1(bdiag[0]+1,&b->a);

385:   b->j    = bj;
386:   b->i    = bi;
387:   b->diag = bdiag;
388:   b->ilen = NULL;
389:   b->imax = NULL;
390:   b->row  = isrow;
391:   b->col  = iscol;
392:   PetscObjectReference((PetscObject)isrow);
393:   PetscObjectReference((PetscObject)iscol);
394:   b->icol = isicol;
395:   PetscMalloc1(n+1,&b->solve_work);

397:   /* In b structure:  Free imax, ilen, old a, old j.  Allocate solve_work, new a, new j */
398:   PetscLogObjectMemory((PetscObject)B,(bdiag[0]+1)*(sizeof(PetscInt)+sizeof(PetscScalar)));
399:   b->maxnz = b->nz = bdiag[0]+1;

401:   B->factortype            = MAT_FACTOR_LU;
402:   B->info.factor_mallocs   = reallocs;
403:   B->info.fill_ratio_given = f;

405:   if (ai[n]) {
406:     B->info.fill_ratio_needed = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]);
407:   } else {
408:     B->info.fill_ratio_needed = 0.0;
409:   }
410: #if defined(PETSC_USE_INFO)
411:   if (ai[n] != 0) {
412:     PetscReal af = B->info.fill_ratio_needed;
413:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);
414:     PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
415:     PetscInfo1(A,"PCFactorSetFill(pc,%g);\n",(double)af);
416:     PetscInfo(A,"for best performance.\n");
417:   } else {
418:     PetscInfo(A,"Empty matrix\n");
419:   }
420: #endif
421:   B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ;
422:   if (a->inode.size) {
423:     B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
424:   }
425:   MatSeqAIJCheckInode_FactorLU(B);
426:   return(0);
427: }

429: /*
430:     Trouble in factorization, should we dump the original matrix?
431: */
432: PetscErrorCode MatFactorDumpMatrix(Mat A)
433: {
435:   PetscBool      flg = PETSC_FALSE;

438:   PetscOptionsGetBool(((PetscObject)A)->options,NULL,"-mat_factor_dump_on_error",&flg,NULL);
439:   if (flg) {
440:     PetscViewer viewer;
441:     char        filename[PETSC_MAX_PATH_LEN];

443:     PetscSNPrintf(filename,PETSC_MAX_PATH_LEN,"matrix_factor_error.%d",PetscGlobalRank);
444:     PetscViewerBinaryOpen(PetscObjectComm((PetscObject)A),filename,FILE_MODE_WRITE,&viewer);
445:     MatView(A,viewer);
446:     PetscViewerDestroy(&viewer);
447:   }
448:   return(0);
449: }

451: PetscErrorCode MatLUFactorNumeric_SeqAIJ(Mat B,Mat A,const MatFactorInfo *info)
452: {
453:   Mat             C     =B;
454:   Mat_SeqAIJ      *a    =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)C->data;
455:   IS              isrow = b->row,isicol = b->icol;
456:   PetscErrorCode  ierr;
457:   const PetscInt  *r,*ic,*ics;
458:   const PetscInt  n=A->rmap->n,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bdiag=b->diag;
459:   PetscInt        i,j,k,nz,nzL,row,*pj;
460:   const PetscInt  *ajtmp,*bjtmp;
461:   MatScalar       *rtmp,*pc,multiplier,*pv;
462:   const MatScalar *aa=a->a,*v;
463:   PetscBool       row_identity,col_identity;
464:   FactorShiftCtx  sctx;
465:   const PetscInt  *ddiag;
466:   PetscReal       rs;
467:   MatScalar       d;

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

473:   if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
474:     ddiag          = a->diag;
475:     sctx.shift_top = info->zeropivot;
476:     for (i=0; i<n; i++) {
477:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
478:       d  = (aa)[ddiag[i]];
479:       rs = -PetscAbsScalar(d) - PetscRealPart(d);
480:       v  = aa+ai[i];
481:       nz = ai[i+1] - ai[i];
482:       for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
483:       if (rs>sctx.shift_top) sctx.shift_top = rs;
484:     }
485:     sctx.shift_top *= 1.1;
486:     sctx.nshift_max = 5;
487:     sctx.shift_lo   = 0.;
488:     sctx.shift_hi   = 1.;
489:   }

491:   ISGetIndices(isrow,&r);
492:   ISGetIndices(isicol,&ic);
493:   PetscMalloc1(n+1,&rtmp);
494:   ics  = ic;

496:   do {
497:     sctx.newshift = PETSC_FALSE;
498:     for (i=0; i<n; i++) {
499:       /* zero rtmp */
500:       /* L part */
501:       nz    = bi[i+1] - bi[i];
502:       bjtmp = bj + bi[i];
503:       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

505:       /* U part */
506:       nz    = bdiag[i]-bdiag[i+1];
507:       bjtmp = bj + bdiag[i+1]+1;
508:       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

510:       /* load in initial (unfactored row) */
511:       nz    = ai[r[i]+1] - ai[r[i]];
512:       ajtmp = aj + ai[r[i]];
513:       v     = aa + ai[r[i]];
514:       for (j=0; j<nz; j++) {
515:         rtmp[ics[ajtmp[j]]] = v[j];
516:       }
517:       /* ZeropivotApply() */
518:       rtmp[i] += sctx.shift_amount;  /* shift the diagonal of the matrix */

520:       /* elimination */
521:       bjtmp = bj + bi[i];
522:       row   = *bjtmp++;
523:       nzL   = bi[i+1] - bi[i];
524:       for (k=0; k < nzL; k++) {
525:         pc = rtmp + row;
526:         if (*pc != 0.0) {
527:           pv         = b->a + bdiag[row];
528:           multiplier = *pc * (*pv);
529:           *pc        = multiplier;

531:           pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
532:           pv = b->a + bdiag[row+1]+1;
533:           nz = bdiag[row]-bdiag[row+1]-1; /* num of entries in U(row,:) excluding diag */

535:           for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
536:           PetscLogFlops(1+2.0*nz);
537:         }
538:         row = *bjtmp++;
539:       }

541:       /* finished row so stick it into b->a */
542:       rs = 0.0;
543:       /* L part */
544:       pv = b->a + bi[i];
545:       pj = b->j + bi[i];
546:       nz = bi[i+1] - bi[i];
547:       for (j=0; j<nz; j++) {
548:         pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
549:       }

551:       /* U part */
552:       pv = b->a + bdiag[i+1]+1;
553:       pj = b->j + bdiag[i+1]+1;
554:       nz = bdiag[i] - bdiag[i+1]-1;
555:       for (j=0; j<nz; j++) {
556:         pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
557:       }

559:       sctx.rs = rs;
560:       sctx.pv = rtmp[i];
561:       MatPivotCheck(B,A,info,&sctx,i);
562:       if (sctx.newshift) break; /* break for-loop */
563:       rtmp[i] = sctx.pv; /* sctx.pv might be updated in the case of MAT_SHIFT_INBLOCKS */

565:       /* Mark diagonal and invert diagonal for simplier triangular solves */
566:       pv  = b->a + bdiag[i];
567:       *pv = 1.0/rtmp[i];

569:     } /* endof for (i=0; i<n; i++) { */

571:     /* MatPivotRefine() */
572:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
573:       /*
574:        * if no shift in this attempt & shifting & started shifting & can refine,
575:        * then try lower shift
576:        */
577:       sctx.shift_hi       = sctx.shift_fraction;
578:       sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
579:       sctx.shift_amount   = sctx.shift_fraction * sctx.shift_top;
580:       sctx.newshift       = PETSC_TRUE;
581:       sctx.nshift++;
582:     }
583:   } while (sctx.newshift);

585:   PetscFree(rtmp);
586:   ISRestoreIndices(isicol,&ic);
587:   ISRestoreIndices(isrow,&r);

589:   ISIdentity(isrow,&row_identity);
590:   ISIdentity(isicol,&col_identity);
591:   if (b->inode.size) {
592:     C->ops->solve = MatSolve_SeqAIJ_Inode;
593:   } else if (row_identity && col_identity) {
594:     C->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
595:   } else {
596:     C->ops->solve = MatSolve_SeqAIJ;
597:   }
598:   C->ops->solveadd          = MatSolveAdd_SeqAIJ;
599:   C->ops->solvetranspose    = MatSolveTranspose_SeqAIJ;
600:   C->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ;
601:   C->ops->matsolve          = MatMatSolve_SeqAIJ;
602:   C->assembled              = PETSC_TRUE;
603:   C->preallocated           = PETSC_TRUE;

605:   PetscLogFlops(C->cmap->n);

607:   /* MatShiftView(A,info,&sctx) */
608:   if (sctx.nshift) {
609:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
610:       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);
611:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
612:       PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
613:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
614:       PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %g\n",sctx.nshift,(double)info->shiftamount);
615:     }
616:   }
617:   return(0);
618: }

620: PetscErrorCode MatLUFactorNumeric_SeqAIJ_inplace(Mat B,Mat A,const MatFactorInfo *info)
621: {
622:   Mat             C     =B;
623:   Mat_SeqAIJ      *a    =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)C->data;
624:   IS              isrow = b->row,isicol = b->icol;
625:   PetscErrorCode  ierr;
626:   const PetscInt  *r,*ic,*ics;
627:   PetscInt        nz,row,i,j,n=A->rmap->n,diag;
628:   const PetscInt  *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
629:   const PetscInt  *ajtmp,*bjtmp,*diag_offset = b->diag,*pj;
630:   MatScalar       *pv,*rtmp,*pc,multiplier,d;
631:   const MatScalar *v,*aa=a->a;
632:   PetscReal       rs=0.0;
633:   FactorShiftCtx  sctx;
634:   const PetscInt  *ddiag;
635:   PetscBool       row_identity, col_identity;

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

641:   if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
642:     ddiag          = a->diag;
643:     sctx.shift_top = info->zeropivot;
644:     for (i=0; i<n; i++) {
645:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
646:       d  = (aa)[ddiag[i]];
647:       rs = -PetscAbsScalar(d) - PetscRealPart(d);
648:       v  = aa+ai[i];
649:       nz = ai[i+1] - ai[i];
650:       for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
651:       if (rs>sctx.shift_top) sctx.shift_top = rs;
652:     }
653:     sctx.shift_top *= 1.1;
654:     sctx.nshift_max = 5;
655:     sctx.shift_lo   = 0.;
656:     sctx.shift_hi   = 1.;
657:   }

659:   ISGetIndices(isrow,&r);
660:   ISGetIndices(isicol,&ic);
661:   PetscMalloc1(n+1,&rtmp);
662:   ics  = ic;

664:   do {
665:     sctx.newshift = PETSC_FALSE;
666:     for (i=0; i<n; i++) {
667:       nz    = bi[i+1] - bi[i];
668:       bjtmp = bj + bi[i];
669:       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

671:       /* load in initial (unfactored row) */
672:       nz    = ai[r[i]+1] - ai[r[i]];
673:       ajtmp = aj + ai[r[i]];
674:       v     = aa + ai[r[i]];
675:       for (j=0; j<nz; j++) {
676:         rtmp[ics[ajtmp[j]]] = v[j];
677:       }
678:       rtmp[ics[r[i]]] += sctx.shift_amount; /* shift the diagonal of the matrix */

680:       row = *bjtmp++;
681:       while  (row < i) {
682:         pc = rtmp + row;
683:         if (*pc != 0.0) {
684:           pv         = b->a + diag_offset[row];
685:           pj         = b->j + diag_offset[row] + 1;
686:           multiplier = *pc / *pv++;
687:           *pc        = multiplier;
688:           nz         = bi[row+1] - diag_offset[row] - 1;
689:           for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
690:           PetscLogFlops(1+2.0*nz);
691:         }
692:         row = *bjtmp++;
693:       }
694:       /* finished row so stick it into b->a */
695:       pv   = b->a + bi[i];
696:       pj   = b->j + bi[i];
697:       nz   = bi[i+1] - bi[i];
698:       diag = diag_offset[i] - bi[i];
699:       rs   = 0.0;
700:       for (j=0; j<nz; j++) {
701:         pv[j] = rtmp[pj[j]];
702:         rs   += PetscAbsScalar(pv[j]);
703:       }
704:       rs -= PetscAbsScalar(pv[diag]);

706:       sctx.rs = rs;
707:       sctx.pv = pv[diag];
708:       MatPivotCheck(B,A,info,&sctx,i);
709:       if (sctx.newshift) break;
710:       pv[diag] = sctx.pv;
711:     }

713:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
714:       /*
715:        * if no shift in this attempt & shifting & started shifting & can refine,
716:        * then try lower shift
717:        */
718:       sctx.shift_hi       = sctx.shift_fraction;
719:       sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
720:       sctx.shift_amount   = sctx.shift_fraction * sctx.shift_top;
721:       sctx.newshift       = PETSC_TRUE;
722:       sctx.nshift++;
723:     }
724:   } while (sctx.newshift);

726:   /* invert diagonal entries for simplier triangular solves */
727:   for (i=0; i<n; i++) {
728:     b->a[diag_offset[i]] = 1.0/b->a[diag_offset[i]];
729:   }
730:   PetscFree(rtmp);
731:   ISRestoreIndices(isicol,&ic);
732:   ISRestoreIndices(isrow,&r);

734:   ISIdentity(isrow,&row_identity);
735:   ISIdentity(isicol,&col_identity);
736:   if (row_identity && col_identity) {
737:     C->ops->solve = MatSolve_SeqAIJ_NaturalOrdering_inplace;
738:   } else {
739:     C->ops->solve = MatSolve_SeqAIJ_inplace;
740:   }
741:   C->ops->solveadd          = MatSolveAdd_SeqAIJ_inplace;
742:   C->ops->solvetranspose    = MatSolveTranspose_SeqAIJ_inplace;
743:   C->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ_inplace;
744:   C->ops->matsolve          = MatMatSolve_SeqAIJ_inplace;

746:   C->assembled    = PETSC_TRUE;
747:   C->preallocated = PETSC_TRUE;

749:   PetscLogFlops(C->cmap->n);
750:   if (sctx.nshift) {
751:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
752:       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);
753:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
754:       PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
755:     }
756:   }
757:   (C)->ops->solve          = MatSolve_SeqAIJ_inplace;
758:   (C)->ops->solvetranspose = MatSolveTranspose_SeqAIJ_inplace;

760:   MatSeqAIJCheckInode(C);
761:   return(0);
762: }

764: /*
765:    This routine implements inplace ILU(0) with row or/and column permutations.
766:    Input:
767:      A - original matrix
768:    Output;
769:      A - a->i (rowptr) is same as original rowptr, but factored i-the row is stored in rowperm[i]
770:          a->j (col index) is permuted by the inverse of colperm, then sorted
771:          a->a reordered accordingly with a->j
772:          a->diag (ptr to diagonal elements) is updated.
773: */
774: PetscErrorCode MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(Mat B,Mat A,const MatFactorInfo *info)
775: {
776:   Mat_SeqAIJ      *a    =(Mat_SeqAIJ*)A->data;
777:   IS              isrow = a->row,isicol = a->icol;
778:   PetscErrorCode  ierr;
779:   const PetscInt  *r,*ic,*ics;
780:   PetscInt        i,j,n=A->rmap->n,*ai=a->i,*aj=a->j;
781:   PetscInt        *ajtmp,nz,row;
782:   PetscInt        *diag = a->diag,nbdiag,*pj;
783:   PetscScalar     *rtmp,*pc,multiplier,d;
784:   MatScalar       *pv,*v;
785:   PetscReal       rs;
786:   FactorShiftCtx  sctx;
787:   const MatScalar *aa=a->a,*vtmp;

790:   if (A != B) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"input and output matrix must have same address");

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

795:   if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
796:     const PetscInt *ddiag = a->diag;
797:     sctx.shift_top = info->zeropivot;
798:     for (i=0; i<n; i++) {
799:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
800:       d    = (aa)[ddiag[i]];
801:       rs   = -PetscAbsScalar(d) - PetscRealPart(d);
802:       vtmp = aa+ai[i];
803:       nz   = ai[i+1] - ai[i];
804:       for (j=0; j<nz; j++) rs += PetscAbsScalar(vtmp[j]);
805:       if (rs>sctx.shift_top) sctx.shift_top = rs;
806:     }
807:     sctx.shift_top *= 1.1;
808:     sctx.nshift_max = 5;
809:     sctx.shift_lo   = 0.;
810:     sctx.shift_hi   = 1.;
811:   }

813:   ISGetIndices(isrow,&r);
814:   ISGetIndices(isicol,&ic);
815:   PetscMalloc1(n+1,&rtmp);
816:   PetscArrayzero(rtmp,n+1);
817:   ics  = ic;

819: #if defined(MV)
820:   sctx.shift_top      = 0.;
821:   sctx.nshift_max     = 0;
822:   sctx.shift_lo       = 0.;
823:   sctx.shift_hi       = 0.;
824:   sctx.shift_fraction = 0.;

826:   if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
827:     sctx.shift_top = 0.;
828:     for (i=0; i<n; i++) {
829:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
830:       d  = (a->a)[diag[i]];
831:       rs = -PetscAbsScalar(d) - PetscRealPart(d);
832:       v  = a->a+ai[i];
833:       nz = ai[i+1] - ai[i];
834:       for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
835:       if (rs>sctx.shift_top) sctx.shift_top = rs;
836:     }
837:     if (sctx.shift_top < info->zeropivot) sctx.shift_top = info->zeropivot;
838:     sctx.shift_top *= 1.1;
839:     sctx.nshift_max = 5;
840:     sctx.shift_lo   = 0.;
841:     sctx.shift_hi   = 1.;
842:   }

844:   sctx.shift_amount = 0.;
845:   sctx.nshift       = 0;
846: #endif

848:   do {
849:     sctx.newshift = PETSC_FALSE;
850:     for (i=0; i<n; i++) {
851:       /* load in initial unfactored row */
852:       nz    = ai[r[i]+1] - ai[r[i]];
853:       ajtmp = aj + ai[r[i]];
854:       v     = a->a + ai[r[i]];
855:       /* sort permuted ajtmp and values v accordingly */
856:       for (j=0; j<nz; j++) ajtmp[j] = ics[ajtmp[j]];
857:       PetscSortIntWithScalarArray(nz,ajtmp,v);

859:       diag[r[i]] = ai[r[i]];
860:       for (j=0; j<nz; j++) {
861:         rtmp[ajtmp[j]] = v[j];
862:         if (ajtmp[j] < i) diag[r[i]]++; /* update a->diag */
863:       }
864:       rtmp[r[i]] += sctx.shift_amount; /* shift the diagonal of the matrix */

866:       row = *ajtmp++;
867:       while  (row < i) {
868:         pc = rtmp + row;
869:         if (*pc != 0.0) {
870:           pv = a->a + diag[r[row]];
871:           pj = aj + diag[r[row]] + 1;

873:           multiplier = *pc / *pv++;
874:           *pc        = multiplier;
875:           nz         = ai[r[row]+1] - diag[r[row]] - 1;
876:           for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
877:           PetscLogFlops(1+2.0*nz);
878:         }
879:         row = *ajtmp++;
880:       }
881:       /* finished row so overwrite it onto a->a */
882:       pv     = a->a + ai[r[i]];
883:       pj     = aj + ai[r[i]];
884:       nz     = ai[r[i]+1] - ai[r[i]];
885:       nbdiag = diag[r[i]] - ai[r[i]]; /* num of entries before the diagonal */

887:       rs = 0.0;
888:       for (j=0; j<nz; j++) {
889:         pv[j] = rtmp[pj[j]];
890:         if (j != nbdiag) rs += PetscAbsScalar(pv[j]);
891:       }

893:       sctx.rs = rs;
894:       sctx.pv = pv[nbdiag];
895:       MatPivotCheck(B,A,info,&sctx,i);
896:       if (sctx.newshift) break;
897:       pv[nbdiag] = sctx.pv;
898:     }

900:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
901:       /*
902:        * if no shift in this attempt & shifting & started shifting & can refine,
903:        * then try lower shift
904:        */
905:       sctx.shift_hi       = sctx.shift_fraction;
906:       sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
907:       sctx.shift_amount   = sctx.shift_fraction * sctx.shift_top;
908:       sctx.newshift       = PETSC_TRUE;
909:       sctx.nshift++;
910:     }
911:   } while (sctx.newshift);

913:   /* invert diagonal entries for simplier triangular solves */
914:   for (i=0; i<n; i++) {
915:     a->a[diag[r[i]]] = 1.0/a->a[diag[r[i]]];
916:   }

918:   PetscFree(rtmp);
919:   ISRestoreIndices(isicol,&ic);
920:   ISRestoreIndices(isrow,&r);

922:   A->ops->solve             = MatSolve_SeqAIJ_InplaceWithPerm;
923:   A->ops->solveadd          = MatSolveAdd_SeqAIJ_inplace;
924:   A->ops->solvetranspose    = MatSolveTranspose_SeqAIJ_inplace;
925:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ_inplace;

927:   A->assembled    = PETSC_TRUE;
928:   A->preallocated = PETSC_TRUE;

930:   PetscLogFlops(A->cmap->n);
931:   if (sctx.nshift) {
932:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
933:       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);
934:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
935:       PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
936:     }
937:   }
938:   return(0);
939: }

941: /* ----------------------------------------------------------- */
942: PetscErrorCode MatLUFactor_SeqAIJ(Mat A,IS row,IS col,const MatFactorInfo *info)
943: {
945:   Mat            C;

948:   MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&C);
949:   MatLUFactorSymbolic(C,A,row,col,info);
950:   MatLUFactorNumeric(C,A,info);

952:   A->ops->solve          = C->ops->solve;
953:   A->ops->solvetranspose = C->ops->solvetranspose;

955:   MatHeaderMerge(A,&C);
956:   PetscLogObjectParent((PetscObject)A,(PetscObject)((Mat_SeqAIJ*)(A->data))->icol);
957:   return(0);
958: }
959: /* ----------------------------------------------------------- */


962: PetscErrorCode MatSolve_SeqAIJ_inplace(Mat A,Vec bb,Vec xx)
963: {
964:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
965:   IS                iscol = a->col,isrow = a->row;
966:   PetscErrorCode    ierr;
967:   PetscInt          i, n = A->rmap->n,*vi,*ai = a->i,*aj = a->j;
968:   PetscInt          nz;
969:   const PetscInt    *rout,*cout,*r,*c;
970:   PetscScalar       *x,*tmp,*tmps,sum;
971:   const PetscScalar *b;
972:   const MatScalar   *aa = a->a,*v;

975:   if (!n) return(0);

977:   VecGetArrayRead(bb,&b);
978:   VecGetArrayWrite(xx,&x);
979:   tmp  = a->solve_work;

981:   ISGetIndices(isrow,&rout); r = rout;
982:   ISGetIndices(iscol,&cout); c = cout + (n-1);

984:   /* forward solve the lower triangular */
985:   tmp[0] = b[*r++];
986:   tmps   = tmp;
987:   for (i=1; i<n; i++) {
988:     v   = aa + ai[i];
989:     vi  = aj + ai[i];
990:     nz  = a->diag[i] - ai[i];
991:     sum = b[*r++];
992:     PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
993:     tmp[i] = sum;
994:   }

996:   /* backward solve the upper triangular */
997:   for (i=n-1; i>=0; i--) {
998:     v   = aa + a->diag[i] + 1;
999:     vi  = aj + a->diag[i] + 1;
1000:     nz  = ai[i+1] - a->diag[i] - 1;
1001:     sum = tmp[i];
1002:     PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
1003:     x[*c--] = tmp[i] = sum*aa[a->diag[i]];
1004:   }

1006:   ISRestoreIndices(isrow,&rout);
1007:   ISRestoreIndices(iscol,&cout);
1008:   VecRestoreArrayRead(bb,&b);
1009:   VecRestoreArrayWrite(xx,&x);
1010:   PetscLogFlops(2.0*a->nz - A->cmap->n);
1011:   return(0);
1012: }

1014: PetscErrorCode MatMatSolve_SeqAIJ_inplace(Mat A,Mat B,Mat X)
1015: {
1016:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1017:   IS                iscol = a->col,isrow = a->row;
1018:   PetscErrorCode    ierr;
1019:   PetscInt          i, n = A->rmap->n,*vi,*ai = a->i,*aj = a->j;
1020:   PetscInt          nz,neq,ldb,ldx;
1021:   const PetscInt    *rout,*cout,*r,*c;
1022:   PetscScalar       *x,*tmp = a->solve_work,*tmps,sum;
1023:   const PetscScalar *b,*aa = a->a,*v;
1024:   PetscBool         isdense;

1027:   if (!n) return(0);
1028:   PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&isdense);
1029:   if (!isdense) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"B matrix must be a SeqDense matrix");
1030:   if (X != B) {
1031:     PetscObjectTypeCompare((PetscObject)X,MATSEQDENSE,&isdense);
1032:     if (!isdense) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"X matrix must be a SeqDense matrix");
1033:   }
1034:   MatDenseGetArrayRead(B,&b);
1035:   MatDenseGetLDA(B,&ldb);
1036:   MatDenseGetArray(X,&x);
1037:   MatDenseGetLDA(X,&ldx);
1038:   ISGetIndices(isrow,&rout); r = rout;
1039:   ISGetIndices(iscol,&cout); c = cout;
1040:   for (neq=0; neq<B->cmap->n; neq++) {
1041:     /* forward solve the lower triangular */
1042:     tmp[0] = b[r[0]];
1043:     tmps   = tmp;
1044:     for (i=1; i<n; i++) {
1045:       v   = aa + ai[i];
1046:       vi  = aj + ai[i];
1047:       nz  = a->diag[i] - ai[i];
1048:       sum = b[r[i]];
1049:       PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
1050:       tmp[i] = sum;
1051:     }
1052:     /* backward solve the upper triangular */
1053:     for (i=n-1; i>=0; i--) {
1054:       v   = aa + a->diag[i] + 1;
1055:       vi  = aj + a->diag[i] + 1;
1056:       nz  = ai[i+1] - a->diag[i] - 1;
1057:       sum = tmp[i];
1058:       PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
1059:       x[c[i]] = tmp[i] = sum*aa[a->diag[i]];
1060:     }
1061:     b += ldb;
1062:     x += ldx;
1063:   }
1064:   ISRestoreIndices(isrow,&rout);
1065:   ISRestoreIndices(iscol,&cout);
1066:   MatDenseRestoreArrayRead(B,&b);
1067:   MatDenseRestoreArray(X,&x);
1068:   PetscLogFlops(B->cmap->n*(2.0*a->nz - n));
1069:   return(0);
1070: }

1072: PetscErrorCode MatMatSolve_SeqAIJ(Mat A,Mat B,Mat X)
1073: {
1074:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1075:   IS                iscol = a->col,isrow = a->row;
1076:   PetscErrorCode    ierr;
1077:   PetscInt          i, n = A->rmap->n,*vi,*ai = a->i,*aj = a->j,*adiag = a->diag;
1078:   PetscInt          nz,neq,ldb,ldx;
1079:   const PetscInt    *rout,*cout,*r,*c;
1080:   PetscScalar       *x,*tmp = a->solve_work,sum;
1081:   const PetscScalar *b,*aa = a->a,*v;
1082:   PetscBool         isdense;

1085:   if (!n) return(0);
1086:   PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&isdense);
1087:   if (!isdense) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"B matrix must be a SeqDense matrix");
1088:   if (X != B) {
1089:     PetscObjectTypeCompare((PetscObject)X,MATSEQDENSE,&isdense);
1090:     if (!isdense) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"X matrix must be a SeqDense matrix");
1091:   }
1092:   MatDenseGetArrayRead(B,&b);
1093:   MatDenseGetLDA(B,&ldb);
1094:   MatDenseGetArray(X,&x);
1095:   MatDenseGetLDA(X,&ldx);
1096:   ISGetIndices(isrow,&rout); r = rout;
1097:   ISGetIndices(iscol,&cout); c = cout;
1098:   for (neq=0; neq<B->cmap->n; neq++) {
1099:     /* forward solve the lower triangular */
1100:     tmp[0] = b[r[0]];
1101:     v      = aa;
1102:     vi     = aj;
1103:     for (i=1; i<n; i++) {
1104:       nz  = ai[i+1] - ai[i];
1105:       sum = b[r[i]];
1106:       PetscSparseDenseMinusDot(sum,tmp,v,vi,nz);
1107:       tmp[i] = sum;
1108:       v     += nz; vi += nz;
1109:     }
1110:     /* backward solve the upper triangular */
1111:     for (i=n-1; i>=0; i--) {
1112:       v   = aa + adiag[i+1]+1;
1113:       vi  = aj + adiag[i+1]+1;
1114:       nz  = adiag[i]-adiag[i+1]-1;
1115:       sum = tmp[i];
1116:       PetscSparseDenseMinusDot(sum,tmp,v,vi,nz);
1117:       x[c[i]] = tmp[i] = sum*v[nz]; /* v[nz] = aa[adiag[i]] */
1118:     }
1119:     b += ldb;
1120:     x += ldx;
1121:   }
1122:   ISRestoreIndices(isrow,&rout);
1123:   ISRestoreIndices(iscol,&cout);
1124:   MatDenseRestoreArrayRead(B,&b);
1125:   MatDenseRestoreArray(X,&x);
1126:   PetscLogFlops(B->cmap->n*(2.0*a->nz - n));
1127:   return(0);
1128: }

1130: PetscErrorCode MatSolve_SeqAIJ_InplaceWithPerm(Mat A,Vec bb,Vec xx)
1131: {
1132:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1133:   IS                iscol = a->col,isrow = a->row;
1134:   PetscErrorCode    ierr;
1135:   const PetscInt    *r,*c,*rout,*cout;
1136:   PetscInt          i, n = A->rmap->n,*vi,*ai = a->i,*aj = a->j;
1137:   PetscInt          nz,row;
1138:   PetscScalar       *x,*tmp,*tmps,sum;
1139:   const PetscScalar *b;
1140:   const MatScalar   *aa = a->a,*v;

1143:   if (!n) return(0);

1145:   VecGetArrayRead(bb,&b);
1146:   VecGetArrayWrite(xx,&x);
1147:   tmp  = a->solve_work;

1149:   ISGetIndices(isrow,&rout); r = rout;
1150:   ISGetIndices(iscol,&cout); c = cout + (n-1);

1152:   /* forward solve the lower triangular */
1153:   tmp[0] = b[*r++];
1154:   tmps   = tmp;
1155:   for (row=1; row<n; row++) {
1156:     i   = rout[row]; /* permuted row */
1157:     v   = aa + ai[i];
1158:     vi  = aj + ai[i];
1159:     nz  = a->diag[i] - ai[i];
1160:     sum = b[*r++];
1161:     PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
1162:     tmp[row] = sum;
1163:   }

1165:   /* backward solve the upper triangular */
1166:   for (row=n-1; row>=0; row--) {
1167:     i   = rout[row]; /* permuted row */
1168:     v   = aa + a->diag[i] + 1;
1169:     vi  = aj + a->diag[i] + 1;
1170:     nz  = ai[i+1] - a->diag[i] - 1;
1171:     sum = tmp[row];
1172:     PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
1173:     x[*c--] = tmp[row] = sum*aa[a->diag[i]];
1174:   }

1176:   ISRestoreIndices(isrow,&rout);
1177:   ISRestoreIndices(iscol,&cout);
1178:   VecRestoreArrayRead(bb,&b);
1179:   VecRestoreArrayWrite(xx,&x);
1180:   PetscLogFlops(2.0*a->nz - A->cmap->n);
1181:   return(0);
1182: }

1184: /* ----------------------------------------------------------- */
1185: #include <../src/mat/impls/aij/seq/ftn-kernels/fsolve.h>
1186: PetscErrorCode MatSolve_SeqAIJ_NaturalOrdering_inplace(Mat A,Vec bb,Vec xx)
1187: {
1188:   Mat_SeqAIJ        *a = (Mat_SeqAIJ*)A->data;
1189:   PetscErrorCode    ierr;
1190:   PetscInt          n   = A->rmap->n;
1191:   const PetscInt    *ai = a->i,*aj = a->j,*adiag = a->diag;
1192:   PetscScalar       *x;
1193:   const PetscScalar *b;
1194:   const MatScalar   *aa = a->a;
1195: #if !defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
1196:   PetscInt        adiag_i,i,nz,ai_i;
1197:   const PetscInt  *vi;
1198:   const MatScalar *v;
1199:   PetscScalar     sum;
1200: #endif

1203:   if (!n) return(0);

1205:   VecGetArrayRead(bb,&b);
1206:   VecGetArrayWrite(xx,&x);

1208: #if defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
1209:   fortransolveaij_(&n,x,ai,aj,adiag,aa,b);
1210: #else
1211:   /* forward solve the lower triangular */
1212:   x[0] = b[0];
1213:   for (i=1; i<n; i++) {
1214:     ai_i = ai[i];
1215:     v    = aa + ai_i;
1216:     vi   = aj + ai_i;
1217:     nz   = adiag[i] - ai_i;
1218:     sum  = b[i];
1219:     PetscSparseDenseMinusDot(sum,x,v,vi,nz);
1220:     x[i] = sum;
1221:   }

1223:   /* backward solve the upper triangular */
1224:   for (i=n-1; i>=0; i--) {
1225:     adiag_i = adiag[i];
1226:     v       = aa + adiag_i + 1;
1227:     vi      = aj + adiag_i + 1;
1228:     nz      = ai[i+1] - adiag_i - 1;
1229:     sum     = x[i];
1230:     PetscSparseDenseMinusDot(sum,x,v,vi,nz);
1231:     x[i] = sum*aa[adiag_i];
1232:   }
1233: #endif
1234:   PetscLogFlops(2.0*a->nz - A->cmap->n);
1235:   VecRestoreArrayRead(bb,&b);
1236:   VecRestoreArrayWrite(xx,&x);
1237:   return(0);
1238: }

1240: PetscErrorCode MatSolveAdd_SeqAIJ_inplace(Mat A,Vec bb,Vec yy,Vec xx)
1241: {
1242:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1243:   IS                iscol = a->col,isrow = a->row;
1244:   PetscErrorCode    ierr;
1245:   PetscInt          i, n = A->rmap->n,j;
1246:   PetscInt          nz;
1247:   const PetscInt    *rout,*cout,*r,*c,*vi,*ai = a->i,*aj = a->j;
1248:   PetscScalar       *x,*tmp,sum;
1249:   const PetscScalar *b;
1250:   const MatScalar   *aa = a->a,*v;

1253:   if (yy != xx) {VecCopy(yy,xx);}

1255:   VecGetArrayRead(bb,&b);
1256:   VecGetArray(xx,&x);
1257:   tmp  = a->solve_work;

1259:   ISGetIndices(isrow,&rout); r = rout;
1260:   ISGetIndices(iscol,&cout); c = cout + (n-1);

1262:   /* forward solve the lower triangular */
1263:   tmp[0] = b[*r++];
1264:   for (i=1; i<n; i++) {
1265:     v   = aa + ai[i];
1266:     vi  = aj + ai[i];
1267:     nz  = a->diag[i] - ai[i];
1268:     sum = b[*r++];
1269:     for (j=0; j<nz; j++) sum -= v[j]*tmp[vi[j]];
1270:     tmp[i] = sum;
1271:   }

1273:   /* backward solve the upper triangular */
1274:   for (i=n-1; i>=0; i--) {
1275:     v   = aa + a->diag[i] + 1;
1276:     vi  = aj + a->diag[i] + 1;
1277:     nz  = ai[i+1] - a->diag[i] - 1;
1278:     sum = tmp[i];
1279:     for (j=0; j<nz; j++) sum -= v[j]*tmp[vi[j]];
1280:     tmp[i]   = sum*aa[a->diag[i]];
1281:     x[*c--] += tmp[i];
1282:   }

1284:   ISRestoreIndices(isrow,&rout);
1285:   ISRestoreIndices(iscol,&cout);
1286:   VecRestoreArrayRead(bb,&b);
1287:   VecRestoreArray(xx,&x);
1288:   PetscLogFlops(2.0*a->nz);
1289:   return(0);
1290: }

1292: PetscErrorCode MatSolveAdd_SeqAIJ(Mat A,Vec bb,Vec yy,Vec xx)
1293: {
1294:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1295:   IS                iscol = a->col,isrow = a->row;
1296:   PetscErrorCode    ierr;
1297:   PetscInt          i, n = A->rmap->n,j;
1298:   PetscInt          nz;
1299:   const PetscInt    *rout,*cout,*r,*c,*vi,*ai = a->i,*aj = a->j,*adiag = a->diag;
1300:   PetscScalar       *x,*tmp,sum;
1301:   const PetscScalar *b;
1302:   const MatScalar   *aa = a->a,*v;

1305:   if (yy != xx) {VecCopy(yy,xx);}

1307:   VecGetArrayRead(bb,&b);
1308:   VecGetArray(xx,&x);
1309:   tmp  = a->solve_work;

1311:   ISGetIndices(isrow,&rout); r = rout;
1312:   ISGetIndices(iscol,&cout); c = cout;

1314:   /* forward solve the lower triangular */
1315:   tmp[0] = b[r[0]];
1316:   v      = aa;
1317:   vi     = aj;
1318:   for (i=1; i<n; i++) {
1319:     nz  = ai[i+1] - ai[i];
1320:     sum = b[r[i]];
1321:     for (j=0; j<nz; j++) sum -= v[j]*tmp[vi[j]];
1322:     tmp[i] = sum;
1323:     v     += nz;
1324:     vi    += nz;
1325:   }

1327:   /* backward solve the upper triangular */
1328:   v  = aa + adiag[n-1];
1329:   vi = aj + adiag[n-1];
1330:   for (i=n-1; i>=0; i--) {
1331:     nz  = adiag[i] - adiag[i+1] - 1;
1332:     sum = tmp[i];
1333:     for (j=0; j<nz; j++) sum -= v[j]*tmp[vi[j]];
1334:     tmp[i]   = sum*v[nz];
1335:     x[c[i]] += tmp[i];
1336:     v       += nz+1; vi += nz+1;
1337:   }

1339:   ISRestoreIndices(isrow,&rout);
1340:   ISRestoreIndices(iscol,&cout);
1341:   VecRestoreArrayRead(bb,&b);
1342:   VecRestoreArray(xx,&x);
1343:   PetscLogFlops(2.0*a->nz);
1344:   return(0);
1345: }

1347: PetscErrorCode MatSolveTranspose_SeqAIJ_inplace(Mat A,Vec bb,Vec xx)
1348: {
1349:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1350:   IS                iscol = a->col,isrow = a->row;
1351:   PetscErrorCode    ierr;
1352:   const PetscInt    *rout,*cout,*r,*c,*diag = a->diag,*ai = a->i,*aj = a->j,*vi;
1353:   PetscInt          i,n = A->rmap->n,j;
1354:   PetscInt          nz;
1355:   PetscScalar       *x,*tmp,s1;
1356:   const MatScalar   *aa = a->a,*v;
1357:   const PetscScalar *b;

1360:   VecGetArrayRead(bb,&b);
1361:   VecGetArrayWrite(xx,&x);
1362:   tmp  = a->solve_work;

1364:   ISGetIndices(isrow,&rout); r = rout;
1365:   ISGetIndices(iscol,&cout); c = cout;

1367:   /* copy the b into temp work space according to permutation */
1368:   for (i=0; i<n; i++) tmp[i] = b[c[i]];

1370:   /* forward solve the U^T */
1371:   for (i=0; i<n; i++) {
1372:     v   = aa + diag[i];
1373:     vi  = aj + diag[i] + 1;
1374:     nz  = ai[i+1] - diag[i] - 1;
1375:     s1  = tmp[i];
1376:     s1 *= (*v++);  /* multiply by inverse of diagonal entry */
1377:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1378:     tmp[i] = s1;
1379:   }

1381:   /* backward solve the L^T */
1382:   for (i=n-1; i>=0; i--) {
1383:     v  = aa + diag[i] - 1;
1384:     vi = aj + diag[i] - 1;
1385:     nz = diag[i] - ai[i];
1386:     s1 = tmp[i];
1387:     for (j=0; j>-nz; j--) tmp[vi[j]] -= s1*v[j];
1388:   }

1390:   /* copy tmp into x according to permutation */
1391:   for (i=0; i<n; i++) x[r[i]] = tmp[i];

1393:   ISRestoreIndices(isrow,&rout);
1394:   ISRestoreIndices(iscol,&cout);
1395:   VecRestoreArrayRead(bb,&b);
1396:   VecRestoreArrayWrite(xx,&x);

1398:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1399:   return(0);
1400: }

1402: PetscErrorCode MatSolveTranspose_SeqAIJ(Mat A,Vec bb,Vec xx)
1403: {
1404:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1405:   IS                iscol = a->col,isrow = a->row;
1406:   PetscErrorCode    ierr;
1407:   const PetscInt    *rout,*cout,*r,*c,*adiag = a->diag,*ai = a->i,*aj = a->j,*vi;
1408:   PetscInt          i,n = A->rmap->n,j;
1409:   PetscInt          nz;
1410:   PetscScalar       *x,*tmp,s1;
1411:   const MatScalar   *aa = a->a,*v;
1412:   const PetscScalar *b;

1415:   VecGetArrayRead(bb,&b);
1416:   VecGetArrayWrite(xx,&x);
1417:   tmp  = a->solve_work;

1419:   ISGetIndices(isrow,&rout); r = rout;
1420:   ISGetIndices(iscol,&cout); c = cout;

1422:   /* copy the b into temp work space according to permutation */
1423:   for (i=0; i<n; i++) tmp[i] = b[c[i]];

1425:   /* forward solve the U^T */
1426:   for (i=0; i<n; i++) {
1427:     v   = aa + adiag[i+1] + 1;
1428:     vi  = aj + adiag[i+1] + 1;
1429:     nz  = adiag[i] - adiag[i+1] - 1;
1430:     s1  = tmp[i];
1431:     s1 *= v[nz];  /* multiply by inverse of diagonal entry */
1432:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1433:     tmp[i] = s1;
1434:   }

1436:   /* backward solve the L^T */
1437:   for (i=n-1; i>=0; i--) {
1438:     v  = aa + ai[i];
1439:     vi = aj + ai[i];
1440:     nz = ai[i+1] - ai[i];
1441:     s1 = tmp[i];
1442:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1443:   }

1445:   /* copy tmp into x according to permutation */
1446:   for (i=0; i<n; i++) x[r[i]] = tmp[i];

1448:   ISRestoreIndices(isrow,&rout);
1449:   ISRestoreIndices(iscol,&cout);
1450:   VecRestoreArrayRead(bb,&b);
1451:   VecRestoreArrayWrite(xx,&x);

1453:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1454:   return(0);
1455: }

1457: PetscErrorCode MatSolveTransposeAdd_SeqAIJ_inplace(Mat A,Vec bb,Vec zz,Vec xx)
1458: {
1459:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1460:   IS                iscol = a->col,isrow = a->row;
1461:   PetscErrorCode    ierr;
1462:   const PetscInt    *rout,*cout,*r,*c,*diag = a->diag,*ai = a->i,*aj = a->j,*vi;
1463:   PetscInt          i,n = A->rmap->n,j;
1464:   PetscInt          nz;
1465:   PetscScalar       *x,*tmp,s1;
1466:   const MatScalar   *aa = a->a,*v;
1467:   const PetscScalar *b;

1470:   if (zz != xx) {VecCopy(zz,xx);}
1471:   VecGetArrayRead(bb,&b);
1472:   VecGetArray(xx,&x);
1473:   tmp  = a->solve_work;

1475:   ISGetIndices(isrow,&rout); r = rout;
1476:   ISGetIndices(iscol,&cout); c = cout;

1478:   /* copy the b into temp work space according to permutation */
1479:   for (i=0; i<n; i++) tmp[i] = b[c[i]];

1481:   /* forward solve the U^T */
1482:   for (i=0; i<n; i++) {
1483:     v   = aa + diag[i];
1484:     vi  = aj + diag[i] + 1;
1485:     nz  = ai[i+1] - diag[i] - 1;
1486:     s1  = tmp[i];
1487:     s1 *= (*v++);  /* multiply by inverse of diagonal entry */
1488:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1489:     tmp[i] = s1;
1490:   }

1492:   /* backward solve the L^T */
1493:   for (i=n-1; i>=0; i--) {
1494:     v  = aa + diag[i] - 1;
1495:     vi = aj + diag[i] - 1;
1496:     nz = diag[i] - ai[i];
1497:     s1 = tmp[i];
1498:     for (j=0; j>-nz; j--) tmp[vi[j]] -= s1*v[j];
1499:   }

1501:   /* copy tmp into x according to permutation */
1502:   for (i=0; i<n; i++) x[r[i]] += tmp[i];

1504:   ISRestoreIndices(isrow,&rout);
1505:   ISRestoreIndices(iscol,&cout);
1506:   VecRestoreArrayRead(bb,&b);
1507:   VecRestoreArray(xx,&x);

1509:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1510:   return(0);
1511: }

1513: PetscErrorCode MatSolveTransposeAdd_SeqAIJ(Mat A,Vec bb,Vec zz,Vec xx)
1514: {
1515:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1516:   IS                iscol = a->col,isrow = a->row;
1517:   PetscErrorCode    ierr;
1518:   const PetscInt    *rout,*cout,*r,*c,*adiag = a->diag,*ai = a->i,*aj = a->j,*vi;
1519:   PetscInt          i,n = A->rmap->n,j;
1520:   PetscInt          nz;
1521:   PetscScalar       *x,*tmp,s1;
1522:   const MatScalar   *aa = a->a,*v;
1523:   const PetscScalar *b;

1526:   if (zz != xx) {VecCopy(zz,xx);}
1527:   VecGetArrayRead(bb,&b);
1528:   VecGetArray(xx,&x);
1529:   tmp  = a->solve_work;

1531:   ISGetIndices(isrow,&rout); r = rout;
1532:   ISGetIndices(iscol,&cout); c = cout;

1534:   /* copy the b into temp work space according to permutation */
1535:   for (i=0; i<n; i++) tmp[i] = b[c[i]];

1537:   /* forward solve the U^T */
1538:   for (i=0; i<n; i++) {
1539:     v   = aa + adiag[i+1] + 1;
1540:     vi  = aj + adiag[i+1] + 1;
1541:     nz  = adiag[i] - adiag[i+1] - 1;
1542:     s1  = tmp[i];
1543:     s1 *= v[nz];  /* multiply by inverse of diagonal entry */
1544:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1545:     tmp[i] = s1;
1546:   }


1549:   /* backward solve the L^T */
1550:   for (i=n-1; i>=0; i--) {
1551:     v  = aa + ai[i];
1552:     vi = aj + ai[i];
1553:     nz = ai[i+1] - ai[i];
1554:     s1 = tmp[i];
1555:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1556:   }

1558:   /* copy tmp into x according to permutation */
1559:   for (i=0; i<n; i++) x[r[i]] += tmp[i];

1561:   ISRestoreIndices(isrow,&rout);
1562:   ISRestoreIndices(iscol,&cout);
1563:   VecRestoreArrayRead(bb,&b);
1564:   VecRestoreArray(xx,&x);

1566:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1567:   return(0);
1568: }

1570: /* ----------------------------------------------------------------*/

1572: /*
1573:    ilu() under revised new data structure.
1574:    Factored arrays bj and ba are stored as
1575:      L(0,:), L(1,:), ...,L(n-1,:),  U(n-1,:),...,U(i,:),U(i-1,:),...,U(0,:)

1577:    bi=fact->i is an array of size n+1, in which
1578:    bi+
1579:      bi[i]:  points to 1st entry of L(i,:),i=0,...,n-1
1580:      bi[n]:  points to L(n-1,n-1)+1

1582:   bdiag=fact->diag is an array of size n+1,in which
1583:      bdiag[i]: points to diagonal of U(i,:), i=0,...,n-1
1584:      bdiag[n]: points to entry of U(n-1,0)-1

1586:    U(i,:) contains bdiag[i] as its last entry, i.e.,
1587:     U(i,:) = (u[i,i+1],...,u[i,n-1],diag[i])
1588: */
1589: PetscErrorCode MatILUFactorSymbolic_SeqAIJ_ilu0(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
1590: {
1591:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data,*b;
1593:   const PetscInt n=A->rmap->n,*ai=a->i,*aj,*adiag=a->diag;
1594:   PetscInt       i,j,k=0,nz,*bi,*bj,*bdiag;
1595:   IS             isicol;

1598:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
1599:   MatDuplicateNoCreate_SeqAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_FALSE);
1600:   b    = (Mat_SeqAIJ*)(fact)->data;

1602:   /* allocate matrix arrays for new data structure */
1603:   PetscMalloc3(ai[n]+1,&b->a,ai[n]+1,&b->j,n+1,&b->i);
1604:   PetscLogObjectMemory((PetscObject)fact,ai[n]*(sizeof(PetscScalar)+sizeof(PetscInt))+(n+1)*sizeof(PetscInt));

1606:   b->singlemalloc = PETSC_TRUE;
1607:   if (!b->diag) {
1608:     PetscMalloc1(n+1,&b->diag);
1609:     PetscLogObjectMemory((PetscObject)fact,(n+1)*sizeof(PetscInt));
1610:   }
1611:   bdiag = b->diag;

1613:   if (n > 0) {
1614:     PetscArrayzero(b->a,ai[n]);
1615:   }

1617:   /* set bi and bj with new data structure */
1618:   bi = b->i;
1619:   bj = b->j;

1621:   /* L part */
1622:   bi[0] = 0;
1623:   for (i=0; i<n; i++) {
1624:     nz      = adiag[i] - ai[i];
1625:     bi[i+1] = bi[i] + nz;
1626:     aj      = a->j + ai[i];
1627:     for (j=0; j<nz; j++) {
1628:       /*   *bj = aj[j]; bj++; */
1629:       bj[k++] = aj[j];
1630:     }
1631:   }

1633:   /* U part */
1634:   bdiag[n] = bi[n]-1;
1635:   for (i=n-1; i>=0; i--) {
1636:     nz = ai[i+1] - adiag[i] - 1;
1637:     aj = a->j + adiag[i] + 1;
1638:     for (j=0; j<nz; j++) {
1639:       /*      *bj = aj[j]; bj++; */
1640:       bj[k++] = aj[j];
1641:     }
1642:     /* diag[i] */
1643:     /*    *bj = i; bj++; */
1644:     bj[k++]  = i;
1645:     bdiag[i] = bdiag[i+1] + nz + 1;
1646:   }

1648:   fact->factortype             = MAT_FACTOR_ILU;
1649:   fact->info.factor_mallocs    = 0;
1650:   fact->info.fill_ratio_given  = info->fill;
1651:   fact->info.fill_ratio_needed = 1.0;
1652:   fact->ops->lufactornumeric   = MatLUFactorNumeric_SeqAIJ;
1653:   MatSeqAIJCheckInode_FactorLU(fact);

1655:   b       = (Mat_SeqAIJ*)(fact)->data;
1656:   b->row  = isrow;
1657:   b->col  = iscol;
1658:   b->icol = isicol;
1659:   PetscMalloc1(fact->rmap->n+1,&b->solve_work);
1660:   PetscObjectReference((PetscObject)isrow);
1661:   PetscObjectReference((PetscObject)iscol);
1662:   return(0);
1663: }

1665: PetscErrorCode MatILUFactorSymbolic_SeqAIJ(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
1666: {
1667:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data,*b;
1668:   IS                 isicol;
1669:   PetscErrorCode     ierr;
1670:   const PetscInt     *r,*ic;
1671:   PetscInt           n=A->rmap->n,*ai=a->i,*aj=a->j;
1672:   PetscInt           *bi,*cols,nnz,*cols_lvl;
1673:   PetscInt           *bdiag,prow,fm,nzbd,reallocs=0,dcount=0;
1674:   PetscInt           i,levels,diagonal_fill;
1675:   PetscBool          col_identity,row_identity,missing;
1676:   PetscReal          f;
1677:   PetscInt           nlnk,*lnk,*lnk_lvl=NULL;
1678:   PetscBT            lnkbt;
1679:   PetscInt           nzi,*bj,**bj_ptr,**bjlvl_ptr;
1680:   PetscFreeSpaceList free_space    =NULL,current_space=NULL;
1681:   PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL;

1684:   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);
1685:   MatMissingDiagonal(A,&missing,&i);
1686:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);

1688:   levels = (PetscInt)info->levels;
1689:   ISIdentity(isrow,&row_identity);
1690:   ISIdentity(iscol,&col_identity);
1691:   if (!levels && row_identity && col_identity) {
1692:     /* special case: ilu(0) with natural ordering */
1693:     MatILUFactorSymbolic_SeqAIJ_ilu0(fact,A,isrow,iscol,info);
1694:     if (a->inode.size) {
1695:       fact->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
1696:     }
1697:     return(0);
1698:   }

1700:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
1701:   ISGetIndices(isrow,&r);
1702:   ISGetIndices(isicol,&ic);

1704:   /* get new row and diagonal pointers, must be allocated separately because they will be given to the Mat_SeqAIJ and freed separately */
1705:   PetscMalloc1(n+1,&bi);
1706:   PetscMalloc1(n+1,&bdiag);
1707:   bi[0] = bdiag[0] = 0;
1708:   PetscMalloc2(n,&bj_ptr,n,&bjlvl_ptr);

1710:   /* create a linked list for storing column indices of the active row */
1711:   nlnk = n + 1;
1712:   PetscIncompleteLLCreate(n,n,nlnk,lnk,lnk_lvl,lnkbt);

1714:   /* initial FreeSpace size is f*(ai[n]+1) */
1715:   f                 = info->fill;
1716:   diagonal_fill     = (PetscInt)info->diagonal_fill;
1717:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
1718:   current_space     = free_space;
1719:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space_lvl);
1720:   current_space_lvl = free_space_lvl;
1721:   for (i=0; i<n; i++) {
1722:     nzi = 0;
1723:     /* copy current row into linked list */
1724:     nnz = ai[r[i]+1] - ai[r[i]];
1725:     if (!nnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",r[i],i);
1726:     cols   = aj + ai[r[i]];
1727:     lnk[i] = -1; /* marker to indicate if diagonal exists */
1728:     PetscIncompleteLLInit(nnz,cols,n,ic,nlnk,lnk,lnk_lvl,lnkbt);
1729:     nzi   += nlnk;

1731:     /* make sure diagonal entry is included */
1732:     if (diagonal_fill && lnk[i] == -1) {
1733:       fm = n;
1734:       while (lnk[fm] < i) fm = lnk[fm];
1735:       lnk[i]     = lnk[fm]; /* insert diagonal into linked list */
1736:       lnk[fm]    = i;
1737:       lnk_lvl[i] = 0;
1738:       nzi++; dcount++;
1739:     }

1741:     /* add pivot rows into the active row */
1742:     nzbd = 0;
1743:     prow = lnk[n];
1744:     while (prow < i) {
1745:       nnz      = bdiag[prow];
1746:       cols     = bj_ptr[prow] + nnz + 1;
1747:       cols_lvl = bjlvl_ptr[prow] + nnz + 1;
1748:       nnz      = bi[prow+1] - bi[prow] - nnz - 1;
1749:       PetscILULLAddSorted(nnz,cols,levels,cols_lvl,prow,nlnk,lnk,lnk_lvl,lnkbt,prow);
1750:       nzi     += nlnk;
1751:       prow     = lnk[prow];
1752:       nzbd++;
1753:     }
1754:     bdiag[i] = nzbd;
1755:     bi[i+1]  = bi[i] + nzi;
1756:     /* if free space is not available, make more free space */
1757:     if (current_space->local_remaining<nzi) {
1758:       nnz  = PetscIntMultTruncate(2,PetscIntMultTruncate(nzi,n - i)); /* estimated and max additional space needed */
1759:       PetscFreeSpaceGet(nnz,&current_space);
1760:       PetscFreeSpaceGet(nnz,&current_space_lvl);
1761:       reallocs++;
1762:     }

1764:     /* copy data into free_space and free_space_lvl, then initialize lnk */
1765:     PetscIncompleteLLClean(n,n,nzi,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);
1766:     bj_ptr[i]    = current_space->array;
1767:     bjlvl_ptr[i] = current_space_lvl->array;

1769:     /* make sure the active row i has diagonal entry */
1770:     if (*(bj_ptr[i]+bdiag[i]) != i) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Row %D has missing diagonal in factored matrix\ntry running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill",i);

1772:     current_space->array               += nzi;
1773:     current_space->local_used          += nzi;
1774:     current_space->local_remaining     -= nzi;
1775:     current_space_lvl->array           += nzi;
1776:     current_space_lvl->local_used      += nzi;
1777:     current_space_lvl->local_remaining -= nzi;
1778:   }

1780:   ISRestoreIndices(isrow,&r);
1781:   ISRestoreIndices(isicol,&ic);
1782:   /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
1783:   PetscMalloc1(bi[n]+1,&bj);
1784:   PetscFreeSpaceContiguous_LU(&free_space,bj,n,bi,bdiag);

1786:   PetscIncompleteLLDestroy(lnk,lnkbt);
1787:   PetscFreeSpaceDestroy(free_space_lvl);
1788:   PetscFree2(bj_ptr,bjlvl_ptr);

1790: #if defined(PETSC_USE_INFO)
1791:   {
1792:     PetscReal af = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]);
1793:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);
1794:     PetscInfo1(A,"Run with -[sub_]pc_factor_fill %g or use \n",(double)af);
1795:     PetscInfo1(A,"PCFactorSetFill([sub]pc,%g);\n",(double)af);
1796:     PetscInfo(A,"for best performance.\n");
1797:     if (diagonal_fill) {
1798:       PetscInfo1(A,"Detected and replaced %D missing diagonals\n",dcount);
1799:     }
1800:   }
1801: #endif
1802:   /* put together the new matrix */
1803:   MatSeqAIJSetPreallocation_SeqAIJ(fact,MAT_SKIP_ALLOCATION,NULL);
1804:   PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol);
1805:   b    = (Mat_SeqAIJ*)(fact)->data;

1807:   b->free_a       = PETSC_TRUE;
1808:   b->free_ij      = PETSC_TRUE;
1809:   b->singlemalloc = PETSC_FALSE;

1811:   PetscMalloc1(bdiag[0]+1,&b->a);

1813:   b->j    = bj;
1814:   b->i    = bi;
1815:   b->diag = bdiag;
1816:   b->ilen = NULL;
1817:   b->imax = NULL;
1818:   b->row  = isrow;
1819:   b->col  = iscol;
1820:   PetscObjectReference((PetscObject)isrow);
1821:   PetscObjectReference((PetscObject)iscol);
1822:   b->icol = isicol;

1824:   PetscMalloc1(n+1,&b->solve_work);
1825:   /* In b structure:  Free imax, ilen, old a, old j.
1826:      Allocate bdiag, solve_work, new a, new j */
1827:   PetscLogObjectMemory((PetscObject)fact,(bdiag[0]+1)*(sizeof(PetscInt)+sizeof(PetscScalar)));
1828:   b->maxnz = b->nz = bdiag[0]+1;

1830:   (fact)->info.factor_mallocs    = reallocs;
1831:   (fact)->info.fill_ratio_given  = f;
1832:   (fact)->info.fill_ratio_needed = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]);
1833:   (fact)->ops->lufactornumeric   = MatLUFactorNumeric_SeqAIJ;
1834:   if (a->inode.size) {
1835:     (fact)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
1836:   }
1837:   MatSeqAIJCheckInode_FactorLU(fact);
1838:   return(0);
1839: }

1841: PetscErrorCode MatILUFactorSymbolic_SeqAIJ_inplace(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
1842: {
1843:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data,*b;
1844:   IS                 isicol;
1845:   PetscErrorCode     ierr;
1846:   const PetscInt     *r,*ic;
1847:   PetscInt           n=A->rmap->n,*ai=a->i,*aj=a->j;
1848:   PetscInt           *bi,*cols,nnz,*cols_lvl;
1849:   PetscInt           *bdiag,prow,fm,nzbd,reallocs=0,dcount=0;
1850:   PetscInt           i,levels,diagonal_fill;
1851:   PetscBool          col_identity,row_identity;
1852:   PetscReal          f;
1853:   PetscInt           nlnk,*lnk,*lnk_lvl=NULL;
1854:   PetscBT            lnkbt;
1855:   PetscInt           nzi,*bj,**bj_ptr,**bjlvl_ptr;
1856:   PetscFreeSpaceList free_space    =NULL,current_space=NULL;
1857:   PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL;
1858:   PetscBool          missing;

1861:   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);
1862:   MatMissingDiagonal(A,&missing,&i);
1863:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);

1865:   f             = info->fill;
1866:   levels        = (PetscInt)info->levels;
1867:   diagonal_fill = (PetscInt)info->diagonal_fill;

1869:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);

1871:   ISIdentity(isrow,&row_identity);
1872:   ISIdentity(iscol,&col_identity);
1873:   if (!levels && row_identity && col_identity) { /* special case: ilu(0) with natural ordering */
1874:     MatDuplicateNoCreate_SeqAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_TRUE);

1876:     (fact)->ops->lufactornumeric =  MatLUFactorNumeric_SeqAIJ_inplace;
1877:     if (a->inode.size) {
1878:       (fact)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode_inplace;
1879:     }
1880:     fact->factortype               = MAT_FACTOR_ILU;
1881:     (fact)->info.factor_mallocs    = 0;
1882:     (fact)->info.fill_ratio_given  = info->fill;
1883:     (fact)->info.fill_ratio_needed = 1.0;

1885:     b    = (Mat_SeqAIJ*)(fact)->data;
1886:     b->row  = isrow;
1887:     b->col  = iscol;
1888:     b->icol = isicol;
1889:     PetscMalloc1((fact)->rmap->n+1,&b->solve_work);
1890:     PetscObjectReference((PetscObject)isrow);
1891:     PetscObjectReference((PetscObject)iscol);
1892:     return(0);
1893:   }

1895:   ISGetIndices(isrow,&r);
1896:   ISGetIndices(isicol,&ic);

1898:   /* get new row and diagonal pointers, must be allocated separately because they will be given to the Mat_SeqAIJ and freed separately */
1899:   PetscMalloc1(n+1,&bi);
1900:   PetscMalloc1(n+1,&bdiag);
1901:   bi[0] = bdiag[0] = 0;

1903:   PetscMalloc2(n,&bj_ptr,n,&bjlvl_ptr);

1905:   /* create a linked list for storing column indices of the active row */
1906:   nlnk = n + 1;
1907:   PetscIncompleteLLCreate(n,n,nlnk,lnk,lnk_lvl,lnkbt);

1909:   /* initial FreeSpace size is f*(ai[n]+1) */
1910:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
1911:   current_space     = free_space;
1912:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space_lvl);
1913:   current_space_lvl = free_space_lvl;

1915:   for (i=0; i<n; i++) {
1916:     nzi = 0;
1917:     /* copy current row into linked list */
1918:     nnz = ai[r[i]+1] - ai[r[i]];
1919:     if (!nnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",r[i],i);
1920:     cols   = aj + ai[r[i]];
1921:     lnk[i] = -1; /* marker to indicate if diagonal exists */
1922:     PetscIncompleteLLInit(nnz,cols,n,ic,nlnk,lnk,lnk_lvl,lnkbt);
1923:     nzi   += nlnk;

1925:     /* make sure diagonal entry is included */
1926:     if (diagonal_fill && lnk[i] == -1) {
1927:       fm = n;
1928:       while (lnk[fm] < i) fm = lnk[fm];
1929:       lnk[i]     = lnk[fm]; /* insert diagonal into linked list */
1930:       lnk[fm]    = i;
1931:       lnk_lvl[i] = 0;
1932:       nzi++; dcount++;
1933:     }

1935:     /* add pivot rows into the active row */
1936:     nzbd = 0;
1937:     prow = lnk[n];
1938:     while (prow < i) {
1939:       nnz      = bdiag[prow];
1940:       cols     = bj_ptr[prow] + nnz + 1;
1941:       cols_lvl = bjlvl_ptr[prow] + nnz + 1;
1942:       nnz      = bi[prow+1] - bi[prow] - nnz - 1;
1943:       PetscILULLAddSorted(nnz,cols,levels,cols_lvl,prow,nlnk,lnk,lnk_lvl,lnkbt,prow);
1944:       nzi     += nlnk;
1945:       prow     = lnk[prow];
1946:       nzbd++;
1947:     }
1948:     bdiag[i] = nzbd;
1949:     bi[i+1]  = bi[i] + nzi;

1951:     /* if free space is not available, make more free space */
1952:     if (current_space->local_remaining<nzi) {
1953:       nnz  = PetscIntMultTruncate(nzi,n - i); /* estimated and max additional space needed */
1954:       PetscFreeSpaceGet(nnz,&current_space);
1955:       PetscFreeSpaceGet(nnz,&current_space_lvl);
1956:       reallocs++;
1957:     }

1959:     /* copy data into free_space and free_space_lvl, then initialize lnk */
1960:     PetscIncompleteLLClean(n,n,nzi,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);
1961:     bj_ptr[i]    = current_space->array;
1962:     bjlvl_ptr[i] = current_space_lvl->array;

1964:     /* make sure the active row i has diagonal entry */
1965:     if (*(bj_ptr[i]+bdiag[i]) != i) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Row %D has missing diagonal in factored matrix\ntry running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill",i);

1967:     current_space->array               += nzi;
1968:     current_space->local_used          += nzi;
1969:     current_space->local_remaining     -= nzi;
1970:     current_space_lvl->array           += nzi;
1971:     current_space_lvl->local_used      += nzi;
1972:     current_space_lvl->local_remaining -= nzi;
1973:   }

1975:   ISRestoreIndices(isrow,&r);
1976:   ISRestoreIndices(isicol,&ic);

1978:   /* destroy list of free space and other temporary arrays */
1979:   PetscMalloc1(bi[n]+1,&bj);
1980:   PetscFreeSpaceContiguous(&free_space,bj); /* copy free_space -> bj */
1981:   PetscIncompleteLLDestroy(lnk,lnkbt);
1982:   PetscFreeSpaceDestroy(free_space_lvl);
1983:   PetscFree2(bj_ptr,bjlvl_ptr);

1985: #if defined(PETSC_USE_INFO)
1986:   {
1987:     PetscReal af = ((PetscReal)bi[n])/((PetscReal)ai[n]);
1988:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);
1989:     PetscInfo1(A,"Run with -[sub_]pc_factor_fill %g or use \n",(double)af);
1990:     PetscInfo1(A,"PCFactorSetFill([sub]pc,%g);\n",(double)af);
1991:     PetscInfo(A,"for best performance.\n");
1992:     if (diagonal_fill) {
1993:       PetscInfo1(A,"Detected and replaced %D missing diagonals\n",dcount);
1994:     }
1995:   }
1996: #endif

1998:   /* put together the new matrix */
1999:   MatSeqAIJSetPreallocation_SeqAIJ(fact,MAT_SKIP_ALLOCATION,NULL);
2000:   PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol);
2001:   b    = (Mat_SeqAIJ*)(fact)->data;

2003:   b->free_a       = PETSC_TRUE;
2004:   b->free_ij      = PETSC_TRUE;
2005:   b->singlemalloc = PETSC_FALSE;

2007:   PetscMalloc1(bi[n],&b->a);
2008:   b->j = bj;
2009:   b->i = bi;
2010:   for (i=0; i<n; i++) bdiag[i] += bi[i];
2011:   b->diag = bdiag;
2012:   b->ilen = NULL;
2013:   b->imax = NULL;
2014:   b->row  = isrow;
2015:   b->col  = iscol;
2016:   PetscObjectReference((PetscObject)isrow);
2017:   PetscObjectReference((PetscObject)iscol);
2018:   b->icol = isicol;
2019:   PetscMalloc1(n+1,&b->solve_work);
2020:   /* In b structure:  Free imax, ilen, old a, old j.
2021:      Allocate bdiag, solve_work, new a, new j */
2022:   PetscLogObjectMemory((PetscObject)fact,(bi[n]-n) * (sizeof(PetscInt)+sizeof(PetscScalar)));
2023:   b->maxnz = b->nz = bi[n];

2025:   (fact)->info.factor_mallocs    = reallocs;
2026:   (fact)->info.fill_ratio_given  = f;
2027:   (fact)->info.fill_ratio_needed = ((PetscReal)bi[n])/((PetscReal)ai[n]);
2028:   (fact)->ops->lufactornumeric   =  MatLUFactorNumeric_SeqAIJ_inplace;
2029:   if (a->inode.size) {
2030:     (fact)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode_inplace;
2031:   }
2032:   return(0);
2033: }

2035: PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ(Mat B,Mat A,const MatFactorInfo *info)
2036: {
2037:   Mat            C = B;
2038:   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data;
2039:   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
2040:   IS             ip=b->row,iip = b->icol;
2042:   const PetscInt *rip,*riip;
2043:   PetscInt       i,j,mbs=A->rmap->n,*bi=b->i,*bj=b->j,*bdiag=b->diag,*bjtmp;
2044:   PetscInt       *ai=a->i,*aj=a->j;
2045:   PetscInt       k,jmin,jmax,*c2r,*il,col,nexti,ili,nz;
2046:   MatScalar      *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
2047:   PetscBool      perm_identity;
2048:   FactorShiftCtx sctx;
2049:   PetscReal      rs;
2050:   MatScalar      d,*v;

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

2056:   if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
2057:     sctx.shift_top = info->zeropivot;
2058:     for (i=0; i<mbs; i++) {
2059:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
2060:       d  = (aa)[a->diag[i]];
2061:       rs = -PetscAbsScalar(d) - PetscRealPart(d);
2062:       v  = aa+ai[i];
2063:       nz = ai[i+1] - ai[i];
2064:       for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
2065:       if (rs>sctx.shift_top) sctx.shift_top = rs;
2066:     }
2067:     sctx.shift_top *= 1.1;
2068:     sctx.nshift_max = 5;
2069:     sctx.shift_lo   = 0.;
2070:     sctx.shift_hi   = 1.;
2071:   }

2073:   ISGetIndices(ip,&rip);
2074:   ISGetIndices(iip,&riip);

2076:   /* allocate working arrays
2077:      c2r: linked list, keep track of pivot rows for a given column. c2r[col]: head of the list for a given col
2078:      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
2079:   */
2080:   PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&c2r);

2082:   do {
2083:     sctx.newshift = PETSC_FALSE;

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

2088:     for (k = 0; k<mbs; k++) {
2089:       /* zero rtmp */
2090:       nz    = bi[k+1] - bi[k];
2091:       bjtmp = bj + bi[k];
2092:       for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

2094:       /* load in initial unfactored row */
2095:       bval = ba + bi[k];
2096:       jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
2097:       for (j = jmin; j < jmax; j++) {
2098:         col = riip[aj[j]];
2099:         if (col >= k) { /* only take upper triangular entry */
2100:           rtmp[col] = aa[j];
2101:           *bval++   = 0.0; /* for in-place factorization */
2102:         }
2103:       }
2104:       /* shift the diagonal of the matrix: ZeropivotApply() */
2105:       rtmp[k] += sctx.shift_amount;  /* shift the diagonal of the matrix */

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

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

2114:         /* compute multiplier, update diag(k) and U(i,k) */
2115:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
2116:         uikdi   = -ba[ili]*ba[bdiag[i]]; /* diagonal(k) */
2117:         dk     += uikdi*ba[ili]; /* update diag[k] */
2118:         ba[ili] = uikdi; /* -U(i,k) */

2120:         /* add multiple of row i to k-th row */
2121:         jmin = ili + 1; jmax = bi[i+1];
2122:         if (jmin < jmax) {
2123:           for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
2124:           /* update il and c2r for row i */
2125:           il[i] = jmin;
2126:           j     = bj[jmin]; c2r[i] = c2r[j]; c2r[j] = i;
2127:         }
2128:         i = nexti;
2129:       }

2131:       /* copy data into U(k,:) */
2132:       rs   = 0.0;
2133:       jmin = bi[k]; jmax = bi[k+1]-1;
2134:       if (jmin < jmax) {
2135:         for (j=jmin; j<jmax; j++) {
2136:           col = bj[j]; ba[j] = rtmp[col]; rs += PetscAbsScalar(ba[j]);
2137:         }
2138:         /* add the k-th row into il and c2r */
2139:         il[k] = jmin;
2140:         i     = bj[jmin]; c2r[k] = c2r[i]; c2r[i] = k;
2141:       }

2143:       /* MatPivotCheck() */
2144:       sctx.rs = rs;
2145:       sctx.pv = dk;
2146:       MatPivotCheck(B,A,info,&sctx,i);
2147:       if (sctx.newshift) break;
2148:       dk = sctx.pv;

2150:       ba[bdiag[k]] = 1.0/dk; /* U(k,k) */
2151:     }
2152:   } while (sctx.newshift);

2154:   PetscFree3(rtmp,il,c2r);
2155:   ISRestoreIndices(ip,&rip);
2156:   ISRestoreIndices(iip,&riip);

2158:   ISIdentity(ip,&perm_identity);
2159:   if (perm_identity) {
2160:     B->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering;
2161:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
2162:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
2163:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;
2164:   } else {
2165:     B->ops->solve          = MatSolve_SeqSBAIJ_1;
2166:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1;
2167:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1;
2168:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1;
2169:   }

2171:   C->assembled    = PETSC_TRUE;
2172:   C->preallocated = PETSC_TRUE;

2174:   PetscLogFlops(C->rmap->n);

2176:   /* MatPivotView() */
2177:   if (sctx.nshift) {
2178:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
2179:       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);
2180:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
2181:       PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
2182:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
2183:       PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %g\n",sctx.nshift,(double)info->shiftamount);
2184:     }
2185:   }
2186:   return(0);
2187: }

2189: PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ_inplace(Mat B,Mat A,const MatFactorInfo *info)
2190: {
2191:   Mat            C = B;
2192:   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data;
2193:   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
2194:   IS             ip=b->row,iip = b->icol;
2196:   const PetscInt *rip,*riip;
2197:   PetscInt       i,j,mbs=A->rmap->n,*bi=b->i,*bj=b->j,*bcol,*bjtmp;
2198:   PetscInt       *ai=a->i,*aj=a->j;
2199:   PetscInt       k,jmin,jmax,*jl,*il,col,nexti,ili,nz;
2200:   MatScalar      *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
2201:   PetscBool      perm_identity;
2202:   FactorShiftCtx sctx;
2203:   PetscReal      rs;
2204:   MatScalar      d,*v;

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

2210:   if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
2211:     sctx.shift_top = info->zeropivot;
2212:     for (i=0; i<mbs; i++) {
2213:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
2214:       d  = (aa)[a->diag[i]];
2215:       rs = -PetscAbsScalar(d) - PetscRealPart(d);
2216:       v  = aa+ai[i];
2217:       nz = ai[i+1] - ai[i];
2218:       for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
2219:       if (rs>sctx.shift_top) sctx.shift_top = rs;
2220:     }
2221:     sctx.shift_top *= 1.1;
2222:     sctx.nshift_max = 5;
2223:     sctx.shift_lo   = 0.;
2224:     sctx.shift_hi   = 1.;
2225:   }

2227:   ISGetIndices(ip,&rip);
2228:   ISGetIndices(iip,&riip);

2230:   /* initialization */
2231:   PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&jl);

2233:   do {
2234:     sctx.newshift = PETSC_FALSE;

2236:     for (i=0; i<mbs; i++) jl[i] = mbs;
2237:     il[0] = 0;

2239:     for (k = 0; k<mbs; k++) {
2240:       /* zero rtmp */
2241:       nz    = bi[k+1] - bi[k];
2242:       bjtmp = bj + bi[k];
2243:       for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

2245:       bval = ba + bi[k];
2246:       /* initialize k-th row by the perm[k]-th row of A */
2247:       jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
2248:       for (j = jmin; j < jmax; j++) {
2249:         col = riip[aj[j]];
2250:         if (col >= k) { /* only take upper triangular entry */
2251:           rtmp[col] = aa[j];
2252:           *bval++   = 0.0; /* for in-place factorization */
2253:         }
2254:       }
2255:       /* shift the diagonal of the matrix */
2256:       if (sctx.nshift) rtmp[k] += sctx.shift_amount;

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

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

2265:         /* compute multiplier, update diag(k) and U(i,k) */
2266:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
2267:         uikdi   = -ba[ili]*ba[bi[i]]; /* diagonal(k) */
2268:         dk     += uikdi*ba[ili];
2269:         ba[ili] = uikdi; /* -U(i,k) */

2271:         /* add multiple of row i to k-th row */
2272:         jmin = ili + 1; jmax = bi[i+1];
2273:         if (jmin < jmax) {
2274:           for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
2275:           /* update il and jl for row i */
2276:           il[i] = jmin;
2277:           j     = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
2278:         }
2279:         i = nexti;
2280:       }

2282:       /* shift the diagonals when zero pivot is detected */
2283:       /* compute rs=sum of abs(off-diagonal) */
2284:       rs   = 0.0;
2285:       jmin = bi[k]+1;
2286:       nz   = bi[k+1] - jmin;
2287:       bcol = bj + jmin;
2288:       for (j=0; j<nz; j++) {
2289:         rs += PetscAbsScalar(rtmp[bcol[j]]);
2290:       }

2292:       sctx.rs = rs;
2293:       sctx.pv = dk;
2294:       MatPivotCheck(B,A,info,&sctx,k);
2295:       if (sctx.newshift) break;
2296:       dk = sctx.pv;

2298:       /* copy data into U(k,:) */
2299:       ba[bi[k]] = 1.0/dk; /* U(k,k) */
2300:       jmin      = bi[k]+1; jmax = bi[k+1];
2301:       if (jmin < jmax) {
2302:         for (j=jmin; j<jmax; j++) {
2303:           col = bj[j]; ba[j] = rtmp[col];
2304:         }
2305:         /* add the k-th row into il and jl */
2306:         il[k] = jmin;
2307:         i     = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
2308:       }
2309:     }
2310:   } while (sctx.newshift);

2312:   PetscFree3(rtmp,il,jl);
2313:   ISRestoreIndices(ip,&rip);
2314:   ISRestoreIndices(iip,&riip);

2316:   ISIdentity(ip,&perm_identity);
2317:   if (perm_identity) {
2318:     B->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2319:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2320:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2321:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2322:   } else {
2323:     B->ops->solve          = MatSolve_SeqSBAIJ_1_inplace;
2324:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
2325:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_inplace;
2326:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_inplace;
2327:   }

2329:   C->assembled    = PETSC_TRUE;
2330:   C->preallocated = PETSC_TRUE;

2332:   PetscLogFlops(C->rmap->n);
2333:   if (sctx.nshift) {
2334:     if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
2335:       PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
2336:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
2337:       PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
2338:     }
2339:   }
2340:   return(0);
2341: }

2343: /*
2344:    icc() under revised new data structure.
2345:    Factored arrays bj and ba are stored as
2346:      U(0,:),...,U(i,:),U(n-1,:)

2348:    ui=fact->i is an array of size n+1, in which
2349:    ui+
2350:      ui[i]:  points to 1st entry of U(i,:),i=0,...,n-1
2351:      ui[n]:  points to U(n-1,n-1)+1

2353:   udiag=fact->diag is an array of size n,in which
2354:      udiag[i]: points to diagonal of U(i,:), i=0,...,n-1

2356:    U(i,:) contains udiag[i] as its last entry, i.e.,
2357:     U(i,:) = (u[i,i+1],...,u[i,n-1],diag[i])
2358: */

2360: PetscErrorCode MatICCFactorSymbolic_SeqAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
2361: {
2362:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data;
2363:   Mat_SeqSBAIJ       *b;
2364:   PetscErrorCode     ierr;
2365:   PetscBool          perm_identity,missing;
2366:   PetscInt           reallocs=0,i,*ai=a->i,*aj=a->j,am=A->rmap->n,*ui,*udiag;
2367:   const PetscInt     *rip,*riip;
2368:   PetscInt           jmin,jmax,nzk,k,j,*jl,prow,*il,nextprow;
2369:   PetscInt           nlnk,*lnk,*lnk_lvl=NULL,d;
2370:   PetscInt           ncols,ncols_upper,*cols,*ajtmp,*uj,**uj_ptr,**uj_lvl_ptr;
2371:   PetscReal          fill          =info->fill,levels=info->levels;
2372:   PetscFreeSpaceList free_space    =NULL,current_space=NULL;
2373:   PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL;
2374:   PetscBT            lnkbt;
2375:   IS                 iperm;

2378:   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);
2379:   MatMissingDiagonal(A,&missing,&d);
2380:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
2381:   ISIdentity(perm,&perm_identity);
2382:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);

2384:   PetscMalloc1(am+1,&ui);
2385:   PetscMalloc1(am+1,&udiag);
2386:   ui[0] = 0;

2388:   /* ICC(0) without matrix ordering: simply rearrange column indices */
2389:   if (!levels && perm_identity) {
2390:     for (i=0; i<am; i++) {
2391:       ncols    = ai[i+1] - a->diag[i];
2392:       ui[i+1]  = ui[i] + ncols;
2393:       udiag[i] = ui[i+1] - 1; /* points to the last entry of U(i,:) */
2394:     }
2395:     PetscMalloc1(ui[am]+1,&uj);
2396:     cols = uj;
2397:     for (i=0; i<am; i++) {
2398:       aj    = a->j + a->diag[i] + 1; /* 1st entry of U(i,:) without diagonal */
2399:       ncols = ai[i+1] - a->diag[i] -1;
2400:       for (j=0; j<ncols; j++) *cols++ = aj[j];
2401:       *cols++ = i; /* diagonal is located as the last entry of U(i,:) */
2402:     }
2403:   } else { /* case: levels>0 || (levels=0 && !perm_identity) */
2404:     ISGetIndices(iperm,&riip);
2405:     ISGetIndices(perm,&rip);

2407:     /* initialization */
2408:     PetscMalloc1(am+1,&ajtmp);

2410:     /* jl: linked list for storing indices of the pivot rows
2411:        il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2412:     PetscMalloc4(am,&uj_ptr,am,&uj_lvl_ptr,am,&jl,am,&il);
2413:     for (i=0; i<am; i++) {
2414:       jl[i] = am; il[i] = 0;
2415:     }

2417:     /* create and initialize a linked list for storing column indices of the active row k */
2418:     nlnk = am + 1;
2419:     PetscIncompleteLLCreate(am,am,nlnk,lnk,lnk_lvl,lnkbt);

2421:     /* initial FreeSpace size is fill*(ai[am]+am)/2 */
2422:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,(ai[am]+am)/2),&free_space);
2423:     current_space     = free_space;
2424:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,(ai[am]+am)/2),&free_space_lvl);
2425:     current_space_lvl = free_space_lvl;

2427:     for (k=0; k<am; k++) {  /* for each active row k */
2428:       /* initialize lnk by the column indices of row rip[k] of A */
2429:       nzk   = 0;
2430:       ncols = ai[rip[k]+1] - ai[rip[k]];
2431:       if (!ncols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",rip[k],k);
2432:       ncols_upper = 0;
2433:       for (j=0; j<ncols; j++) {
2434:         i = *(aj + ai[rip[k]] + j); /* unpermuted column index */
2435:         if (riip[i] >= k) { /* only take upper triangular entry */
2436:           ajtmp[ncols_upper] = i;
2437:           ncols_upper++;
2438:         }
2439:       }
2440:       PetscIncompleteLLInit(ncols_upper,ajtmp,am,riip,nlnk,lnk,lnk_lvl,lnkbt);
2441:       nzk += nlnk;

2443:       /* update lnk by computing fill-in for each pivot row to be merged in */
2444:       prow = jl[k]; /* 1st pivot row */

2446:       while (prow < k) {
2447:         nextprow = jl[prow];

2449:         /* merge prow into k-th row */
2450:         jmin  = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
2451:         jmax  = ui[prow+1];
2452:         ncols = jmax-jmin;
2453:         i     = jmin - ui[prow];
2454:         cols  = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
2455:         uj    = uj_lvl_ptr[prow] + i; /* levels of cols */
2456:         j     = *(uj - 1);
2457:         PetscICCLLAddSorted(ncols,cols,levels,uj,am,nlnk,lnk,lnk_lvl,lnkbt,j);
2458:         nzk  += nlnk;

2460:         /* update il and jl for prow */
2461:         if (jmin < jmax) {
2462:           il[prow] = jmin;
2463:           j        = *cols; jl[prow] = jl[j]; jl[j] = prow;
2464:         }
2465:         prow = nextprow;
2466:       }

2468:       /* if free space is not available, make more free space */
2469:       if (current_space->local_remaining<nzk) {
2470:         i    = am - k + 1; /* num of unfactored rows */
2471:         i    = PetscIntMultTruncate(i,PetscMin(nzk, i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
2472:         PetscFreeSpaceGet(i,&current_space);
2473:         PetscFreeSpaceGet(i,&current_space_lvl);
2474:         reallocs++;
2475:       }

2477:       /* copy data into free_space and free_space_lvl, then initialize lnk */
2478:       if (nzk == 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Empty row %D in ICC matrix factor",k);
2479:       PetscIncompleteLLClean(am,am,nzk,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);

2481:       /* add the k-th row into il and jl */
2482:       if (nzk > 1) {
2483:         i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
2484:         jl[k] = jl[i]; jl[i] = k;
2485:         il[k] = ui[k] + 1;
2486:       }
2487:       uj_ptr[k]     = current_space->array;
2488:       uj_lvl_ptr[k] = current_space_lvl->array;

2490:       current_space->array           += nzk;
2491:       current_space->local_used      += nzk;
2492:       current_space->local_remaining -= nzk;

2494:       current_space_lvl->array           += nzk;
2495:       current_space_lvl->local_used      += nzk;
2496:       current_space_lvl->local_remaining -= nzk;

2498:       ui[k+1] = ui[k] + nzk;
2499:     }

2501:     ISRestoreIndices(perm,&rip);
2502:     ISRestoreIndices(iperm,&riip);
2503:     PetscFree4(uj_ptr,uj_lvl_ptr,jl,il);
2504:     PetscFree(ajtmp);

2506:     /* copy free_space into uj and free free_space; set ui, uj, udiag in new datastructure; */
2507:     PetscMalloc1(ui[am]+1,&uj);
2508:     PetscFreeSpaceContiguous_Cholesky(&free_space,uj,am,ui,udiag); /* store matrix factor  */
2509:     PetscIncompleteLLDestroy(lnk,lnkbt);
2510:     PetscFreeSpaceDestroy(free_space_lvl);

2512:   } /* end of case: levels>0 || (levels=0 && !perm_identity) */

2514:   /* put together the new matrix in MATSEQSBAIJ format */
2515:   b               = (Mat_SeqSBAIJ*)(fact)->data;
2516:   b->singlemalloc = PETSC_FALSE;

2518:   PetscMalloc1(ui[am]+1,&b->a);

2520:   b->j             = uj;
2521:   b->i             = ui;
2522:   b->diag          = udiag;
2523:   b->free_diag     = PETSC_TRUE;
2524:   b->ilen          = NULL;
2525:   b->imax          = NULL;
2526:   b->row           = perm;
2527:   b->col           = perm;
2528:   PetscObjectReference((PetscObject)perm);
2529:   PetscObjectReference((PetscObject)perm);
2530:   b->icol          = iperm;
2531:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */

2533:   PetscMalloc1(am+1,&b->solve_work);
2534:   PetscLogObjectMemory((PetscObject)fact,ui[am]*(sizeof(PetscInt)+sizeof(MatScalar)));

2536:   b->maxnz   = b->nz = ui[am];
2537:   b->free_a  = PETSC_TRUE;
2538:   b->free_ij = PETSC_TRUE;

2540:   fact->info.factor_mallocs   = reallocs;
2541:   fact->info.fill_ratio_given = fill;
2542:   if (ai[am] != 0) {
2543:     /* nonzeros in lower triangular part of A (including diagonals) = (ai[am]+am)/2 */
2544:     fact->info.fill_ratio_needed = ((PetscReal)2*ui[am])/(ai[am]+am);
2545:   } else {
2546:     fact->info.fill_ratio_needed = 0.0;
2547:   }
2548: #if defined(PETSC_USE_INFO)
2549:   if (ai[am] != 0) {
2550:     PetscReal af = fact->info.fill_ratio_needed;
2551:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
2552:     PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
2553:     PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
2554:   } else {
2555:     PetscInfo(A,"Empty matrix\n");
2556:   }
2557: #endif
2558:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ;
2559:   return(0);
2560: }

2562: PetscErrorCode MatICCFactorSymbolic_SeqAIJ_inplace(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
2563: {
2564:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data;
2565:   Mat_SeqSBAIJ       *b;
2566:   PetscErrorCode     ierr;
2567:   PetscBool          perm_identity,missing;
2568:   PetscInt           reallocs=0,i,*ai=a->i,*aj=a->j,am=A->rmap->n,*ui,*udiag;
2569:   const PetscInt     *rip,*riip;
2570:   PetscInt           jmin,jmax,nzk,k,j,*jl,prow,*il,nextprow;
2571:   PetscInt           nlnk,*lnk,*lnk_lvl=NULL,d;
2572:   PetscInt           ncols,ncols_upper,*cols,*ajtmp,*uj,**uj_ptr,**uj_lvl_ptr;
2573:   PetscReal          fill          =info->fill,levels=info->levels;
2574:   PetscFreeSpaceList free_space    =NULL,current_space=NULL;
2575:   PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL;
2576:   PetscBT            lnkbt;
2577:   IS                 iperm;

2580:   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);
2581:   MatMissingDiagonal(A,&missing,&d);
2582:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
2583:   ISIdentity(perm,&perm_identity);
2584:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);

2586:   PetscMalloc1(am+1,&ui);
2587:   PetscMalloc1(am+1,&udiag);
2588:   ui[0] = 0;

2590:   /* ICC(0) without matrix ordering: simply copies fill pattern */
2591:   if (!levels && perm_identity) {

2593:     for (i=0; i<am; i++) {
2594:       ui[i+1]  = ui[i] + ai[i+1] - a->diag[i];
2595:       udiag[i] = ui[i];
2596:     }
2597:     PetscMalloc1(ui[am]+1,&uj);
2598:     cols = uj;
2599:     for (i=0; i<am; i++) {
2600:       aj    = a->j + a->diag[i];
2601:       ncols = ui[i+1] - ui[i];
2602:       for (j=0; j<ncols; j++) *cols++ = *aj++;
2603:     }
2604:   } else { /* case: levels>0 || (levels=0 && !perm_identity) */
2605:     ISGetIndices(iperm,&riip);
2606:     ISGetIndices(perm,&rip);

2608:     /* initialization */
2609:     PetscMalloc1(am+1,&ajtmp);

2611:     /* jl: linked list for storing indices of the pivot rows
2612:        il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2613:     PetscMalloc4(am,&uj_ptr,am,&uj_lvl_ptr,am,&jl,am,&il);
2614:     for (i=0; i<am; i++) {
2615:       jl[i] = am; il[i] = 0;
2616:     }

2618:     /* create and initialize a linked list for storing column indices of the active row k */
2619:     nlnk = am + 1;
2620:     PetscIncompleteLLCreate(am,am,nlnk,lnk,lnk_lvl,lnkbt);

2622:     /* initial FreeSpace size is fill*(ai[am]+1) */
2623:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,ai[am]+1),&free_space);
2624:     current_space     = free_space;
2625:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,ai[am]+1),&free_space_lvl);
2626:     current_space_lvl = free_space_lvl;

2628:     for (k=0; k<am; k++) {  /* for each active row k */
2629:       /* initialize lnk by the column indices of row rip[k] of A */
2630:       nzk   = 0;
2631:       ncols = ai[rip[k]+1] - ai[rip[k]];
2632:       if (!ncols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",rip[k],k);
2633:       ncols_upper = 0;
2634:       for (j=0; j<ncols; j++) {
2635:         i = *(aj + ai[rip[k]] + j); /* unpermuted column index */
2636:         if (riip[i] >= k) { /* only take upper triangular entry */
2637:           ajtmp[ncols_upper] = i;
2638:           ncols_upper++;
2639:         }
2640:       }
2641:       PetscIncompleteLLInit(ncols_upper,ajtmp,am,riip,nlnk,lnk,lnk_lvl,lnkbt);
2642:       nzk += nlnk;

2644:       /* update lnk by computing fill-in for each pivot row to be merged in */
2645:       prow = jl[k]; /* 1st pivot row */

2647:       while (prow < k) {
2648:         nextprow = jl[prow];

2650:         /* merge prow into k-th row */
2651:         jmin  = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
2652:         jmax  = ui[prow+1];
2653:         ncols = jmax-jmin;
2654:         i     = jmin - ui[prow];
2655:         cols  = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
2656:         uj    = uj_lvl_ptr[prow] + i; /* levels of cols */
2657:         j     = *(uj - 1);
2658:         PetscICCLLAddSorted(ncols,cols,levels,uj,am,nlnk,lnk,lnk_lvl,lnkbt,j);
2659:         nzk  += nlnk;

2661:         /* update il and jl for prow */
2662:         if (jmin < jmax) {
2663:           il[prow] = jmin;
2664:           j        = *cols; jl[prow] = jl[j]; jl[j] = prow;
2665:         }
2666:         prow = nextprow;
2667:       }

2669:       /* if free space is not available, make more free space */
2670:       if (current_space->local_remaining<nzk) {
2671:         i    = am - k + 1; /* num of unfactored rows */
2672:         i    = PetscIntMultTruncate(i,PetscMin(nzk, i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
2673:         PetscFreeSpaceGet(i,&current_space);
2674:         PetscFreeSpaceGet(i,&current_space_lvl);
2675:         reallocs++;
2676:       }

2678:       /* copy data into free_space and free_space_lvl, then initialize lnk */
2679:       if (!nzk) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Empty row %D in ICC matrix factor",k);
2680:       PetscIncompleteLLClean(am,am,nzk,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);

2682:       /* add the k-th row into il and jl */
2683:       if (nzk > 1) {
2684:         i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
2685:         jl[k] = jl[i]; jl[i] = k;
2686:         il[k] = ui[k] + 1;
2687:       }
2688:       uj_ptr[k]     = current_space->array;
2689:       uj_lvl_ptr[k] = current_space_lvl->array;

2691:       current_space->array           += nzk;
2692:       current_space->local_used      += nzk;
2693:       current_space->local_remaining -= nzk;

2695:       current_space_lvl->array           += nzk;
2696:       current_space_lvl->local_used      += nzk;
2697:       current_space_lvl->local_remaining -= nzk;

2699:       ui[k+1] = ui[k] + nzk;
2700:     }

2702: #if defined(PETSC_USE_INFO)
2703:     if (ai[am] != 0) {
2704:       PetscReal af = (PetscReal)ui[am]/((PetscReal)ai[am]);
2705:       PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
2706:       PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
2707:       PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
2708:     } else {
2709:       PetscInfo(A,"Empty matrix\n");
2710:     }
2711: #endif

2713:     ISRestoreIndices(perm,&rip);
2714:     ISRestoreIndices(iperm,&riip);
2715:     PetscFree4(uj_ptr,uj_lvl_ptr,jl,il);
2716:     PetscFree(ajtmp);

2718:     /* destroy list of free space and other temporary array(s) */
2719:     PetscMalloc1(ui[am]+1,&uj);
2720:     PetscFreeSpaceContiguous(&free_space,uj);
2721:     PetscIncompleteLLDestroy(lnk,lnkbt);
2722:     PetscFreeSpaceDestroy(free_space_lvl);

2724:   } /* end of case: levels>0 || (levels=0 && !perm_identity) */

2726:   /* put together the new matrix in MATSEQSBAIJ format */

2728:   b               = (Mat_SeqSBAIJ*)fact->data;
2729:   b->singlemalloc = PETSC_FALSE;

2731:   PetscMalloc1(ui[am]+1,&b->a);

2733:   b->j         = uj;
2734:   b->i         = ui;
2735:   b->diag      = udiag;
2736:   b->free_diag = PETSC_TRUE;
2737:   b->ilen      = NULL;
2738:   b->imax      = NULL;
2739:   b->row       = perm;
2740:   b->col       = perm;

2742:   PetscObjectReference((PetscObject)perm);
2743:   PetscObjectReference((PetscObject)perm);

2745:   b->icol          = iperm;
2746:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
2747:   PetscMalloc1(am+1,&b->solve_work);
2748:   PetscLogObjectMemory((PetscObject)fact,(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));
2749:   b->maxnz         = b->nz = ui[am];
2750:   b->free_a        = PETSC_TRUE;
2751:   b->free_ij       = PETSC_TRUE;

2753:   fact->info.factor_mallocs   = reallocs;
2754:   fact->info.fill_ratio_given = fill;
2755:   if (ai[am] != 0) {
2756:     fact->info.fill_ratio_needed = ((PetscReal)ui[am])/((PetscReal)ai[am]);
2757:   } else {
2758:     fact->info.fill_ratio_needed = 0.0;
2759:   }
2760:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_inplace;
2761:   return(0);
2762: }

2764: PetscErrorCode MatCholeskyFactorSymbolic_SeqAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
2765: {
2766:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data;
2767:   Mat_SeqSBAIJ       *b;
2768:   PetscErrorCode     ierr;
2769:   PetscBool          perm_identity,missing;
2770:   PetscReal          fill = info->fill;
2771:   const PetscInt     *rip,*riip;
2772:   PetscInt           i,am=A->rmap->n,*ai=a->i,*aj=a->j,reallocs=0,prow;
2773:   PetscInt           *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
2774:   PetscInt           nlnk,*lnk,ncols,ncols_upper,*cols,*uj,**ui_ptr,*uj_ptr,*udiag;
2775:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
2776:   PetscBT            lnkbt;
2777:   IS                 iperm;

2780:   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);
2781:   MatMissingDiagonal(A,&missing,&i);
2782:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);

2784:   /* check whether perm is the identity mapping */
2785:   ISIdentity(perm,&perm_identity);
2786:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);
2787:   ISGetIndices(iperm,&riip);
2788:   ISGetIndices(perm,&rip);

2790:   /* initialization */
2791:   PetscMalloc1(am+1,&ui);
2792:   PetscMalloc1(am+1,&udiag);
2793:   ui[0] = 0;

2795:   /* jl: linked list for storing indices of the pivot rows
2796:      il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2797:   PetscMalloc4(am,&ui_ptr,am,&jl,am,&il,am,&cols);
2798:   for (i=0; i<am; i++) {
2799:     jl[i] = am; il[i] = 0;
2800:   }

2802:   /* create and initialize a linked list for storing column indices of the active row k */
2803:   nlnk = am + 1;
2804:   PetscLLCreate(am,am,nlnk,lnk,lnkbt);

2806:   /* initial FreeSpace size is fill*(ai[am]+am)/2 */
2807:   PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,(ai[am]+am)/2),&free_space);
2808:   current_space = free_space;

2810:   for (k=0; k<am; k++) {  /* for each active row k */
2811:     /* initialize lnk by the column indices of row rip[k] of A */
2812:     nzk   = 0;
2813:     ncols = ai[rip[k]+1] - ai[rip[k]];
2814:     if (!ncols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",rip[k],k);
2815:     ncols_upper = 0;
2816:     for (j=0; j<ncols; j++) {
2817:       i = riip[*(aj + ai[rip[k]] + j)];
2818:       if (i >= k) { /* only take upper triangular entry */
2819:         cols[ncols_upper] = i;
2820:         ncols_upper++;
2821:       }
2822:     }
2823:     PetscLLAdd(ncols_upper,cols,am,nlnk,lnk,lnkbt);
2824:     nzk += nlnk;

2826:     /* update lnk by computing fill-in for each pivot row to be merged in */
2827:     prow = jl[k]; /* 1st pivot row */

2829:     while (prow < k) {
2830:       nextprow = jl[prow];
2831:       /* merge prow into k-th row */
2832:       jmin   = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
2833:       jmax   = ui[prow+1];
2834:       ncols  = jmax-jmin;
2835:       uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:am-1) */
2836:       PetscLLAddSorted(ncols,uj_ptr,am,nlnk,lnk,lnkbt);
2837:       nzk   += nlnk;

2839:       /* update il and jl for prow */
2840:       if (jmin < jmax) {
2841:         il[prow] = jmin;
2842:         j        = *uj_ptr;
2843:         jl[prow] = jl[j];
2844:         jl[j]    = prow;
2845:       }
2846:       prow = nextprow;
2847:     }

2849:     /* if free space is not available, make more free space */
2850:     if (current_space->local_remaining<nzk) {
2851:       i    = am - k + 1; /* num of unfactored rows */
2852:       i    = PetscIntMultTruncate(i,PetscMin(nzk,i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
2853:       PetscFreeSpaceGet(i,&current_space);
2854:       reallocs++;
2855:     }

2857:     /* copy data into free space, then initialize lnk */
2858:     PetscLLClean(am,am,nzk,lnk,current_space->array,lnkbt);

2860:     /* add the k-th row into il and jl */
2861:     if (nzk > 1) {
2862:       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
2863:       jl[k] = jl[i]; jl[i] = k;
2864:       il[k] = ui[k] + 1;
2865:     }
2866:     ui_ptr[k] = current_space->array;

2868:     current_space->array           += nzk;
2869:     current_space->local_used      += nzk;
2870:     current_space->local_remaining -= nzk;

2872:     ui[k+1] = ui[k] + nzk;
2873:   }

2875:   ISRestoreIndices(perm,&rip);
2876:   ISRestoreIndices(iperm,&riip);
2877:   PetscFree4(ui_ptr,jl,il,cols);

2879:   /* copy free_space into uj and free free_space; set ui, uj, udiag in new datastructure; */
2880:   PetscMalloc1(ui[am]+1,&uj);
2881:   PetscFreeSpaceContiguous_Cholesky(&free_space,uj,am,ui,udiag); /* store matrix factor */
2882:   PetscLLDestroy(lnk,lnkbt);

2884:   /* put together the new matrix in MATSEQSBAIJ format */

2886:   b               = (Mat_SeqSBAIJ*)fact->data;
2887:   b->singlemalloc = PETSC_FALSE;
2888:   b->free_a       = PETSC_TRUE;
2889:   b->free_ij      = PETSC_TRUE;

2891:   PetscMalloc1(ui[am]+1,&b->a);

2893:   b->j         = uj;
2894:   b->i         = ui;
2895:   b->diag      = udiag;
2896:   b->free_diag = PETSC_TRUE;
2897:   b->ilen      = NULL;
2898:   b->imax      = NULL;
2899:   b->row       = perm;
2900:   b->col       = perm;

2902:   PetscObjectReference((PetscObject)perm);
2903:   PetscObjectReference((PetscObject)perm);

2905:   b->icol          = iperm;
2906:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */

2908:   PetscMalloc1(am+1,&b->solve_work);
2909:   PetscLogObjectMemory((PetscObject)fact,ui[am]*(sizeof(PetscInt)+sizeof(MatScalar)));

2911:   b->maxnz = b->nz = ui[am];

2913:   fact->info.factor_mallocs   = reallocs;
2914:   fact->info.fill_ratio_given = fill;
2915:   if (ai[am] != 0) {
2916:     /* nonzeros in lower triangular part of A (including diagonals) = (ai[am]+am)/2 */
2917:     fact->info.fill_ratio_needed = ((PetscReal)2*ui[am])/(ai[am]+am);
2918:   } else {
2919:     fact->info.fill_ratio_needed = 0.0;
2920:   }
2921: #if defined(PETSC_USE_INFO)
2922:   if (ai[am] != 0) {
2923:     PetscReal af = fact->info.fill_ratio_needed;
2924:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
2925:     PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
2926:     PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
2927:   } else {
2928:     PetscInfo(A,"Empty matrix\n");
2929:   }
2930: #endif
2931:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ;
2932:   return(0);
2933: }

2935: PetscErrorCode MatCholeskyFactorSymbolic_SeqAIJ_inplace(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
2936: {
2937:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data;
2938:   Mat_SeqSBAIJ       *b;
2939:   PetscErrorCode     ierr;
2940:   PetscBool          perm_identity,missing;
2941:   PetscReal          fill = info->fill;
2942:   const PetscInt     *rip,*riip;
2943:   PetscInt           i,am=A->rmap->n,*ai=a->i,*aj=a->j,reallocs=0,prow;
2944:   PetscInt           *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
2945:   PetscInt           nlnk,*lnk,ncols,ncols_upper,*cols,*uj,**ui_ptr,*uj_ptr;
2946:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
2947:   PetscBT            lnkbt;
2948:   IS                 iperm;

2951:   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);
2952:   MatMissingDiagonal(A,&missing,&i);
2953:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);

2955:   /* check whether perm is the identity mapping */
2956:   ISIdentity(perm,&perm_identity);
2957:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);
2958:   ISGetIndices(iperm,&riip);
2959:   ISGetIndices(perm,&rip);

2961:   /* initialization */
2962:   PetscMalloc1(am+1,&ui);
2963:   ui[0] = 0;

2965:   /* jl: linked list for storing indices of the pivot rows
2966:      il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2967:   PetscMalloc4(am,&ui_ptr,am,&jl,am,&il,am,&cols);
2968:   for (i=0; i<am; i++) {
2969:     jl[i] = am; il[i] = 0;
2970:   }

2972:   /* create and initialize a linked list for storing column indices of the active row k */
2973:   nlnk = am + 1;
2974:   PetscLLCreate(am,am,nlnk,lnk,lnkbt);

2976:   /* initial FreeSpace size is fill*(ai[am]+1) */
2977:   PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,ai[am]+1),&free_space);
2978:   current_space = free_space;

2980:   for (k=0; k<am; k++) {  /* for each active row k */
2981:     /* initialize lnk by the column indices of row rip[k] of A */
2982:     nzk   = 0;
2983:     ncols = ai[rip[k]+1] - ai[rip[k]];
2984:     if (!ncols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",rip[k],k);
2985:     ncols_upper = 0;
2986:     for (j=0; j<ncols; j++) {
2987:       i = riip[*(aj + ai[rip[k]] + j)];
2988:       if (i >= k) { /* only take upper triangular entry */
2989:         cols[ncols_upper] = i;
2990:         ncols_upper++;
2991:       }
2992:     }
2993:     PetscLLAdd(ncols_upper,cols,am,nlnk,lnk,lnkbt);
2994:     nzk += nlnk;

2996:     /* update lnk by computing fill-in for each pivot row to be merged in */
2997:     prow = jl[k]; /* 1st pivot row */

2999:     while (prow < k) {
3000:       nextprow = jl[prow];
3001:       /* merge prow into k-th row */
3002:       jmin   = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
3003:       jmax   = ui[prow+1];
3004:       ncols  = jmax-jmin;
3005:       uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:am-1) */
3006:       PetscLLAddSorted(ncols,uj_ptr,am,nlnk,lnk,lnkbt);
3007:       nzk   += nlnk;

3009:       /* update il and jl for prow */
3010:       if (jmin < jmax) {
3011:         il[prow] = jmin;
3012:         j        = *uj_ptr; jl[prow] = jl[j]; jl[j] = prow;
3013:       }
3014:       prow = nextprow;
3015:     }

3017:     /* if free space is not available, make more free space */
3018:     if (current_space->local_remaining<nzk) {
3019:       i    = am - k + 1; /* num of unfactored rows */
3020:       i    = PetscMin(i*nzk, i*(i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
3021:       PetscFreeSpaceGet(i,&current_space);
3022:       reallocs++;
3023:     }

3025:     /* copy data into free space, then initialize lnk */
3026:     PetscLLClean(am,am,nzk,lnk,current_space->array,lnkbt);

3028:     /* add the k-th row into il and jl */
3029:     if (nzk-1 > 0) {
3030:       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
3031:       jl[k] = jl[i]; jl[i] = k;
3032:       il[k] = ui[k] + 1;
3033:     }
3034:     ui_ptr[k] = current_space->array;

3036:     current_space->array           += nzk;
3037:     current_space->local_used      += nzk;
3038:     current_space->local_remaining -= nzk;

3040:     ui[k+1] = ui[k] + nzk;
3041:   }

3043: #if defined(PETSC_USE_INFO)
3044:   if (ai[am] != 0) {
3045:     PetscReal af = (PetscReal)(ui[am])/((PetscReal)ai[am]);
3046:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
3047:     PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
3048:     PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
3049:   } else {
3050:     PetscInfo(A,"Empty matrix\n");
3051:   }
3052: #endif

3054:   ISRestoreIndices(perm,&rip);
3055:   ISRestoreIndices(iperm,&riip);
3056:   PetscFree4(ui_ptr,jl,il,cols);

3058:   /* destroy list of free space and other temporary array(s) */
3059:   PetscMalloc1(ui[am]+1,&uj);
3060:   PetscFreeSpaceContiguous(&free_space,uj);
3061:   PetscLLDestroy(lnk,lnkbt);

3063:   /* put together the new matrix in MATSEQSBAIJ format */

3065:   b               = (Mat_SeqSBAIJ*)fact->data;
3066:   b->singlemalloc = PETSC_FALSE;
3067:   b->free_a       = PETSC_TRUE;
3068:   b->free_ij      = PETSC_TRUE;

3070:   PetscMalloc1(ui[am]+1,&b->a);

3072:   b->j    = uj;
3073:   b->i    = ui;
3074:   b->diag = NULL;
3075:   b->ilen = NULL;
3076:   b->imax = NULL;
3077:   b->row  = perm;
3078:   b->col  = perm;

3080:   PetscObjectReference((PetscObject)perm);
3081:   PetscObjectReference((PetscObject)perm);

3083:   b->icol          = iperm;
3084:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */

3086:   PetscMalloc1(am+1,&b->solve_work);
3087:   PetscLogObjectMemory((PetscObject)fact,(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));
3088:   b->maxnz = b->nz = ui[am];

3090:   fact->info.factor_mallocs   = reallocs;
3091:   fact->info.fill_ratio_given = fill;
3092:   if (ai[am] != 0) {
3093:     fact->info.fill_ratio_needed = ((PetscReal)ui[am])/((PetscReal)ai[am]);
3094:   } else {
3095:     fact->info.fill_ratio_needed = 0.0;
3096:   }
3097:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_inplace;
3098:   return(0);
3099: }

3101: PetscErrorCode MatSolve_SeqAIJ_NaturalOrdering(Mat A,Vec bb,Vec xx)
3102: {
3103:   Mat_SeqAIJ        *a = (Mat_SeqAIJ*)A->data;
3104:   PetscErrorCode    ierr;
3105:   PetscInt          n   = A->rmap->n;
3106:   const PetscInt    *ai = a->i,*aj = a->j,*adiag = a->diag,*vi;
3107:   PetscScalar       *x,sum;
3108:   const PetscScalar *b;
3109:   const MatScalar   *aa = a->a,*v;
3110:   PetscInt          i,nz;

3113:   if (!n) return(0);

3115:   VecGetArrayRead(bb,&b);
3116:   VecGetArrayWrite(xx,&x);

3118:   /* forward solve the lower triangular */
3119:   x[0] = b[0];
3120:   v    = aa;
3121:   vi   = aj;
3122:   for (i=1; i<n; i++) {
3123:     nz  = ai[i+1] - ai[i];
3124:     sum = b[i];
3125:     PetscSparseDenseMinusDot(sum,x,v,vi,nz);
3126:     v   += nz;
3127:     vi  += nz;
3128:     x[i] = sum;
3129:   }

3131:   /* backward solve the upper triangular */
3132:   for (i=n-1; i>=0; i--) {
3133:     v   = aa + adiag[i+1] + 1;
3134:     vi  = aj + adiag[i+1] + 1;
3135:     nz  = adiag[i] - adiag[i+1]-1;
3136:     sum = x[i];
3137:     PetscSparseDenseMinusDot(sum,x,v,vi,nz);
3138:     x[i] = sum*v[nz]; /* x[i]=aa[adiag[i]]*sum; v++; */
3139:   }

3141:   PetscLogFlops(2.0*a->nz - A->cmap->n);
3142:   VecRestoreArrayRead(bb,&b);
3143:   VecRestoreArrayWrite(xx,&x);
3144:   return(0);
3145: }

3147: PetscErrorCode MatSolve_SeqAIJ(Mat A,Vec bb,Vec xx)
3148: {
3149:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
3150:   IS                iscol = a->col,isrow = a->row;
3151:   PetscErrorCode    ierr;
3152:   PetscInt          i,n=A->rmap->n,*vi,*ai=a->i,*aj=a->j,*adiag = a->diag,nz;
3153:   const PetscInt    *rout,*cout,*r,*c;
3154:   PetscScalar       *x,*tmp,sum;
3155:   const PetscScalar *b;
3156:   const MatScalar   *aa = a->a,*v;

3159:   if (!n) return(0);

3161:   VecGetArrayRead(bb,&b);
3162:   VecGetArrayWrite(xx,&x);
3163:   tmp  = a->solve_work;

3165:   ISGetIndices(isrow,&rout); r = rout;
3166:   ISGetIndices(iscol,&cout); c = cout;

3168:   /* forward solve the lower triangular */
3169:   tmp[0] = b[r[0]];
3170:   v      = aa;
3171:   vi     = aj;
3172:   for (i=1; i<n; i++) {
3173:     nz  = ai[i+1] - ai[i];
3174:     sum = b[r[i]];
3175:     PetscSparseDenseMinusDot(sum,tmp,v,vi,nz);
3176:     tmp[i] = sum;
3177:     v     += nz; vi += nz;
3178:   }

3180:   /* backward solve the upper triangular */
3181:   for (i=n-1; i>=0; i--) {
3182:     v   = aa + adiag[i+1]+1;
3183:     vi  = aj + adiag[i+1]+1;
3184:     nz  = adiag[i]-adiag[i+1]-1;
3185:     sum = tmp[i];
3186:     PetscSparseDenseMinusDot(sum,tmp,v,vi,nz);
3187:     x[c[i]] = tmp[i] = sum*v[nz]; /* v[nz] = aa[adiag[i]] */
3188:   }

3190:   ISRestoreIndices(isrow,&rout);
3191:   ISRestoreIndices(iscol,&cout);
3192:   VecRestoreArrayRead(bb,&b);
3193:   VecRestoreArrayWrite(xx,&x);
3194:   PetscLogFlops(2.0*a->nz - A->cmap->n);
3195:   return(0);
3196: }

3198: /*
3199:     This will get a new name and become a varient of MatILUFactor_SeqAIJ() there is no longer separate functions in the matrix function table for dt factors
3200: */
3201: PetscErrorCode MatILUDTFactor_SeqAIJ(Mat A,IS isrow,IS iscol,const MatFactorInfo *info,Mat *fact)
3202: {
3203:   Mat            B = *fact;
3204:   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data,*b;
3205:   IS             isicol;
3207:   const PetscInt *r,*ic;
3208:   PetscInt       i,n=A->rmap->n,*ai=a->i,*aj=a->j,*ajtmp,*adiag;
3209:   PetscInt       *bi,*bj,*bdiag,*bdiag_rev;
3210:   PetscInt       row,nzi,nzi_bl,nzi_bu,*im,nzi_al,nzi_au;
3211:   PetscInt       nlnk,*lnk;
3212:   PetscBT        lnkbt;
3213:   PetscBool      row_identity,icol_identity;
3214:   MatScalar      *aatmp,*pv,*batmp,*ba,*rtmp,*pc,multiplier,*vtmp,diag_tmp;
3215:   const PetscInt *ics;
3216:   PetscInt       j,nz,*pj,*bjtmp,k,ncut,*jtmp;
3217:   PetscReal      dt     =info->dt,shift=info->shiftamount;
3218:   PetscInt       dtcount=(PetscInt)info->dtcount,nnz_max;
3219:   PetscBool      missing;

3222:   if (dt      == PETSC_DEFAULT) dt = 0.005;
3223:   if (dtcount == PETSC_DEFAULT) dtcount = (PetscInt)(1.5*a->rmax);

3225:   /* ------- symbolic factorization, can be reused ---------*/
3226:   MatMissingDiagonal(A,&missing,&i);
3227:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);
3228:   adiag=a->diag;

3230:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);

3232:   /* bdiag is location of diagonal in factor */
3233:   PetscMalloc1(n+1,&bdiag);     /* becomes b->diag */
3234:   PetscMalloc1(n+1,&bdiag_rev); /* temporary */

3236:   /* allocate row pointers bi */
3237:   PetscMalloc1(2*n+2,&bi);

3239:   /* allocate bj and ba; max num of nonzero entries is (ai[n]+2*n*dtcount+2) */
3240:   if (dtcount > n-1) dtcount = n-1; /* diagonal is excluded */
3241:   nnz_max = ai[n]+2*n*dtcount+2;

3243:   PetscMalloc1(nnz_max+1,&bj);
3244:   PetscMalloc1(nnz_max+1,&ba);

3246:   /* put together the new matrix */
3247:   MatSeqAIJSetPreallocation_SeqAIJ(B,MAT_SKIP_ALLOCATION,NULL);
3248:   PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);
3249:   b    = (Mat_SeqAIJ*)B->data;

3251:   b->free_a       = PETSC_TRUE;
3252:   b->free_ij      = PETSC_TRUE;
3253:   b->singlemalloc = PETSC_FALSE;

3255:   b->a    = ba;
3256:   b->j    = bj;
3257:   b->i    = bi;
3258:   b->diag = bdiag;
3259:   b->ilen = NULL;
3260:   b->imax = NULL;
3261:   b->row  = isrow;
3262:   b->col  = iscol;
3263:   PetscObjectReference((PetscObject)isrow);
3264:   PetscObjectReference((PetscObject)iscol);
3265:   b->icol = isicol;

3267:   PetscMalloc1(n+1,&b->solve_work);
3268:   PetscLogObjectMemory((PetscObject)B,nnz_max*(sizeof(PetscInt)+sizeof(MatScalar)));
3269:   b->maxnz = nnz_max;

3271:   B->factortype            = MAT_FACTOR_ILUDT;
3272:   B->info.factor_mallocs   = 0;
3273:   B->info.fill_ratio_given = ((PetscReal)nnz_max)/((PetscReal)ai[n]);
3274:   /* ------- end of symbolic factorization ---------*/

3276:   ISGetIndices(isrow,&r);
3277:   ISGetIndices(isicol,&ic);
3278:   ics  = ic;

3280:   /* linked list for storing column indices of the active row */
3281:   nlnk = n + 1;
3282:   PetscLLCreate(n,n,nlnk,lnk,lnkbt);

3284:   /* im: used by PetscLLAddSortedLU(); jtmp: working array for column indices of active row */
3285:   PetscMalloc2(n,&im,n,&jtmp);
3286:   /* rtmp, vtmp: working arrays for sparse and contiguous row entries of active row */
3287:   PetscMalloc2(n,&rtmp,n,&vtmp);
3288:   PetscArrayzero(rtmp,n);

3290:   bi[0]        = 0;
3291:   bdiag[0]     = nnz_max-1; /* location of diag[0] in factor B */
3292:   bdiag_rev[n] = bdiag[0];
3293:   bi[2*n+1]    = bdiag[0]+1; /* endof bj and ba array */
3294:   for (i=0; i<n; i++) {
3295:     /* copy initial fill into linked list */
3296:     nzi = ai[r[i]+1] - ai[r[i]];
3297:     if (!nzi) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",r[i],i);
3298:     nzi_al = adiag[r[i]] - ai[r[i]];
3299:     nzi_au = ai[r[i]+1] - adiag[r[i]] -1;
3300:     ajtmp  = aj + ai[r[i]];
3301:     PetscLLAddPerm(nzi,ajtmp,ic,n,nlnk,lnk,lnkbt);

3303:     /* load in initial (unfactored row) */
3304:     aatmp = a->a + ai[r[i]];
3305:     for (j=0; j<nzi; j++) {
3306:       rtmp[ics[*ajtmp++]] = *aatmp++;
3307:     }

3309:     /* add pivot rows into linked list */
3310:     row = lnk[n];
3311:     while (row < i) {
3312:       nzi_bl = bi[row+1] - bi[row] + 1;
3313:       bjtmp  = bj + bdiag[row+1]+1; /* points to 1st column next to the diagonal in U */
3314:       PetscLLAddSortedLU(bjtmp,row,nlnk,lnk,lnkbt,i,nzi_bl,im);
3315:       nzi   += nlnk;
3316:       row    = lnk[row];
3317:     }

3319:     /* copy data from lnk into jtmp, then initialize lnk */
3320:     PetscLLClean(n,n,nzi,lnk,jtmp,lnkbt);

3322:     /* numerical factorization */
3323:     bjtmp = jtmp;
3324:     row   = *bjtmp++; /* 1st pivot row */
3325:     while (row < i) {
3326:       pc         = rtmp + row;
3327:       pv         = ba + bdiag[row]; /* 1./(diag of the pivot row) */
3328:       multiplier = (*pc) * (*pv);
3329:       *pc        = multiplier;
3330:       if (PetscAbsScalar(*pc) > dt) { /* apply tolerance dropping rule */
3331:         pj = bj + bdiag[row+1] + 1;         /* point to 1st entry of U(row,:) */
3332:         pv = ba + bdiag[row+1] + 1;
3333:         nz = bdiag[row] - bdiag[row+1] - 1;         /* num of entries in U(row,:), excluding diagonal */
3334:         for (j=0; j<nz; j++) rtmp[*pj++] -= multiplier * (*pv++);
3335:         PetscLogFlops(1+2.0*nz);
3336:       }
3337:       row = *bjtmp++;
3338:     }

3340:     /* copy sparse rtmp into contiguous vtmp; separate L and U part */
3341:     diag_tmp = rtmp[i];  /* save diagonal value - may not needed?? */
3342:     nzi_bl   = 0; j = 0;
3343:     while (jtmp[j] < i) { /* Note: jtmp is sorted */
3344:       vtmp[j] = rtmp[jtmp[j]]; rtmp[jtmp[j]]=0.0;
3345:       nzi_bl++; j++;
3346:     }
3347:     nzi_bu = nzi - nzi_bl -1;
3348:     while (j < nzi) {
3349:       vtmp[j] = rtmp[jtmp[j]]; rtmp[jtmp[j]]=0.0;
3350:       j++;
3351:     }

3353:     bjtmp = bj + bi[i];
3354:     batmp = ba + bi[i];
3355:     /* apply level dropping rule to L part */
3356:     ncut = nzi_al + dtcount;
3357:     if (ncut < nzi_bl) {
3358:       PetscSortSplit(ncut,nzi_bl,vtmp,jtmp);
3359:       PetscSortIntWithScalarArray(ncut,jtmp,vtmp);
3360:     } else {
3361:       ncut = nzi_bl;
3362:     }
3363:     for (j=0; j<ncut; j++) {
3364:       bjtmp[j] = jtmp[j];
3365:       batmp[j] = vtmp[j];
3366:     }
3367:     bi[i+1] = bi[i] + ncut;
3368:     nzi     = ncut + 1;

3370:     /* apply level dropping rule to U part */
3371:     ncut = nzi_au + dtcount;
3372:     if (ncut < nzi_bu) {
3373:       PetscSortSplit(ncut,nzi_bu,vtmp+nzi_bl+1,jtmp+nzi_bl+1);
3374:       PetscSortIntWithScalarArray(ncut,jtmp+nzi_bl+1,vtmp+nzi_bl+1);
3375:     } else {
3376:       ncut = nzi_bu;
3377:     }
3378:     nzi += ncut;

3380:     /* mark bdiagonal */
3381:     bdiag[i+1]       = bdiag[i] - (ncut + 1);
3382:     bdiag_rev[n-i-1] = bdiag[i+1];
3383:     bi[2*n - i]      = bi[2*n - i +1] - (ncut + 1);
3384:     bjtmp            = bj + bdiag[i];
3385:     batmp            = ba + bdiag[i];
3386:     *bjtmp           = i;
3387:     *batmp           = diag_tmp; /* rtmp[i]; */
3388:     if (*batmp == 0.0) {
3389:       *batmp = dt+shift;
3390:     }
3391:     *batmp = 1.0/(*batmp); /* invert diagonal entries for simplier triangular solves */

3393:     bjtmp = bj + bdiag[i+1]+1;
3394:     batmp = ba + bdiag[i+1]+1;
3395:     for (k=0; k<ncut; k++) {
3396:       bjtmp[k] = jtmp[nzi_bl+1+k];
3397:       batmp[k] = vtmp[nzi_bl+1+k];
3398:     }

3400:     im[i] = nzi;   /* used by PetscLLAddSortedLU() */
3401:   } /* for (i=0; i<n; i++) */
3402:   if (bi[n] >= bdiag[n]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"end of L array %d cannot >= the beginning of U array %d",bi[n],bdiag[n]);

3404:   ISRestoreIndices(isrow,&r);
3405:   ISRestoreIndices(isicol,&ic);

3407:   PetscLLDestroy(lnk,lnkbt);
3408:   PetscFree2(im,jtmp);
3409:   PetscFree2(rtmp,vtmp);
3410:   PetscFree(bdiag_rev);

3412:   PetscLogFlops(B->cmap->n);
3413:   b->maxnz = b->nz = bi[n] + bdiag[0] - bdiag[n];

3415:   ISIdentity(isrow,&row_identity);
3416:   ISIdentity(isicol,&icol_identity);
3417:   if (row_identity && icol_identity) {
3418:     B->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
3419:   } else {
3420:     B->ops->solve = MatSolve_SeqAIJ;
3421:   }

3423:   B->ops->solveadd          = NULL;
3424:   B->ops->solvetranspose    = NULL;
3425:   B->ops->solvetransposeadd = NULL;
3426:   B->ops->matsolve          = NULL;
3427:   B->assembled              = PETSC_TRUE;
3428:   B->preallocated           = PETSC_TRUE;
3429:   return(0);
3430: }

3432: /* a wraper of MatILUDTFactor_SeqAIJ() */
3433: /*
3434:     This will get a new name and become a varient of MatILUFactor_SeqAIJ() there is no longer separate functions in the matrix function table for dt factors
3435: */

3437: PetscErrorCode  MatILUDTFactorSymbolic_SeqAIJ(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
3438: {

3442:   MatILUDTFactor_SeqAIJ(A,row,col,info,&fact);
3443:   return(0);
3444: }

3446: /*
3447:    same as MatLUFactorNumeric_SeqAIJ(), except using contiguous array matrix factors
3448:    - intend to replace existing MatLUFactorNumeric_SeqAIJ()
3449: */
3450: /*
3451:     This will get a new name and become a varient of MatILUFactor_SeqAIJ() there is no longer separate functions in the matrix function table for dt factors
3452: */

3454: PetscErrorCode  MatILUDTFactorNumeric_SeqAIJ(Mat fact,Mat A,const MatFactorInfo *info)
3455: {
3456:   Mat            C     =fact;
3457:   Mat_SeqAIJ     *a    =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)C->data;
3458:   IS             isrow = b->row,isicol = b->icol;
3460:   const PetscInt *r,*ic,*ics;
3461:   PetscInt       i,j,k,n=A->rmap->n,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
3462:   PetscInt       *ajtmp,*bjtmp,nz,nzl,nzu,row,*bdiag = b->diag,*pj;
3463:   MatScalar      *rtmp,*pc,multiplier,*v,*pv,*aa=a->a;
3464:   PetscReal      dt=info->dt,shift=info->shiftamount;
3465:   PetscBool      row_identity, col_identity;

3468:   ISGetIndices(isrow,&r);
3469:   ISGetIndices(isicol,&ic);
3470:   PetscMalloc1(n+1,&rtmp);
3471:   ics  = ic;

3473:   for (i=0; i<n; i++) {
3474:     /* initialize rtmp array */
3475:     nzl   = bi[i+1] - bi[i];       /* num of nozeros in L(i,:) */
3476:     bjtmp = bj + bi[i];
3477:     for  (j=0; j<nzl; j++) rtmp[*bjtmp++] = 0.0;
3478:     rtmp[i] = 0.0;
3479:     nzu     = bdiag[i] - bdiag[i+1]; /* num of nozeros in U(i,:) */
3480:     bjtmp   = bj + bdiag[i+1] + 1;
3481:     for  (j=0; j<nzu; j++) rtmp[*bjtmp++] = 0.0;

3483:     /* load in initial unfactored row of A */
3484:     nz    = ai[r[i]+1] - ai[r[i]];
3485:     ajtmp = aj + ai[r[i]];
3486:     v     = aa + ai[r[i]];
3487:     for (j=0; j<nz; j++) {
3488:       rtmp[ics[*ajtmp++]] = v[j];
3489:     }

3491:     /* numerical factorization */
3492:     bjtmp = bj + bi[i]; /* point to 1st entry of L(i,:) */
3493:     nzl   = bi[i+1] - bi[i]; /* num of entries in L(i,:) */
3494:     k     = 0;
3495:     while (k < nzl) {
3496:       row        = *bjtmp++;
3497:       pc         = rtmp + row;
3498:       pv         = b->a + bdiag[row]; /* 1./(diag of the pivot row) */
3499:       multiplier = (*pc) * (*pv);
3500:       *pc        = multiplier;
3501:       if (PetscAbsScalar(multiplier) > dt) {
3502:         pj = bj + bdiag[row+1] + 1;         /* point to 1st entry of U(row,:) */
3503:         pv = b->a + bdiag[row+1] + 1;
3504:         nz = bdiag[row] - bdiag[row+1] - 1;         /* num of entries in U(row,:), excluding diagonal */
3505:         for (j=0; j<nz; j++) rtmp[*pj++] -= multiplier * (*pv++);
3506:         PetscLogFlops(1+2.0*nz);
3507:       }
3508:       k++;
3509:     }

3511:     /* finished row so stick it into b->a */
3512:     /* L-part */
3513:     pv  = b->a + bi[i];
3514:     pj  = bj + bi[i];
3515:     nzl = bi[i+1] - bi[i];
3516:     for (j=0; j<nzl; j++) {
3517:       pv[j] = rtmp[pj[j]];
3518:     }

3520:     /* diagonal: invert diagonal entries for simplier triangular solves */
3521:     if (rtmp[i] == 0.0) rtmp[i] = dt+shift;
3522:     b->a[bdiag[i]] = 1.0/rtmp[i];

3524:     /* U-part */
3525:     pv  = b->a + bdiag[i+1] + 1;
3526:     pj  = bj + bdiag[i+1] + 1;
3527:     nzu = bdiag[i] - bdiag[i+1] - 1;
3528:     for (j=0; j<nzu; j++) {
3529:       pv[j] = rtmp[pj[j]];
3530:     }
3531:   }

3533:   PetscFree(rtmp);
3534:   ISRestoreIndices(isicol,&ic);
3535:   ISRestoreIndices(isrow,&r);

3537:   ISIdentity(isrow,&row_identity);
3538:   ISIdentity(isicol,&col_identity);
3539:   if (row_identity && col_identity) {
3540:     C->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
3541:   } else {
3542:     C->ops->solve = MatSolve_SeqAIJ;
3543:   }
3544:   C->ops->solveadd          = NULL;
3545:   C->ops->solvetranspose    = NULL;
3546:   C->ops->solvetransposeadd = NULL;
3547:   C->ops->matsolve          = NULL;
3548:   C->assembled              = PETSC_TRUE;
3549:   C->preallocated           = PETSC_TRUE;

3551:   PetscLogFlops(C->cmap->n);
3552:   return(0);
3553: }