Actual source code: aijfact.c

petsc-master 2019-10-20
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 && (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian 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:   return(0);
122: }

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

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

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

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

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

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

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

162:   /* initial FreeSpace size is f*(ai[n]+1) */
163:   f             = info->fill;
164:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
165:   current_space = free_space;

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

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

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

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

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

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

224:   ISRestoreIndices(isrow,&r);
225:   ISRestoreIndices(isicol,&ic);

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

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

238:   b->free_a       = PETSC_TRUE;
239:   b->free_ij      = PETSC_TRUE;
240:   b->singlemalloc = PETSC_FALSE;

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

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

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

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

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

291:   if (A->rmap->N != A->cmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"matrix must be square");
292:   MatMissingDiagonal(A,&missing,&i);
293:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);
294: 
295:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
296:   ISGetIndices(isrow,&r);
297:   ISGetIndices(isicol,&ic);

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

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

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

310:   /* initial FreeSpace size is f*(ai[n]+1) */
311:   f             = info->fill;
312:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
313:   current_space = free_space;

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

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

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

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

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

356:     bi_ptr[i]                       = current_space->array;
357:     current_space->array           += nzi;
358:     current_space->local_used      += nzi;
359:     current_space->local_remaining -= nzi;
360:   }

362:   ISRestoreIndices(isrow,&r);
363:   ISRestoreIndices(isicol,&ic);

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

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

376:   b->free_a       = PETSC_TRUE;
377:   b->free_ij      = PETSC_TRUE;
378:   b->singlemalloc = PETSC_FALSE;

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

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

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

398:   B->factortype            = MAT_FACTOR_LU;
399:   B->info.factor_mallocs   = reallocs;
400:   B->info.fill_ratio_given = f;

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

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

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

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

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

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

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

488:   ISGetIndices(isrow,&r);
489:   ISGetIndices(isicol,&ic);
490:   PetscMalloc1(n+1,&rtmp);
491:   ics  = ic;

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

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

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

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

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

532:           for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
533:           PetscLogFlops(1+2*nz);
534:         }
535:         row = *bjtmp++;
536:       }

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

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

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

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

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

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

582:   PetscFree(rtmp);
583:   ISRestoreIndices(isicol,&ic);
584:   ISRestoreIndices(isrow,&r);

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

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

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

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

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

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

656:   ISGetIndices(isrow,&r);
657:   ISGetIndices(isicol,&ic);
658:   PetscMalloc1(n+1,&rtmp);
659:   ics  = ic;

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

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

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

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

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

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

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

743:   C->assembled    = PETSC_TRUE;
744:   C->preallocated = PETSC_TRUE;

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

757:   MatSeqAIJCheckInode(C);
758:   return(0);
759: }

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

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

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

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

810:   ISGetIndices(isrow,&r);
811:   ISGetIndices(isicol,&ic);
812:   PetscMalloc1(n+1,&rtmp);
813:   PetscArrayzero(rtmp,n+1);
814:   ics  = ic;

816: #if defined(MV)
817:   sctx.shift_top      = 0.;
818:   sctx.nshift_max     = 0;
819:   sctx.shift_lo       = 0.;
820:   sctx.shift_hi       = 0.;
821:   sctx.shift_fraction = 0.;

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

841:   sctx.shift_amount = 0.;
842:   sctx.nshift       = 0;
843: #endif

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

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

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

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

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

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

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

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

915:   PetscFree(rtmp);
916:   ISRestoreIndices(isicol,&ic);
917:   ISRestoreIndices(isrow,&r);

919:   A->ops->solve             = MatSolve_SeqAIJ_InplaceWithPerm;
920:   A->ops->solveadd          = MatSolveAdd_SeqAIJ_inplace;
921:   A->ops->solvetranspose    = MatSolveTranspose_SeqAIJ_inplace;
922:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ_inplace;

924:   A->assembled    = PETSC_TRUE;
925:   A->preallocated = PETSC_TRUE;

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

938: /* ----------------------------------------------------------- */
939: PetscErrorCode MatLUFactor_SeqAIJ(Mat A,IS row,IS col,const MatFactorInfo *info)
940: {
942:   Mat            C;

945:   MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&C);
946:   MatLUFactorSymbolic(C,A,row,col,info);
947:   MatLUFactorNumeric(C,A,info);

949:   A->ops->solve          = C->ops->solve;
950:   A->ops->solvetranspose = C->ops->solvetranspose;

952:   MatHeaderMerge(A,&C);
953:   PetscLogObjectParent((PetscObject)A,(PetscObject)((Mat_SeqAIJ*)(A->data))->icol);
954:   return(0);
955: }
956: /* ----------------------------------------------------------- */


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

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

974:   VecGetArrayRead(bb,&b);
975:   VecGetArrayWrite(xx,&x);
976:   tmp  = a->solve_work;

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

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

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

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

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

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

1027:   PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&bisdense);
1028:   if (!bisdense) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"B matrix must be a SeqDense matrix");
1029:   if (X != B) {
1030:     PetscObjectTypeCompare((PetscObject)X,MATSEQDENSE,&xisdense);
1031:     if (!xisdense) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"X matrix must be a SeqDense matrix");
1032:   }

1034:   MatDenseGetArrayRead(B,&b);
1035:   MatDenseGetArray(X,&x);

1037:   tmp  = a->solve_work;
1038:   ISGetIndices(isrow,&rout); r = rout;
1039:   ISGetIndices(iscol,&cout); c = cout;

1041:   for (neq=0; neq<B->cmap->n; neq++) {
1042:     /* forward solve the lower triangular */
1043:     tmp[0] = b[r[0]];
1044:     tmps   = tmp;
1045:     for (i=1; i<n; i++) {
1046:       v   = aa + ai[i];
1047:       vi  = aj + ai[i];
1048:       nz  = a->diag[i] - ai[i];
1049:       sum = b[r[i]];
1050:       PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
1051:       tmp[i] = sum;
1052:     }
1053:     /* backward solve the upper triangular */
1054:     for (i=n-1; i>=0; i--) {
1055:       v   = aa + a->diag[i] + 1;
1056:       vi  = aj + a->diag[i] + 1;
1057:       nz  = ai[i+1] - a->diag[i] - 1;
1058:       sum = tmp[i];
1059:       PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
1060:       x[c[i]] = tmp[i] = sum*aa[a->diag[i]];
1061:     }

1063:     b += n;
1064:     x += n;
1065:   }
1066:   ISRestoreIndices(isrow,&rout);
1067:   ISRestoreIndices(iscol,&cout);
1068:   MatDenseRestoreArrayRead(B,&b);
1069:   MatDenseRestoreArray(X,&x);
1070:   PetscLogFlops(B->cmap->n*(2.0*a->nz - n));
1071:   return(0);
1072: }

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

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

1090:   PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&bisdense);
1091:   if (!bisdense) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"B matrix must be a SeqDense matrix");
1092:   if (X != B) {
1093:     PetscObjectTypeCompare((PetscObject)X,MATSEQDENSE,&xisdense);
1094:     if (!xisdense) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"X matrix must be a SeqDense matrix");
1095:   }

1097:   MatDenseGetArrayRead(B,&b);
1098:   MatDenseGetArray(X,&x);

1100:   tmp  = a->solve_work;
1101:   ISGetIndices(isrow,&rout); r = rout;
1102:   ISGetIndices(iscol,&cout); c = cout;

1104:   for (neq=0; neq<B->cmap->n; neq++) {
1105:     /* forward solve the lower triangular */
1106:     tmp[0] = b[r[0]];
1107:     v      = aa;
1108:     vi     = aj;
1109:     for (i=1; i<n; i++) {
1110:       nz  = ai[i+1] - ai[i];
1111:       sum = b[r[i]];
1112:       PetscSparseDenseMinusDot(sum,tmp,v,vi,nz);
1113:       tmp[i] = sum;
1114:       v     += nz; vi += nz;
1115:     }

1117:     /* backward solve the upper triangular */
1118:     for (i=n-1; i>=0; i--) {
1119:       v   = aa + adiag[i+1]+1;
1120:       vi  = aj + adiag[i+1]+1;
1121:       nz  = adiag[i]-adiag[i+1]-1;
1122:       sum = tmp[i];
1123:       PetscSparseDenseMinusDot(sum,tmp,v,vi,nz);
1124:       x[c[i]] = tmp[i] = sum*v[nz]; /* v[nz] = aa[adiag[i]] */
1125:     }

1127:     b += n;
1128:     x += n;
1129:   }
1130:   ISRestoreIndices(isrow,&rout);
1131:   ISRestoreIndices(iscol,&cout);
1132:   MatDenseRestoreArrayRead(B,&b);
1133:   MatDenseRestoreArray(X,&x);
1134:   PetscLogFlops(B->cmap->n*(2.0*a->nz - n));
1135:   return(0);
1136: }

1138: PetscErrorCode MatSolve_SeqAIJ_InplaceWithPerm(Mat A,Vec bb,Vec xx)
1139: {
1140:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1141:   IS                iscol = a->col,isrow = a->row;
1142:   PetscErrorCode    ierr;
1143:   const PetscInt    *r,*c,*rout,*cout;
1144:   PetscInt          i, n = A->rmap->n,*vi,*ai = a->i,*aj = a->j;
1145:   PetscInt          nz,row;
1146:   PetscScalar       *x,*tmp,*tmps,sum;
1147:   const PetscScalar *b;
1148:   const MatScalar   *aa = a->a,*v;

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

1153:   VecGetArrayRead(bb,&b);
1154:   VecGetArrayWrite(xx,&x);
1155:   tmp  = a->solve_work;

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

1160:   /* forward solve the lower triangular */
1161:   tmp[0] = b[*r++];
1162:   tmps   = tmp;
1163:   for (row=1; row<n; row++) {
1164:     i   = rout[row]; /* permuted row */
1165:     v   = aa + ai[i];
1166:     vi  = aj + ai[i];
1167:     nz  = a->diag[i] - ai[i];
1168:     sum = b[*r++];
1169:     PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
1170:     tmp[row] = sum;
1171:   }

1173:   /* backward solve the upper triangular */
1174:   for (row=n-1; row>=0; row--) {
1175:     i   = rout[row]; /* permuted row */
1176:     v   = aa + a->diag[i] + 1;
1177:     vi  = aj + a->diag[i] + 1;
1178:     nz  = ai[i+1] - a->diag[i] - 1;
1179:     sum = tmp[row];
1180:     PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
1181:     x[*c--] = tmp[row] = sum*aa[a->diag[i]];
1182:   }

1184:   ISRestoreIndices(isrow,&rout);
1185:   ISRestoreIndices(iscol,&cout);
1186:   VecRestoreArrayRead(bb,&b);
1187:   VecRestoreArrayWrite(xx,&x);
1188:   PetscLogFlops(2.0*a->nz - A->cmap->n);
1189:   return(0);
1190: }

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

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

1213:   VecGetArrayRead(bb,&b);
1214:   VecGetArrayWrite(xx,&x);

1216: #if defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
1217:   fortransolveaij_(&n,x,ai,aj,adiag,aa,b);
1218: #else
1219:   /* forward solve the lower triangular */
1220:   x[0] = b[0];
1221:   for (i=1; i<n; i++) {
1222:     ai_i = ai[i];
1223:     v    = aa + ai_i;
1224:     vi   = aj + ai_i;
1225:     nz   = adiag[i] - ai_i;
1226:     sum  = b[i];
1227:     PetscSparseDenseMinusDot(sum,x,v,vi,nz);
1228:     x[i] = sum;
1229:   }

1231:   /* backward solve the upper triangular */
1232:   for (i=n-1; i>=0; i--) {
1233:     adiag_i = adiag[i];
1234:     v       = aa + adiag_i + 1;
1235:     vi      = aj + adiag_i + 1;
1236:     nz      = ai[i+1] - adiag_i - 1;
1237:     sum     = x[i];
1238:     PetscSparseDenseMinusDot(sum,x,v,vi,nz);
1239:     x[i] = sum*aa[adiag_i];
1240:   }
1241: #endif
1242:   PetscLogFlops(2.0*a->nz - A->cmap->n);
1243:   VecRestoreArrayRead(bb,&b);
1244:   VecRestoreArrayWrite(xx,&x);
1245:   return(0);
1246: }

1248: PetscErrorCode MatSolveAdd_SeqAIJ_inplace(Mat A,Vec bb,Vec yy,Vec xx)
1249: {
1250:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1251:   IS                iscol = a->col,isrow = a->row;
1252:   PetscErrorCode    ierr;
1253:   PetscInt          i, n = A->rmap->n,j;
1254:   PetscInt          nz;
1255:   const PetscInt    *rout,*cout,*r,*c,*vi,*ai = a->i,*aj = a->j;
1256:   PetscScalar       *x,*tmp,sum;
1257:   const PetscScalar *b;
1258:   const MatScalar   *aa = a->a,*v;

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

1263:   VecGetArrayRead(bb,&b);
1264:   VecGetArray(xx,&x);
1265:   tmp  = a->solve_work;

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

1270:   /* forward solve the lower triangular */
1271:   tmp[0] = b[*r++];
1272:   for (i=1; i<n; i++) {
1273:     v   = aa + ai[i];
1274:     vi  = aj + ai[i];
1275:     nz  = a->diag[i] - ai[i];
1276:     sum = b[*r++];
1277:     for (j=0; j<nz; j++) sum -= v[j]*tmp[vi[j]];
1278:     tmp[i] = sum;
1279:   }

1281:   /* backward solve the upper triangular */
1282:   for (i=n-1; i>=0; i--) {
1283:     v   = aa + a->diag[i] + 1;
1284:     vi  = aj + a->diag[i] + 1;
1285:     nz  = ai[i+1] - a->diag[i] - 1;
1286:     sum = tmp[i];
1287:     for (j=0; j<nz; j++) sum -= v[j]*tmp[vi[j]];
1288:     tmp[i]   = sum*aa[a->diag[i]];
1289:     x[*c--] += tmp[i];
1290:   }

1292:   ISRestoreIndices(isrow,&rout);
1293:   ISRestoreIndices(iscol,&cout);
1294:   VecRestoreArrayRead(bb,&b);
1295:   VecRestoreArray(xx,&x);
1296:   PetscLogFlops(2.0*a->nz);
1297:   return(0);
1298: }

1300: PetscErrorCode MatSolveAdd_SeqAIJ(Mat A,Vec bb,Vec yy,Vec xx)
1301: {
1302:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1303:   IS                iscol = a->col,isrow = a->row;
1304:   PetscErrorCode    ierr;
1305:   PetscInt          i, n = A->rmap->n,j;
1306:   PetscInt          nz;
1307:   const PetscInt    *rout,*cout,*r,*c,*vi,*ai = a->i,*aj = a->j,*adiag = a->diag;
1308:   PetscScalar       *x,*tmp,sum;
1309:   const PetscScalar *b;
1310:   const MatScalar   *aa = a->a,*v;

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

1315:   VecGetArrayRead(bb,&b);
1316:   VecGetArray(xx,&x);
1317:   tmp  = a->solve_work;

1319:   ISGetIndices(isrow,&rout); r = rout;
1320:   ISGetIndices(iscol,&cout); c = cout;

1322:   /* forward solve the lower triangular */
1323:   tmp[0] = b[r[0]];
1324:   v      = aa;
1325:   vi     = aj;
1326:   for (i=1; i<n; i++) {
1327:     nz  = ai[i+1] - ai[i];
1328:     sum = b[r[i]];
1329:     for (j=0; j<nz; j++) sum -= v[j]*tmp[vi[j]];
1330:     tmp[i] = sum;
1331:     v     += nz;
1332:     vi    += nz;
1333:   }

1335:   /* backward solve the upper triangular */
1336:   v  = aa + adiag[n-1];
1337:   vi = aj + adiag[n-1];
1338:   for (i=n-1; i>=0; i--) {
1339:     nz  = adiag[i] - adiag[i+1] - 1;
1340:     sum = tmp[i];
1341:     for (j=0; j<nz; j++) sum -= v[j]*tmp[vi[j]];
1342:     tmp[i]   = sum*v[nz];
1343:     x[c[i]] += tmp[i];
1344:     v       += nz+1; vi += nz+1;
1345:   }

1347:   ISRestoreIndices(isrow,&rout);
1348:   ISRestoreIndices(iscol,&cout);
1349:   VecRestoreArrayRead(bb,&b);
1350:   VecRestoreArray(xx,&x);
1351:   PetscLogFlops(2.0*a->nz);
1352:   return(0);
1353: }

1355: PetscErrorCode MatSolveTranspose_SeqAIJ_inplace(Mat A,Vec bb,Vec xx)
1356: {
1357:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1358:   IS                iscol = a->col,isrow = a->row;
1359:   PetscErrorCode    ierr;
1360:   const PetscInt    *rout,*cout,*r,*c,*diag = a->diag,*ai = a->i,*aj = a->j,*vi;
1361:   PetscInt          i,n = A->rmap->n,j;
1362:   PetscInt          nz;
1363:   PetscScalar       *x,*tmp,s1;
1364:   const MatScalar   *aa = a->a,*v;
1365:   const PetscScalar *b;

1368:   VecGetArrayRead(bb,&b);
1369:   VecGetArrayWrite(xx,&x);
1370:   tmp  = a->solve_work;

1372:   ISGetIndices(isrow,&rout); r = rout;
1373:   ISGetIndices(iscol,&cout); c = cout;

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

1378:   /* forward solve the U^T */
1379:   for (i=0; i<n; i++) {
1380:     v   = aa + diag[i];
1381:     vi  = aj + diag[i] + 1;
1382:     nz  = ai[i+1] - diag[i] - 1;
1383:     s1  = tmp[i];
1384:     s1 *= (*v++);  /* multiply by inverse of diagonal entry */
1385:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1386:     tmp[i] = s1;
1387:   }

1389:   /* backward solve the L^T */
1390:   for (i=n-1; i>=0; i--) {
1391:     v  = aa + diag[i] - 1;
1392:     vi = aj + diag[i] - 1;
1393:     nz = diag[i] - ai[i];
1394:     s1 = tmp[i];
1395:     for (j=0; j>-nz; j--) tmp[vi[j]] -= s1*v[j];
1396:   }

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

1401:   ISRestoreIndices(isrow,&rout);
1402:   ISRestoreIndices(iscol,&cout);
1403:   VecRestoreArrayRead(bb,&b);
1404:   VecRestoreArrayWrite(xx,&x);

1406:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1407:   return(0);
1408: }

1410: PetscErrorCode MatSolveTranspose_SeqAIJ(Mat A,Vec bb,Vec xx)
1411: {
1412:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1413:   IS                iscol = a->col,isrow = a->row;
1414:   PetscErrorCode    ierr;
1415:   const PetscInt    *rout,*cout,*r,*c,*adiag = a->diag,*ai = a->i,*aj = a->j,*vi;
1416:   PetscInt          i,n = A->rmap->n,j;
1417:   PetscInt          nz;
1418:   PetscScalar       *x,*tmp,s1;
1419:   const MatScalar   *aa = a->a,*v;
1420:   const PetscScalar *b;

1423:   VecGetArrayRead(bb,&b);
1424:   VecGetArrayWrite(xx,&x);
1425:   tmp  = a->solve_work;

1427:   ISGetIndices(isrow,&rout); r = rout;
1428:   ISGetIndices(iscol,&cout); c = cout;

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

1433:   /* forward solve the U^T */
1434:   for (i=0; i<n; i++) {
1435:     v   = aa + adiag[i+1] + 1;
1436:     vi  = aj + adiag[i+1] + 1;
1437:     nz  = adiag[i] - adiag[i+1] - 1;
1438:     s1  = tmp[i];
1439:     s1 *= v[nz];  /* multiply by inverse of diagonal entry */
1440:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1441:     tmp[i] = s1;
1442:   }

1444:   /* backward solve the L^T */
1445:   for (i=n-1; i>=0; i--) {
1446:     v  = aa + ai[i];
1447:     vi = aj + ai[i];
1448:     nz = ai[i+1] - ai[i];
1449:     s1 = tmp[i];
1450:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1451:   }

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

1456:   ISRestoreIndices(isrow,&rout);
1457:   ISRestoreIndices(iscol,&cout);
1458:   VecRestoreArrayRead(bb,&b);
1459:   VecRestoreArrayWrite(xx,&x);

1461:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1462:   return(0);
1463: }

1465: PetscErrorCode MatSolveTransposeAdd_SeqAIJ_inplace(Mat A,Vec bb,Vec zz,Vec xx)
1466: {
1467:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1468:   IS                iscol = a->col,isrow = a->row;
1469:   PetscErrorCode    ierr;
1470:   const PetscInt    *rout,*cout,*r,*c,*diag = a->diag,*ai = a->i,*aj = a->j,*vi;
1471:   PetscInt          i,n = A->rmap->n,j;
1472:   PetscInt          nz;
1473:   PetscScalar       *x,*tmp,s1;
1474:   const MatScalar   *aa = a->a,*v;
1475:   const PetscScalar *b;

1478:   if (zz != xx) {VecCopy(zz,xx);}
1479:   VecGetArrayRead(bb,&b);
1480:   VecGetArray(xx,&x);
1481:   tmp  = a->solve_work;

1483:   ISGetIndices(isrow,&rout); r = rout;
1484:   ISGetIndices(iscol,&cout); c = cout;

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

1489:   /* forward solve the U^T */
1490:   for (i=0; i<n; i++) {
1491:     v   = aa + diag[i];
1492:     vi  = aj + diag[i] + 1;
1493:     nz  = ai[i+1] - diag[i] - 1;
1494:     s1  = tmp[i];
1495:     s1 *= (*v++);  /* multiply by inverse of diagonal entry */
1496:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1497:     tmp[i] = s1;
1498:   }

1500:   /* backward solve the L^T */
1501:   for (i=n-1; i>=0; i--) {
1502:     v  = aa + diag[i] - 1;
1503:     vi = aj + diag[i] - 1;
1504:     nz = diag[i] - ai[i];
1505:     s1 = tmp[i];
1506:     for (j=0; j>-nz; j--) tmp[vi[j]] -= s1*v[j];
1507:   }

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

1512:   ISRestoreIndices(isrow,&rout);
1513:   ISRestoreIndices(iscol,&cout);
1514:   VecRestoreArrayRead(bb,&b);
1515:   VecRestoreArray(xx,&x);

1517:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1518:   return(0);
1519: }

1521: PetscErrorCode MatSolveTransposeAdd_SeqAIJ(Mat A,Vec bb,Vec zz,Vec xx)
1522: {
1523:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1524:   IS                iscol = a->col,isrow = a->row;
1525:   PetscErrorCode    ierr;
1526:   const PetscInt    *rout,*cout,*r,*c,*adiag = a->diag,*ai = a->i,*aj = a->j,*vi;
1527:   PetscInt          i,n = A->rmap->n,j;
1528:   PetscInt          nz;
1529:   PetscScalar       *x,*tmp,s1;
1530:   const MatScalar   *aa = a->a,*v;
1531:   const PetscScalar *b;

1534:   if (zz != xx) {VecCopy(zz,xx);}
1535:   VecGetArrayRead(bb,&b);
1536:   VecGetArray(xx,&x);
1537:   tmp  = a->solve_work;

1539:   ISGetIndices(isrow,&rout); r = rout;
1540:   ISGetIndices(iscol,&cout); c = cout;

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

1545:   /* forward solve the U^T */
1546:   for (i=0; i<n; i++) {
1547:     v   = aa + adiag[i+1] + 1;
1548:     vi  = aj + adiag[i+1] + 1;
1549:     nz  = adiag[i] - adiag[i+1] - 1;
1550:     s1  = tmp[i];
1551:     s1 *= v[nz];  /* multiply by inverse of diagonal entry */
1552:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1553:     tmp[i] = s1;
1554:   }


1557:   /* backward solve the L^T */
1558:   for (i=n-1; i>=0; i--) {
1559:     v  = aa + ai[i];
1560:     vi = aj + ai[i];
1561:     nz = ai[i+1] - ai[i];
1562:     s1 = tmp[i];
1563:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1564:   }

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

1569:   ISRestoreIndices(isrow,&rout);
1570:   ISRestoreIndices(iscol,&cout);
1571:   VecRestoreArrayRead(bb,&b);
1572:   VecRestoreArray(xx,&x);

1574:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1575:   return(0);
1576: }

1578: /* ----------------------------------------------------------------*/

1580: /*
1581:    ilu() under revised new data structure.
1582:    Factored arrays bj and ba are stored as
1583:      L(0,:), L(1,:), ...,L(n-1,:),  U(n-1,:),...,U(i,:),U(i-1,:),...,U(0,:)

1585:    bi=fact->i is an array of size n+1, in which
1586:    bi+
1587:      bi[i]:  points to 1st entry of L(i,:),i=0,...,n-1
1588:      bi[n]:  points to L(n-1,n-1)+1

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

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

1606:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
1607:   MatDuplicateNoCreate_SeqAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_FALSE);
1608:   b    = (Mat_SeqAIJ*)(fact)->data;

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

1614:   b->singlemalloc = PETSC_TRUE;
1615:   if (!b->diag) {
1616:     PetscMalloc1(n+1,&b->diag);
1617:     PetscLogObjectMemory((PetscObject)fact,(n+1)*sizeof(PetscInt));
1618:   }
1619:   bdiag = b->diag;

1621:   if (n > 0) {
1622:     PetscArrayzero(b->a,ai[n]);
1623:   }

1625:   /* set bi and bj with new data structure */
1626:   bi = b->i;
1627:   bj = b->j;

1629:   /* L part */
1630:   bi[0] = 0;
1631:   for (i=0; i<n; i++) {
1632:     nz      = adiag[i] - ai[i];
1633:     bi[i+1] = bi[i] + nz;
1634:     aj      = a->j + ai[i];
1635:     for (j=0; j<nz; j++) {
1636:       /*   *bj = aj[j]; bj++; */
1637:       bj[k++] = aj[j];
1638:     }
1639:   }

1641:   /* U part */
1642:   bdiag[n] = bi[n]-1;
1643:   for (i=n-1; i>=0; i--) {
1644:     nz = ai[i+1] - adiag[i] - 1;
1645:     aj = a->j + adiag[i] + 1;
1646:     for (j=0; j<nz; j++) {
1647:       /*      *bj = aj[j]; bj++; */
1648:       bj[k++] = aj[j];
1649:     }
1650:     /* diag[i] */
1651:     /*    *bj = i; bj++; */
1652:     bj[k++]  = i;
1653:     bdiag[i] = bdiag[i+1] + nz + 1;
1654:   }

1656:   fact->factortype             = MAT_FACTOR_ILU;
1657:   fact->info.factor_mallocs    = 0;
1658:   fact->info.fill_ratio_given  = info->fill;
1659:   fact->info.fill_ratio_needed = 1.0;
1660:   fact->ops->lufactornumeric   = MatLUFactorNumeric_SeqAIJ;
1661:   MatSeqAIJCheckInode_FactorLU(fact);

1663:   b       = (Mat_SeqAIJ*)(fact)->data;
1664:   b->row  = isrow;
1665:   b->col  = iscol;
1666:   b->icol = isicol;
1667:   PetscMalloc1(fact->rmap->n+1,&b->solve_work);
1668:   PetscObjectReference((PetscObject)isrow);
1669:   PetscObjectReference((PetscObject)iscol);
1670:   return(0);
1671: }

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

1692:   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);
1693:   MatMissingDiagonal(A,&missing,&i);
1694:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);
1695: 
1696:   levels = (PetscInt)info->levels;
1697:   ISIdentity(isrow,&row_identity);
1698:   ISIdentity(iscol,&col_identity);
1699:   if (!levels && row_identity && col_identity) {
1700:     /* special case: ilu(0) with natural ordering */
1701:     MatILUFactorSymbolic_SeqAIJ_ilu0(fact,A,isrow,iscol,info);
1702:     if (a->inode.size) {
1703:       fact->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
1704:     }
1705:     return(0);
1706:   }

1708:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
1709:   ISGetIndices(isrow,&r);
1710:   ISGetIndices(isicol,&ic);

1712:   /* get new row and diagonal pointers, must be allocated separately because they will be given to the Mat_SeqAIJ and freed separately */
1713:   PetscMalloc1(n+1,&bi);
1714:   PetscMalloc1(n+1,&bdiag);
1715:   bi[0] = bdiag[0] = 0;
1716:   PetscMalloc2(n,&bj_ptr,n,&bjlvl_ptr);

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

1722:   /* initial FreeSpace size is f*(ai[n]+1) */
1723:   f                 = info->fill;
1724:   diagonal_fill     = (PetscInt)info->diagonal_fill;
1725:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
1726:   current_space     = free_space;
1727:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space_lvl);
1728:   current_space_lvl = free_space_lvl;
1729:   for (i=0; i<n; i++) {
1730:     nzi = 0;
1731:     /* copy current row into linked list */
1732:     nnz = ai[r[i]+1] - ai[r[i]];
1733:     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);
1734:     cols   = aj + ai[r[i]];
1735:     lnk[i] = -1; /* marker to indicate if diagonal exists */
1736:     PetscIncompleteLLInit(nnz,cols,n,ic,nlnk,lnk,lnk_lvl,lnkbt);
1737:     nzi   += nlnk;

1739:     /* make sure diagonal entry is included */
1740:     if (diagonal_fill && lnk[i] == -1) {
1741:       fm = n;
1742:       while (lnk[fm] < i) fm = lnk[fm];
1743:       lnk[i]     = lnk[fm]; /* insert diagonal into linked list */
1744:       lnk[fm]    = i;
1745:       lnk_lvl[i] = 0;
1746:       nzi++; dcount++;
1747:     }

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

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

1777:     /* make sure the active row i has diagonal entry */
1778:     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);

1780:     current_space->array               += nzi;
1781:     current_space->local_used          += nzi;
1782:     current_space->local_remaining     -= nzi;
1783:     current_space_lvl->array           += nzi;
1784:     current_space_lvl->local_used      += nzi;
1785:     current_space_lvl->local_remaining -= nzi;
1786:   }

1788:   ISRestoreIndices(isrow,&r);
1789:   ISRestoreIndices(isicol,&ic);
1790:   /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
1791:   PetscMalloc1(bi[n]+1,&bj);
1792:   PetscFreeSpaceContiguous_LU(&free_space,bj,n,bi,bdiag);

1794:   PetscIncompleteLLDestroy(lnk,lnkbt);
1795:   PetscFreeSpaceDestroy(free_space_lvl);
1796:   PetscFree2(bj_ptr,bjlvl_ptr);

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

1815:   b->free_a       = PETSC_TRUE;
1816:   b->free_ij      = PETSC_TRUE;
1817:   b->singlemalloc = PETSC_FALSE;

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

1821:   b->j    = bj;
1822:   b->i    = bi;
1823:   b->diag = bdiag;
1824:   b->ilen = 0;
1825:   b->imax = 0;
1826:   b->row  = isrow;
1827:   b->col  = iscol;
1828:   PetscObjectReference((PetscObject)isrow);
1829:   PetscObjectReference((PetscObject)iscol);
1830:   b->icol = isicol;

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

1838:   (fact)->info.factor_mallocs    = reallocs;
1839:   (fact)->info.fill_ratio_given  = f;
1840:   (fact)->info.fill_ratio_needed = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]);
1841:   (fact)->ops->lufactornumeric   = MatLUFactorNumeric_SeqAIJ;
1842:   if (a->inode.size) {
1843:     (fact)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
1844:   }
1845:   MatSeqAIJCheckInode_FactorLU(fact);
1846:   return(0);
1847: }

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

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

1873:   f             = info->fill;
1874:   levels        = (PetscInt)info->levels;
1875:   diagonal_fill = (PetscInt)info->diagonal_fill;

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

1879:   ISIdentity(isrow,&row_identity);
1880:   ISIdentity(iscol,&col_identity);
1881:   if (!levels && row_identity && col_identity) { /* special case: ilu(0) with natural ordering */
1882:     MatDuplicateNoCreate_SeqAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_TRUE);

1884:     (fact)->ops->lufactornumeric =  MatLUFactorNumeric_SeqAIJ_inplace;
1885:     if (a->inode.size) {
1886:       (fact)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode_inplace;
1887:     }
1888:     fact->factortype               = MAT_FACTOR_ILU;
1889:     (fact)->info.factor_mallocs    = 0;
1890:     (fact)->info.fill_ratio_given  = info->fill;
1891:     (fact)->info.fill_ratio_needed = 1.0;

1893:     b    = (Mat_SeqAIJ*)(fact)->data;
1894:     b->row  = isrow;
1895:     b->col  = iscol;
1896:     b->icol = isicol;
1897:     PetscMalloc1((fact)->rmap->n+1,&b->solve_work);
1898:     PetscObjectReference((PetscObject)isrow);
1899:     PetscObjectReference((PetscObject)iscol);
1900:     return(0);
1901:   }

1903:   ISGetIndices(isrow,&r);
1904:   ISGetIndices(isicol,&ic);

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

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

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

1917:   /* initial FreeSpace size is f*(ai[n]+1) */
1918:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
1919:   current_space     = free_space;
1920:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space_lvl);
1921:   current_space_lvl = free_space_lvl;

1923:   for (i=0; i<n; i++) {
1924:     nzi = 0;
1925:     /* copy current row into linked list */
1926:     nnz = ai[r[i]+1] - ai[r[i]];
1927:     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);
1928:     cols   = aj + ai[r[i]];
1929:     lnk[i] = -1; /* marker to indicate if diagonal exists */
1930:     PetscIncompleteLLInit(nnz,cols,n,ic,nlnk,lnk,lnk_lvl,lnkbt);
1931:     nzi   += nlnk;

1933:     /* make sure diagonal entry is included */
1934:     if (diagonal_fill && lnk[i] == -1) {
1935:       fm = n;
1936:       while (lnk[fm] < i) fm = lnk[fm];
1937:       lnk[i]     = lnk[fm]; /* insert diagonal into linked list */
1938:       lnk[fm]    = i;
1939:       lnk_lvl[i] = 0;
1940:       nzi++; dcount++;
1941:     }

1943:     /* add pivot rows into the active row */
1944:     nzbd = 0;
1945:     prow = lnk[n];
1946:     while (prow < i) {
1947:       nnz      = bdiag[prow];
1948:       cols     = bj_ptr[prow] + nnz + 1;
1949:       cols_lvl = bjlvl_ptr[prow] + nnz + 1;
1950:       nnz      = bi[prow+1] - bi[prow] - nnz - 1;
1951:       PetscILULLAddSorted(nnz,cols,levels,cols_lvl,prow,nlnk,lnk,lnk_lvl,lnkbt,prow);
1952:       nzi     += nlnk;
1953:       prow     = lnk[prow];
1954:       nzbd++;
1955:     }
1956:     bdiag[i] = nzbd;
1957:     bi[i+1]  = bi[i] + nzi;

1959:     /* if free space is not available, make more free space */
1960:     if (current_space->local_remaining<nzi) {
1961:       nnz  = PetscIntMultTruncate(nzi,n - i); /* estimated and max additional space needed */
1962:       PetscFreeSpaceGet(nnz,&current_space);
1963:       PetscFreeSpaceGet(nnz,&current_space_lvl);
1964:       reallocs++;
1965:     }

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

1972:     /* make sure the active row i has diagonal entry */
1973:     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);

1975:     current_space->array               += nzi;
1976:     current_space->local_used          += nzi;
1977:     current_space->local_remaining     -= nzi;
1978:     current_space_lvl->array           += nzi;
1979:     current_space_lvl->local_used      += nzi;
1980:     current_space_lvl->local_remaining -= nzi;
1981:   }

1983:   ISRestoreIndices(isrow,&r);
1984:   ISRestoreIndices(isicol,&ic);

1986:   /* destroy list of free space and other temporary arrays */
1987:   PetscMalloc1(bi[n]+1,&bj);
1988:   PetscFreeSpaceContiguous(&free_space,bj); /* copy free_space -> bj */
1989:   PetscIncompleteLLDestroy(lnk,lnkbt);
1990:   PetscFreeSpaceDestroy(free_space_lvl);
1991:   PetscFree2(bj_ptr,bjlvl_ptr);

1993: #if defined(PETSC_USE_INFO)
1994:   {
1995:     PetscReal af = ((PetscReal)bi[n])/((PetscReal)ai[n]);
1996:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);
1997:     PetscInfo1(A,"Run with -[sub_]pc_factor_fill %g or use \n",(double)af);
1998:     PetscInfo1(A,"PCFactorSetFill([sub]pc,%g);\n",(double)af);
1999:     PetscInfo(A,"for best performance.\n");
2000:     if (diagonal_fill) {
2001:       PetscInfo1(A,"Detected and replaced %D missing diagonals\n",dcount);
2002:     }
2003:   }
2004: #endif

2006:   /* put together the new matrix */
2007:   MatSeqAIJSetPreallocation_SeqAIJ(fact,MAT_SKIP_ALLOCATION,NULL);
2008:   PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol);
2009:   b    = (Mat_SeqAIJ*)(fact)->data;

2011:   b->free_a       = PETSC_TRUE;
2012:   b->free_ij      = PETSC_TRUE;
2013:   b->singlemalloc = PETSC_FALSE;

2015:   PetscMalloc1(bi[n],&b->a);
2016:   b->j = bj;
2017:   b->i = bi;
2018:   for (i=0; i<n; i++) bdiag[i] += bi[i];
2019:   b->diag = bdiag;
2020:   b->ilen = 0;
2021:   b->imax = 0;
2022:   b->row  = isrow;
2023:   b->col  = iscol;
2024:   PetscObjectReference((PetscObject)isrow);
2025:   PetscObjectReference((PetscObject)iscol);
2026:   b->icol = isicol;
2027:   PetscMalloc1(n+1,&b->solve_work);
2028:   /* In b structure:  Free imax, ilen, old a, old j.
2029:      Allocate bdiag, solve_work, new a, new j */
2030:   PetscLogObjectMemory((PetscObject)fact,(bi[n]-n) * (sizeof(PetscInt)+sizeof(PetscScalar)));
2031:   b->maxnz = b->nz = bi[n];

2033:   (fact)->info.factor_mallocs    = reallocs;
2034:   (fact)->info.fill_ratio_given  = f;
2035:   (fact)->info.fill_ratio_needed = ((PetscReal)bi[n])/((PetscReal)ai[n]);
2036:   (fact)->ops->lufactornumeric   =  MatLUFactorNumeric_SeqAIJ_inplace;
2037:   if (a->inode.size) {
2038:     (fact)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode_inplace;
2039:   }
2040:   return(0);
2041: }

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

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

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

2081:   ISGetIndices(ip,&rip);
2082:   ISGetIndices(iip,&riip);

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

2090:   do {
2091:     sctx.newshift = PETSC_FALSE;

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

2096:     for (k = 0; k<mbs; k++) {
2097:       /* zero rtmp */
2098:       nz    = bi[k+1] - bi[k];
2099:       bjtmp = bj + bi[k];
2100:       for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

2102:       /* load in initial unfactored row */
2103:       bval = ba + bi[k];
2104:       jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
2105:       for (j = jmin; j < jmax; j++) {
2106:         col = riip[aj[j]];
2107:         if (col >= k) { /* only take upper triangular entry */
2108:           rtmp[col] = aa[j];
2109:           *bval++   = 0.0; /* for in-place factorization */
2110:         }
2111:       }
2112:       /* shift the diagonal of the matrix: ZeropivotApply() */
2113:       rtmp[k] += sctx.shift_amount;  /* shift the diagonal of the matrix */

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

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

2122:         /* compute multiplier, update diag(k) and U(i,k) */
2123:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
2124:         uikdi   = -ba[ili]*ba[bdiag[i]]; /* diagonal(k) */
2125:         dk     += uikdi*ba[ili]; /* update diag[k] */
2126:         ba[ili] = uikdi; /* -U(i,k) */

2128:         /* add multiple of row i to k-th row */
2129:         jmin = ili + 1; jmax = bi[i+1];
2130:         if (jmin < jmax) {
2131:           for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
2132:           /* update il and c2r for row i */
2133:           il[i] = jmin;
2134:           j     = bj[jmin]; c2r[i] = c2r[j]; c2r[j] = i;
2135:         }
2136:         i = nexti;
2137:       }

2139:       /* copy data into U(k,:) */
2140:       rs   = 0.0;
2141:       jmin = bi[k]; jmax = bi[k+1]-1;
2142:       if (jmin < jmax) {
2143:         for (j=jmin; j<jmax; j++) {
2144:           col = bj[j]; ba[j] = rtmp[col]; rs += PetscAbsScalar(ba[j]);
2145:         }
2146:         /* add the k-th row into il and c2r */
2147:         il[k] = jmin;
2148:         i     = bj[jmin]; c2r[k] = c2r[i]; c2r[i] = k;
2149:       }

2151:       /* MatPivotCheck() */
2152:       sctx.rs = rs;
2153:       sctx.pv = dk;
2154:       MatPivotCheck(B,A,info,&sctx,i);
2155:       if (sctx.newshift) break;
2156:       dk = sctx.pv;

2158:       ba[bdiag[k]] = 1.0/dk; /* U(k,k) */
2159:     }
2160:   } while (sctx.newshift);

2162:   PetscFree3(rtmp,il,c2r);
2163:   ISRestoreIndices(ip,&rip);
2164:   ISRestoreIndices(iip,&riip);

2166:   ISIdentity(ip,&perm_identity);
2167:   if (perm_identity) {
2168:     B->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering;
2169:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
2170:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
2171:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;
2172:   } else {
2173:     B->ops->solve          = MatSolve_SeqSBAIJ_1;
2174:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1;
2175:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1;
2176:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1;
2177:   }

2179:   C->assembled    = PETSC_TRUE;
2180:   C->preallocated = PETSC_TRUE;

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

2184:   /* MatPivotView() */
2185:   if (sctx.nshift) {
2186:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
2187:       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);
2188:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
2189:       PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
2190:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
2191:       PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %g\n",sctx.nshift,(double)info->shiftamount);
2192:     }
2193:   }
2194:   return(0);
2195: }

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

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

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

2235:   ISGetIndices(ip,&rip);
2236:   ISGetIndices(iip,&riip);

2238:   /* initialization */
2239:   PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&jl);

2241:   do {
2242:     sctx.newshift = PETSC_FALSE;

2244:     for (i=0; i<mbs; i++) jl[i] = mbs;
2245:     il[0] = 0;

2247:     for (k = 0; k<mbs; k++) {
2248:       /* zero rtmp */
2249:       nz    = bi[k+1] - bi[k];
2250:       bjtmp = bj + bi[k];
2251:       for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

2253:       bval = ba + bi[k];
2254:       /* initialize k-th row by the perm[k]-th row of A */
2255:       jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
2256:       for (j = jmin; j < jmax; j++) {
2257:         col = riip[aj[j]];
2258:         if (col >= k) { /* only take upper triangular entry */
2259:           rtmp[col] = aa[j];
2260:           *bval++   = 0.0; /* for in-place factorization */
2261:         }
2262:       }
2263:       /* shift the diagonal of the matrix */
2264:       if (sctx.nshift) rtmp[k] += sctx.shift_amount;

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

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

2273:         /* compute multiplier, update diag(k) and U(i,k) */
2274:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
2275:         uikdi   = -ba[ili]*ba[bi[i]]; /* diagonal(k) */
2276:         dk     += uikdi*ba[ili];
2277:         ba[ili] = uikdi; /* -U(i,k) */

2279:         /* add multiple of row i to k-th row */
2280:         jmin = ili + 1; jmax = bi[i+1];
2281:         if (jmin < jmax) {
2282:           for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
2283:           /* update il and jl for row i */
2284:           il[i] = jmin;
2285:           j     = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
2286:         }
2287:         i = nexti;
2288:       }

2290:       /* shift the diagonals when zero pivot is detected */
2291:       /* compute rs=sum of abs(off-diagonal) */
2292:       rs   = 0.0;
2293:       jmin = bi[k]+1;
2294:       nz   = bi[k+1] - jmin;
2295:       bcol = bj + jmin;
2296:       for (j=0; j<nz; j++) {
2297:         rs += PetscAbsScalar(rtmp[bcol[j]]);
2298:       }

2300:       sctx.rs = rs;
2301:       sctx.pv = dk;
2302:       MatPivotCheck(B,A,info,&sctx,k);
2303:       if (sctx.newshift) break;
2304:       dk = sctx.pv;

2306:       /* copy data into U(k,:) */
2307:       ba[bi[k]] = 1.0/dk; /* U(k,k) */
2308:       jmin      = bi[k]+1; jmax = bi[k+1];
2309:       if (jmin < jmax) {
2310:         for (j=jmin; j<jmax; j++) {
2311:           col = bj[j]; ba[j] = rtmp[col];
2312:         }
2313:         /* add the k-th row into il and jl */
2314:         il[k] = jmin;
2315:         i     = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
2316:       }
2317:     }
2318:   } while (sctx.newshift);

2320:   PetscFree3(rtmp,il,jl);
2321:   ISRestoreIndices(ip,&rip);
2322:   ISRestoreIndices(iip,&riip);

2324:   ISIdentity(ip,&perm_identity);
2325:   if (perm_identity) {
2326:     B->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2327:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2328:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2329:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2330:   } else {
2331:     B->ops->solve          = MatSolve_SeqSBAIJ_1_inplace;
2332:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
2333:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_inplace;
2334:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_inplace;
2335:   }

2337:   C->assembled    = PETSC_TRUE;
2338:   C->preallocated = PETSC_TRUE;

2340:   PetscLogFlops(C->rmap->n);
2341:   if (sctx.nshift) {
2342:     if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
2343:       PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
2344:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
2345:       PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
2346:     }
2347:   }
2348:   return(0);
2349: }

2351: /*
2352:    icc() under revised new data structure.
2353:    Factored arrays bj and ba are stored as
2354:      U(0,:),...,U(i,:),U(n-1,:)

2356:    ui=fact->i is an array of size n+1, in which
2357:    ui+
2358:      ui[i]:  points to 1st entry of U(i,:),i=0,...,n-1
2359:      ui[n]:  points to U(n-1,n-1)+1

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

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

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

2386:   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);
2387:   MatMissingDiagonal(A,&missing,&d);
2388:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
2389:   ISIdentity(perm,&perm_identity);
2390:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);

2392:   PetscMalloc1(am+1,&ui);
2393:   PetscMalloc1(am+1,&udiag);
2394:   ui[0] = 0;

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

2415:     /* initialization */
2416:     PetscMalloc1(am+1,&ajtmp);

2418:     /* jl: linked list for storing indices of the pivot rows
2419:        il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2420:     PetscMalloc4(am,&uj_ptr,am,&uj_lvl_ptr,am,&jl,am,&il);
2421:     for (i=0; i<am; i++) {
2422:       jl[i] = am; il[i] = 0;
2423:     }

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

2429:     /* initial FreeSpace size is fill*(ai[am]+am)/2 */
2430:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,(ai[am]+am)/2),&free_space);
2431:     current_space     = free_space;
2432:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,(ai[am]+am)/2),&free_space_lvl);
2433:     current_space_lvl = free_space_lvl;

2435:     for (k=0; k<am; k++) {  /* for each active row k */
2436:       /* initialize lnk by the column indices of row rip[k] of A */
2437:       nzk   = 0;
2438:       ncols = ai[rip[k]+1] - ai[rip[k]];
2439:       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);
2440:       ncols_upper = 0;
2441:       for (j=0; j<ncols; j++) {
2442:         i = *(aj + ai[rip[k]] + j); /* unpermuted column index */
2443:         if (riip[i] >= k) { /* only take upper triangular entry */
2444:           ajtmp[ncols_upper] = i;
2445:           ncols_upper++;
2446:         }
2447:       }
2448:       PetscIncompleteLLInit(ncols_upper,ajtmp,am,riip,nlnk,lnk,lnk_lvl,lnkbt);
2449:       nzk += nlnk;

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

2454:       while (prow < k) {
2455:         nextprow = jl[prow];

2457:         /* merge prow into k-th row */
2458:         jmin  = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
2459:         jmax  = ui[prow+1];
2460:         ncols = jmax-jmin;
2461:         i     = jmin - ui[prow];
2462:         cols  = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
2463:         uj    = uj_lvl_ptr[prow] + i; /* levels of cols */
2464:         j     = *(uj - 1);
2465:         PetscICCLLAddSorted(ncols,cols,levels,uj,am,nlnk,lnk,lnk_lvl,lnkbt,j);
2466:         nzk  += nlnk;

2468:         /* update il and jl for prow */
2469:         if (jmin < jmax) {
2470:           il[prow] = jmin;
2471:           j        = *cols; jl[prow] = jl[j]; jl[j] = prow;
2472:         }
2473:         prow = nextprow;
2474:       }

2476:       /* if free space is not available, make more free space */
2477:       if (current_space->local_remaining<nzk) {
2478:         i    = am - k + 1; /* num of unfactored rows */
2479:         i    = PetscIntMultTruncate(i,PetscMin(nzk, i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
2480:         PetscFreeSpaceGet(i,&current_space);
2481:         PetscFreeSpaceGet(i,&current_space_lvl);
2482:         reallocs++;
2483:       }

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

2489:       /* add the k-th row into il and jl */
2490:       if (nzk > 1) {
2491:         i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
2492:         jl[k] = jl[i]; jl[i] = k;
2493:         il[k] = ui[k] + 1;
2494:       }
2495:       uj_ptr[k]     = current_space->array;
2496:       uj_lvl_ptr[k] = current_space_lvl->array;

2498:       current_space->array           += nzk;
2499:       current_space->local_used      += nzk;
2500:       current_space->local_remaining -= nzk;

2502:       current_space_lvl->array           += nzk;
2503:       current_space_lvl->local_used      += nzk;
2504:       current_space_lvl->local_remaining -= nzk;

2506:       ui[k+1] = ui[k] + nzk;
2507:     }

2509:     ISRestoreIndices(perm,&rip);
2510:     ISRestoreIndices(iperm,&riip);
2511:     PetscFree4(uj_ptr,uj_lvl_ptr,jl,il);
2512:     PetscFree(ajtmp);

2514:     /* copy free_space into uj and free free_space; set ui, uj, udiag in new datastructure; */
2515:     PetscMalloc1(ui[am]+1,&uj);
2516:     PetscFreeSpaceContiguous_Cholesky(&free_space,uj,am,ui,udiag); /* store matrix factor  */
2517:     PetscIncompleteLLDestroy(lnk,lnkbt);
2518:     PetscFreeSpaceDestroy(free_space_lvl);

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

2522:   /* put together the new matrix in MATSEQSBAIJ format */
2523:   b               = (Mat_SeqSBAIJ*)(fact)->data;
2524:   b->singlemalloc = PETSC_FALSE;

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

2528:   b->j             = uj;
2529:   b->i             = ui;
2530:   b->diag          = udiag;
2531:   b->free_diag     = PETSC_TRUE;
2532:   b->ilen          = 0;
2533:   b->imax          = 0;
2534:   b->row           = perm;
2535:   b->col           = perm;
2536:   PetscObjectReference((PetscObject)perm);
2537:   PetscObjectReference((PetscObject)perm);
2538:   b->icol          = iperm;
2539:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */

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

2544:   b->maxnz   = b->nz = ui[am];
2545:   b->free_a  = PETSC_TRUE;
2546:   b->free_ij = PETSC_TRUE;

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

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

2588:   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);
2589:   MatMissingDiagonal(A,&missing,&d);
2590:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
2591:   ISIdentity(perm,&perm_identity);
2592:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);

2594:   PetscMalloc1(am+1,&ui);
2595:   PetscMalloc1(am+1,&udiag);
2596:   ui[0] = 0;

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

2601:     for (i=0; i<am; i++) {
2602:       ui[i+1]  = ui[i] + ai[i+1] - a->diag[i];
2603:       udiag[i] = ui[i];
2604:     }
2605:     PetscMalloc1(ui[am]+1,&uj);
2606:     cols = uj;
2607:     for (i=0; i<am; i++) {
2608:       aj    = a->j + a->diag[i];
2609:       ncols = ui[i+1] - ui[i];
2610:       for (j=0; j<ncols; j++) *cols++ = *aj++;
2611:     }
2612:   } else { /* case: levels>0 || (levels=0 && !perm_identity) */
2613:     ISGetIndices(iperm,&riip);
2614:     ISGetIndices(perm,&rip);

2616:     /* initialization */
2617:     PetscMalloc1(am+1,&ajtmp);

2619:     /* jl: linked list for storing indices of the pivot rows
2620:        il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2621:     PetscMalloc4(am,&uj_ptr,am,&uj_lvl_ptr,am,&jl,am,&il);
2622:     for (i=0; i<am; i++) {
2623:       jl[i] = am; il[i] = 0;
2624:     }

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

2630:     /* initial FreeSpace size is fill*(ai[am]+1) */
2631:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,ai[am]+1),&free_space);
2632:     current_space     = free_space;
2633:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,ai[am]+1),&free_space_lvl);
2634:     current_space_lvl = free_space_lvl;

2636:     for (k=0; k<am; k++) {  /* for each active row k */
2637:       /* initialize lnk by the column indices of row rip[k] of A */
2638:       nzk   = 0;
2639:       ncols = ai[rip[k]+1] - ai[rip[k]];
2640:       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);
2641:       ncols_upper = 0;
2642:       for (j=0; j<ncols; j++) {
2643:         i = *(aj + ai[rip[k]] + j); /* unpermuted column index */
2644:         if (riip[i] >= k) { /* only take upper triangular entry */
2645:           ajtmp[ncols_upper] = i;
2646:           ncols_upper++;
2647:         }
2648:       }
2649:       PetscIncompleteLLInit(ncols_upper,ajtmp,am,riip,nlnk,lnk,lnk_lvl,lnkbt);
2650:       nzk += nlnk;

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

2655:       while (prow < k) {
2656:         nextprow = jl[prow];

2658:         /* merge prow into k-th row */
2659:         jmin  = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
2660:         jmax  = ui[prow+1];
2661:         ncols = jmax-jmin;
2662:         i     = jmin - ui[prow];
2663:         cols  = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
2664:         uj    = uj_lvl_ptr[prow] + i; /* levels of cols */
2665:         j     = *(uj - 1);
2666:         PetscICCLLAddSorted(ncols,cols,levels,uj,am,nlnk,lnk,lnk_lvl,lnkbt,j);
2667:         nzk  += nlnk;

2669:         /* update il and jl for prow */
2670:         if (jmin < jmax) {
2671:           il[prow] = jmin;
2672:           j        = *cols; jl[prow] = jl[j]; jl[j] = prow;
2673:         }
2674:         prow = nextprow;
2675:       }

2677:       /* if free space is not available, make more free space */
2678:       if (current_space->local_remaining<nzk) {
2679:         i    = am - k + 1; /* num of unfactored rows */
2680:         i    = PetscIntMultTruncate(i,PetscMin(nzk, i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
2681:         PetscFreeSpaceGet(i,&current_space);
2682:         PetscFreeSpaceGet(i,&current_space_lvl);
2683:         reallocs++;
2684:       }

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

2690:       /* add the k-th row into il and jl */
2691:       if (nzk > 1) {
2692:         i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
2693:         jl[k] = jl[i]; jl[i] = k;
2694:         il[k] = ui[k] + 1;
2695:       }
2696:       uj_ptr[k]     = current_space->array;
2697:       uj_lvl_ptr[k] = current_space_lvl->array;

2699:       current_space->array           += nzk;
2700:       current_space->local_used      += nzk;
2701:       current_space->local_remaining -= nzk;

2703:       current_space_lvl->array           += nzk;
2704:       current_space_lvl->local_used      += nzk;
2705:       current_space_lvl->local_remaining -= nzk;

2707:       ui[k+1] = ui[k] + nzk;
2708:     }

2710: #if defined(PETSC_USE_INFO)
2711:     if (ai[am] != 0) {
2712:       PetscReal af = (PetscReal)ui[am]/((PetscReal)ai[am]);
2713:       PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
2714:       PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
2715:       PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
2716:     } else {
2717:       PetscInfo(A,"Empty matrix.\n");
2718:     }
2719: #endif

2721:     ISRestoreIndices(perm,&rip);
2722:     ISRestoreIndices(iperm,&riip);
2723:     PetscFree4(uj_ptr,uj_lvl_ptr,jl,il);
2724:     PetscFree(ajtmp);

2726:     /* destroy list of free space and other temporary array(s) */
2727:     PetscMalloc1(ui[am]+1,&uj);
2728:     PetscFreeSpaceContiguous(&free_space,uj);
2729:     PetscIncompleteLLDestroy(lnk,lnkbt);
2730:     PetscFreeSpaceDestroy(free_space_lvl);

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

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

2736:   b               = (Mat_SeqSBAIJ*)fact->data;
2737:   b->singlemalloc = PETSC_FALSE;

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

2741:   b->j         = uj;
2742:   b->i         = ui;
2743:   b->diag      = udiag;
2744:   b->free_diag = PETSC_TRUE;
2745:   b->ilen      = 0;
2746:   b->imax      = 0;
2747:   b->row       = perm;
2748:   b->col       = perm;

2750:   PetscObjectReference((PetscObject)perm);
2751:   PetscObjectReference((PetscObject)perm);

2753:   b->icol          = iperm;
2754:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
2755:   PetscMalloc1(am+1,&b->solve_work);
2756:   PetscLogObjectMemory((PetscObject)fact,(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));
2757:   b->maxnz         = b->nz = ui[am];
2758:   b->free_a        = PETSC_TRUE;
2759:   b->free_ij       = PETSC_TRUE;

2761:   fact->info.factor_mallocs   = reallocs;
2762:   fact->info.fill_ratio_given = fill;
2763:   if (ai[am] != 0) {
2764:     fact->info.fill_ratio_needed = ((PetscReal)ui[am])/((PetscReal)ai[am]);
2765:   } else {
2766:     fact->info.fill_ratio_needed = 0.0;
2767:   }
2768:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_inplace;
2769:   return(0);
2770: }

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

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

2792:   /* check whether perm is the identity mapping */
2793:   ISIdentity(perm,&perm_identity);
2794:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);
2795:   ISGetIndices(iperm,&riip);
2796:   ISGetIndices(perm,&rip);

2798:   /* initialization */
2799:   PetscMalloc1(am+1,&ui);
2800:   PetscMalloc1(am+1,&udiag);
2801:   ui[0] = 0;

2803:   /* jl: linked list for storing indices of the pivot rows
2804:      il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2805:   PetscMalloc4(am,&ui_ptr,am,&jl,am,&il,am,&cols);
2806:   for (i=0; i<am; i++) {
2807:     jl[i] = am; il[i] = 0;
2808:   }

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

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

2818:   for (k=0; k<am; k++) {  /* for each active row k */
2819:     /* initialize lnk by the column indices of row rip[k] of A */
2820:     nzk   = 0;
2821:     ncols = ai[rip[k]+1] - ai[rip[k]];
2822:     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);
2823:     ncols_upper = 0;
2824:     for (j=0; j<ncols; j++) {
2825:       i = riip[*(aj + ai[rip[k]] + j)];
2826:       if (i >= k) { /* only take upper triangular entry */
2827:         cols[ncols_upper] = i;
2828:         ncols_upper++;
2829:       }
2830:     }
2831:     PetscLLAdd(ncols_upper,cols,am,nlnk,lnk,lnkbt);
2832:     nzk += nlnk;

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

2837:     while (prow < k) {
2838:       nextprow = jl[prow];
2839:       /* merge prow into k-th row */
2840:       jmin   = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
2841:       jmax   = ui[prow+1];
2842:       ncols  = jmax-jmin;
2843:       uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:am-1) */
2844:       PetscLLAddSorted(ncols,uj_ptr,am,nlnk,lnk,lnkbt);
2845:       nzk   += nlnk;

2847:       /* update il and jl for prow */
2848:       if (jmin < jmax) {
2849:         il[prow] = jmin;
2850:         j        = *uj_ptr;
2851:         jl[prow] = jl[j];
2852:         jl[j]    = prow;
2853:       }
2854:       prow = nextprow;
2855:     }

2857:     /* if free space is not available, make more free space */
2858:     if (current_space->local_remaining<nzk) {
2859:       i    = am - k + 1; /* num of unfactored rows */
2860:       i    = PetscIntMultTruncate(i,PetscMin(nzk,i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
2861:       PetscFreeSpaceGet(i,&current_space);
2862:       reallocs++;
2863:     }

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

2868:     /* add the k-th row into il and jl */
2869:     if (nzk > 1) {
2870:       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
2871:       jl[k] = jl[i]; jl[i] = k;
2872:       il[k] = ui[k] + 1;
2873:     }
2874:     ui_ptr[k] = current_space->array;

2876:     current_space->array           += nzk;
2877:     current_space->local_used      += nzk;
2878:     current_space->local_remaining -= nzk;

2880:     ui[k+1] = ui[k] + nzk;
2881:   }

2883:   ISRestoreIndices(perm,&rip);
2884:   ISRestoreIndices(iperm,&riip);
2885:   PetscFree4(ui_ptr,jl,il,cols);

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

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

2894:   b               = (Mat_SeqSBAIJ*)fact->data;
2895:   b->singlemalloc = PETSC_FALSE;
2896:   b->free_a       = PETSC_TRUE;
2897:   b->free_ij      = PETSC_TRUE;

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

2901:   b->j         = uj;
2902:   b->i         = ui;
2903:   b->diag      = udiag;
2904:   b->free_diag = PETSC_TRUE;
2905:   b->ilen      = 0;
2906:   b->imax      = 0;
2907:   b->row       = perm;
2908:   b->col       = perm;

2910:   PetscObjectReference((PetscObject)perm);
2911:   PetscObjectReference((PetscObject)perm);

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

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

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

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

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

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

2963:   /* check whether perm is the identity mapping */
2964:   ISIdentity(perm,&perm_identity);
2965:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);
2966:   ISGetIndices(iperm,&riip);
2967:   ISGetIndices(perm,&rip);

2969:   /* initialization */
2970:   PetscMalloc1(am+1,&ui);
2971:   ui[0] = 0;

2973:   /* jl: linked list for storing indices of the pivot rows
2974:      il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
2975:   PetscMalloc4(am,&ui_ptr,am,&jl,am,&il,am,&cols);
2976:   for (i=0; i<am; i++) {
2977:     jl[i] = am; il[i] = 0;
2978:   }

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

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

2988:   for (k=0; k<am; k++) {  /* for each active row k */
2989:     /* initialize lnk by the column indices of row rip[k] of A */
2990:     nzk   = 0;
2991:     ncols = ai[rip[k]+1] - ai[rip[k]];
2992:     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);
2993:     ncols_upper = 0;
2994:     for (j=0; j<ncols; j++) {
2995:       i = riip[*(aj + ai[rip[k]] + j)];
2996:       if (i >= k) { /* only take upper triangular entry */
2997:         cols[ncols_upper] = i;
2998:         ncols_upper++;
2999:       }
3000:     }
3001:     PetscLLAdd(ncols_upper,cols,am,nlnk,lnk,lnkbt);
3002:     nzk += nlnk;

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

3007:     while (prow < k) {
3008:       nextprow = jl[prow];
3009:       /* merge prow into k-th row */
3010:       jmin   = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
3011:       jmax   = ui[prow+1];
3012:       ncols  = jmax-jmin;
3013:       uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:am-1) */
3014:       PetscLLAddSorted(ncols,uj_ptr,am,nlnk,lnk,lnkbt);
3015:       nzk   += nlnk;

3017:       /* update il and jl for prow */
3018:       if (jmin < jmax) {
3019:         il[prow] = jmin;
3020:         j        = *uj_ptr; jl[prow] = jl[j]; jl[j] = prow;
3021:       }
3022:       prow = nextprow;
3023:     }

3025:     /* if free space is not available, make more free space */
3026:     if (current_space->local_remaining<nzk) {
3027:       i    = am - k + 1; /* num of unfactored rows */
3028:       i    = PetscMin(i*nzk, i*(i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
3029:       PetscFreeSpaceGet(i,&current_space);
3030:       reallocs++;
3031:     }

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

3036:     /* add the k-th row into il and jl */
3037:     if (nzk-1 > 0) {
3038:       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
3039:       jl[k] = jl[i]; jl[i] = k;
3040:       il[k] = ui[k] + 1;
3041:     }
3042:     ui_ptr[k] = current_space->array;

3044:     current_space->array           += nzk;
3045:     current_space->local_used      += nzk;
3046:     current_space->local_remaining -= nzk;

3048:     ui[k+1] = ui[k] + nzk;
3049:   }

3051: #if defined(PETSC_USE_INFO)
3052:   if (ai[am] != 0) {
3053:     PetscReal af = (PetscReal)(ui[am])/((PetscReal)ai[am]);
3054:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);
3055:     PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
3056:     PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);
3057:   } else {
3058:     PetscInfo(A,"Empty matrix.\n");
3059:   }
3060: #endif

3062:   ISRestoreIndices(perm,&rip);
3063:   ISRestoreIndices(iperm,&riip);
3064:   PetscFree4(ui_ptr,jl,il,cols);

3066:   /* destroy list of free space and other temporary array(s) */
3067:   PetscMalloc1(ui[am]+1,&uj);
3068:   PetscFreeSpaceContiguous(&free_space,uj);
3069:   PetscLLDestroy(lnk,lnkbt);

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

3073:   b               = (Mat_SeqSBAIJ*)fact->data;
3074:   b->singlemalloc = PETSC_FALSE;
3075:   b->free_a       = PETSC_TRUE;
3076:   b->free_ij      = PETSC_TRUE;

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

3080:   b->j    = uj;
3081:   b->i    = ui;
3082:   b->diag = 0;
3083:   b->ilen = 0;
3084:   b->imax = 0;
3085:   b->row  = perm;
3086:   b->col  = perm;

3088:   PetscObjectReference((PetscObject)perm);
3089:   PetscObjectReference((PetscObject)perm);

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

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

3098:   fact->info.factor_mallocs   = reallocs;
3099:   fact->info.fill_ratio_given = fill;
3100:   if (ai[am] != 0) {
3101:     fact->info.fill_ratio_needed = ((PetscReal)ui[am])/((PetscReal)ai[am]);
3102:   } else {
3103:     fact->info.fill_ratio_needed = 0.0;
3104:   }
3105:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_inplace;
3106:   return(0);
3107: }

3109: PetscErrorCode MatSolve_SeqAIJ_NaturalOrdering(Mat A,Vec bb,Vec xx)
3110: {
3111:   Mat_SeqAIJ        *a = (Mat_SeqAIJ*)A->data;
3112:   PetscErrorCode    ierr;
3113:   PetscInt          n   = A->rmap->n;
3114:   const PetscInt    *ai = a->i,*aj = a->j,*adiag = a->diag,*vi;
3115:   PetscScalar       *x,sum;
3116:   const PetscScalar *b;
3117:   const MatScalar   *aa = a->a,*v;
3118:   PetscInt          i,nz;

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

3123:   VecGetArrayRead(bb,&b);
3124:   VecGetArrayWrite(xx,&x);

3126:   /* forward solve the lower triangular */
3127:   x[0] = b[0];
3128:   v    = aa;
3129:   vi   = aj;
3130:   for (i=1; i<n; i++) {
3131:     nz  = ai[i+1] - ai[i];
3132:     sum = b[i];
3133:     PetscSparseDenseMinusDot(sum,x,v,vi,nz);
3134:     v   += nz;
3135:     vi  += nz;
3136:     x[i] = sum;
3137:   }

3139:   /* backward solve the upper triangular */
3140:   for (i=n-1; i>=0; i--) {
3141:     v   = aa + adiag[i+1] + 1;
3142:     vi  = aj + adiag[i+1] + 1;
3143:     nz  = adiag[i] - adiag[i+1]-1;
3144:     sum = x[i];
3145:     PetscSparseDenseMinusDot(sum,x,v,vi,nz);
3146:     x[i] = sum*v[nz]; /* x[i]=aa[adiag[i]]*sum; v++; */
3147:   }

3149:   PetscLogFlops(2.0*a->nz - A->cmap->n);
3150:   VecRestoreArrayRead(bb,&b);
3151:   VecRestoreArrayWrite(xx,&x);
3152:   return(0);
3153: }

3155: PetscErrorCode MatSolve_SeqAIJ(Mat A,Vec bb,Vec xx)
3156: {
3157:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
3158:   IS                iscol = a->col,isrow = a->row;
3159:   PetscErrorCode    ierr;
3160:   PetscInt          i,n=A->rmap->n,*vi,*ai=a->i,*aj=a->j,*adiag = a->diag,nz;
3161:   const PetscInt    *rout,*cout,*r,*c;
3162:   PetscScalar       *x,*tmp,sum;
3163:   const PetscScalar *b;
3164:   const MatScalar   *aa = a->a,*v;

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

3169:   VecGetArrayRead(bb,&b);
3170:   VecGetArrayWrite(xx,&x);
3171:   tmp  = a->solve_work;

3173:   ISGetIndices(isrow,&rout); r = rout;
3174:   ISGetIndices(iscol,&cout); c = cout;

3176:   /* forward solve the lower triangular */
3177:   tmp[0] = b[r[0]];
3178:   v      = aa;
3179:   vi     = aj;
3180:   for (i=1; i<n; i++) {
3181:     nz  = ai[i+1] - ai[i];
3182:     sum = b[r[i]];
3183:     PetscSparseDenseMinusDot(sum,tmp,v,vi,nz);
3184:     tmp[i] = sum;
3185:     v     += nz; vi += nz;
3186:   }

3188:   /* backward solve the upper triangular */
3189:   for (i=n-1; i>=0; i--) {
3190:     v   = aa + adiag[i+1]+1;
3191:     vi  = aj + adiag[i+1]+1;
3192:     nz  = adiag[i]-adiag[i+1]-1;
3193:     sum = tmp[i];
3194:     PetscSparseDenseMinusDot(sum,tmp,v,vi,nz);
3195:     x[c[i]] = tmp[i] = sum*v[nz]; /* v[nz] = aa[adiag[i]] */
3196:   }

3198:   ISRestoreIndices(isrow,&rout);
3199:   ISRestoreIndices(iscol,&cout);
3200:   VecRestoreArrayRead(bb,&b);
3201:   VecRestoreArrayWrite(xx,&x);
3202:   PetscLogFlops(2*a->nz - A->cmap->n);
3203:   return(0);
3204: }

3206: /*
3207:     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
3208: */
3209: PetscErrorCode MatILUDTFactor_SeqAIJ(Mat A,IS isrow,IS iscol,const MatFactorInfo *info,Mat *fact)
3210: {
3211:   Mat            B = *fact;
3212:   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data,*b;
3213:   IS             isicol;
3215:   const PetscInt *r,*ic;
3216:   PetscInt       i,n=A->rmap->n,*ai=a->i,*aj=a->j,*ajtmp,*adiag;
3217:   PetscInt       *bi,*bj,*bdiag,*bdiag_rev;
3218:   PetscInt       row,nzi,nzi_bl,nzi_bu,*im,nzi_al,nzi_au;
3219:   PetscInt       nlnk,*lnk;
3220:   PetscBT        lnkbt;
3221:   PetscBool      row_identity,icol_identity;
3222:   MatScalar      *aatmp,*pv,*batmp,*ba,*rtmp,*pc,multiplier,*vtmp,diag_tmp;
3223:   const PetscInt *ics;
3224:   PetscInt       j,nz,*pj,*bjtmp,k,ncut,*jtmp;
3225:   PetscReal      dt     =info->dt,shift=info->shiftamount;
3226:   PetscInt       dtcount=(PetscInt)info->dtcount,nnz_max;
3227:   PetscBool      missing;

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

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

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

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

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

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

3251:   PetscMalloc1(nnz_max+1,&bj);
3252:   PetscMalloc1(nnz_max+1,&ba);

3254:   /* put together the new matrix */
3255:   MatSeqAIJSetPreallocation_SeqAIJ(B,MAT_SKIP_ALLOCATION,NULL);
3256:   PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);
3257:   b    = (Mat_SeqAIJ*)B->data;

3259:   b->free_a       = PETSC_TRUE;
3260:   b->free_ij      = PETSC_TRUE;
3261:   b->singlemalloc = PETSC_FALSE;

3263:   b->a    = ba;
3264:   b->j    = bj;
3265:   b->i    = bi;
3266:   b->diag = bdiag;
3267:   b->ilen = 0;
3268:   b->imax = 0;
3269:   b->row  = isrow;
3270:   b->col  = iscol;
3271:   PetscObjectReference((PetscObject)isrow);
3272:   PetscObjectReference((PetscObject)iscol);
3273:   b->icol = isicol;

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

3279:   B->factortype            = MAT_FACTOR_ILUDT;
3280:   B->info.factor_mallocs   = 0;
3281:   B->info.fill_ratio_given = ((PetscReal)nnz_max)/((PetscReal)ai[n]);
3282:   /* ------- end of symbolic factorization ---------*/

3284:   ISGetIndices(isrow,&r);
3285:   ISGetIndices(isicol,&ic);
3286:   ics  = ic;

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

3292:   /* im: used by PetscLLAddSortedLU(); jtmp: working array for column indices of active row */
3293:   PetscMalloc2(n,&im,n,&jtmp);
3294:   /* rtmp, vtmp: working arrays for sparse and contiguous row entries of active row */
3295:   PetscMalloc2(n,&rtmp,n,&vtmp);
3296:   PetscArrayzero(rtmp,n);

3298:   bi[0]        = 0;
3299:   bdiag[0]     = nnz_max-1; /* location of diag[0] in factor B */
3300:   bdiag_rev[n] = bdiag[0];
3301:   bi[2*n+1]    = bdiag[0]+1; /* endof bj and ba array */
3302:   for (i=0; i<n; i++) {
3303:     /* copy initial fill into linked list */
3304:     nzi = ai[r[i]+1] - ai[r[i]];
3305:     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);
3306:     nzi_al = adiag[r[i]] - ai[r[i]];
3307:     nzi_au = ai[r[i]+1] - adiag[r[i]] -1;
3308:     ajtmp  = aj + ai[r[i]];
3309:     PetscLLAddPerm(nzi,ajtmp,ic,n,nlnk,lnk,lnkbt);

3311:     /* load in initial (unfactored row) */
3312:     aatmp = a->a + ai[r[i]];
3313:     for (j=0; j<nzi; j++) {
3314:       rtmp[ics[*ajtmp++]] = *aatmp++;
3315:     }

3317:     /* add pivot rows into linked list */
3318:     row = lnk[n];
3319:     while (row < i) {
3320:       nzi_bl = bi[row+1] - bi[row] + 1;
3321:       bjtmp  = bj + bdiag[row+1]+1; /* points to 1st column next to the diagonal in U */
3322:       PetscLLAddSortedLU(bjtmp,row,nlnk,lnk,lnkbt,i,nzi_bl,im);
3323:       nzi   += nlnk;
3324:       row    = lnk[row];
3325:     }

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

3330:     /* numerical factorization */
3331:     bjtmp = jtmp;
3332:     row   = *bjtmp++; /* 1st pivot row */
3333:     while (row < i) {
3334:       pc         = rtmp + row;
3335:       pv         = ba + bdiag[row]; /* 1./(diag of the pivot row) */
3336:       multiplier = (*pc) * (*pv);
3337:       *pc        = multiplier;
3338:       if (PetscAbsScalar(*pc) > dt) { /* apply tolerance dropping rule */
3339:         pj = bj + bdiag[row+1] + 1;         /* point to 1st entry of U(row,:) */
3340:         pv = ba + bdiag[row+1] + 1;
3341:         /* if (multiplier < -1.0 or multiplier >1.0) printf("row/prow %d, %d, multiplier %g\n",i,row,multiplier); */
3342:         nz = bdiag[row] - bdiag[row+1] - 1;         /* num of entries in U(row,:), excluding diagonal */
3343:         for (j=0; j<nz; j++) rtmp[*pj++] -= multiplier * (*pv++);
3344:         PetscLogFlops(1+2*nz);
3345:       }
3346:       row = *bjtmp++;
3347:     }

3349:     /* copy sparse rtmp into contiguous vtmp; separate L and U part */
3350:     diag_tmp = rtmp[i];  /* save diagonal value - may not needed?? */
3351:     nzi_bl   = 0; j = 0;
3352:     while (jtmp[j] < i) { /* Note: jtmp is sorted */
3353:       vtmp[j] = rtmp[jtmp[j]]; rtmp[jtmp[j]]=0.0;
3354:       nzi_bl++; j++;
3355:     }
3356:     nzi_bu = nzi - nzi_bl -1;
3357:     while (j < nzi) {
3358:       vtmp[j] = rtmp[jtmp[j]]; rtmp[jtmp[j]]=0.0;
3359:       j++;
3360:     }

3362:     bjtmp = bj + bi[i];
3363:     batmp = ba + bi[i];
3364:     /* apply level dropping rule to L part */
3365:     ncut = nzi_al + dtcount;
3366:     if (ncut < nzi_bl) {
3367:       PetscSortSplit(ncut,nzi_bl,vtmp,jtmp);
3368:       PetscSortIntWithScalarArray(ncut,jtmp,vtmp);
3369:     } else {
3370:       ncut = nzi_bl;
3371:     }
3372:     for (j=0; j<ncut; j++) {
3373:       bjtmp[j] = jtmp[j];
3374:       batmp[j] = vtmp[j];
3375:       /* printf(" (%d,%g),",bjtmp[j],batmp[j]); */
3376:     }
3377:     bi[i+1] = bi[i] + ncut;
3378:     nzi     = ncut + 1;

3380:     /* apply level dropping rule to U part */
3381:     ncut = nzi_au + dtcount;
3382:     if (ncut < nzi_bu) {
3383:       PetscSortSplit(ncut,nzi_bu,vtmp+nzi_bl+1,jtmp+nzi_bl+1);
3384:       PetscSortIntWithScalarArray(ncut,jtmp+nzi_bl+1,vtmp+nzi_bl+1);
3385:     } else {
3386:       ncut = nzi_bu;
3387:     }
3388:     nzi += ncut;

3390:     /* mark bdiagonal */
3391:     bdiag[i+1]       = bdiag[i] - (ncut + 1);
3392:     bdiag_rev[n-i-1] = bdiag[i+1];
3393:     bi[2*n - i]      = bi[2*n - i +1] - (ncut + 1);
3394:     bjtmp            = bj + bdiag[i];
3395:     batmp            = ba + bdiag[i];
3396:     *bjtmp           = i;
3397:     *batmp           = diag_tmp; /* rtmp[i]; */
3398:     if (*batmp == 0.0) {
3399:       *batmp = dt+shift;
3400:       /* printf(" row %d add shift %g\n",i,shift); */
3401:     }
3402:     *batmp = 1.0/(*batmp); /* invert diagonal entries for simplier triangular solves */
3403:     /* printf(" (%d,%g),",*bjtmp,*batmp); */

3405:     bjtmp = bj + bdiag[i+1]+1;
3406:     batmp = ba + bdiag[i+1]+1;
3407:     for (k=0; k<ncut; k++) {
3408:       bjtmp[k] = jtmp[nzi_bl+1+k];
3409:       batmp[k] = vtmp[nzi_bl+1+k];
3410:       /* printf(" (%d,%g),",bjtmp[k],batmp[k]); */
3411:     }
3412:     /* printf("\n"); */

3414:     im[i] = nzi;   /* used by PetscLLAddSortedLU() */
3415:     /*
3416:     printf("row %d: bi %d, bdiag %d\n",i,bi[i],bdiag[i]);
3417:     printf(" ----------------------------\n");
3418:     */
3419:   } /* for (i=0; i<n; i++) */
3420:     /* printf("end of L %d, beginning of U %d\n",bi[n],bdiag[n]); */
3421:   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]);

3423:   ISRestoreIndices(isrow,&r);
3424:   ISRestoreIndices(isicol,&ic);

3426:   PetscLLDestroy(lnk,lnkbt);
3427:   PetscFree2(im,jtmp);
3428:   PetscFree2(rtmp,vtmp);
3429:   PetscFree(bdiag_rev);

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

3434:   ISIdentity(isrow,&row_identity);
3435:   ISIdentity(isicol,&icol_identity);
3436:   if (row_identity && icol_identity) {
3437:     B->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
3438:   } else {
3439:     B->ops->solve = MatSolve_SeqAIJ;
3440:   }

3442:   B->ops->solveadd          = 0;
3443:   B->ops->solvetranspose    = 0;
3444:   B->ops->solvetransposeadd = 0;
3445:   B->ops->matsolve          = 0;
3446:   B->assembled              = PETSC_TRUE;
3447:   B->preallocated           = PETSC_TRUE;
3448:   return(0);
3449: }

3451: /* a wraper of MatILUDTFactor_SeqAIJ() */
3452: /*
3453:     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
3454: */

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

3461:   MatILUDTFactor_SeqAIJ(A,row,col,info,&fact);
3462:   return(0);
3463: }

3465: /*
3466:    same as MatLUFactorNumeric_SeqAIJ(), except using contiguous array matrix factors
3467:    - intend to replace existing MatLUFactorNumeric_SeqAIJ()
3468: */
3469: /*
3470:     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
3471: */

3473: PetscErrorCode  MatILUDTFactorNumeric_SeqAIJ(Mat fact,Mat A,const MatFactorInfo *info)
3474: {
3475:   Mat            C     =fact;
3476:   Mat_SeqAIJ     *a    =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)C->data;
3477:   IS             isrow = b->row,isicol = b->icol;
3479:   const PetscInt *r,*ic,*ics;
3480:   PetscInt       i,j,k,n=A->rmap->n,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
3481:   PetscInt       *ajtmp,*bjtmp,nz,nzl,nzu,row,*bdiag = b->diag,*pj;
3482:   MatScalar      *rtmp,*pc,multiplier,*v,*pv,*aa=a->a;
3483:   PetscReal      dt=info->dt,shift=info->shiftamount;
3484:   PetscBool      row_identity, col_identity;

3487:   ISGetIndices(isrow,&r);
3488:   ISGetIndices(isicol,&ic);
3489:   PetscMalloc1(n+1,&rtmp);
3490:   ics  = ic;

3492:   for (i=0; i<n; i++) {
3493:     /* initialize rtmp array */
3494:     nzl   = bi[i+1] - bi[i];       /* num of nozeros in L(i,:) */
3495:     bjtmp = bj + bi[i];
3496:     for  (j=0; j<nzl; j++) rtmp[*bjtmp++] = 0.0;
3497:     rtmp[i] = 0.0;
3498:     nzu     = bdiag[i] - bdiag[i+1]; /* num of nozeros in U(i,:) */
3499:     bjtmp   = bj + bdiag[i+1] + 1;
3500:     for  (j=0; j<nzu; j++) rtmp[*bjtmp++] = 0.0;

3502:     /* load in initial unfactored row of A */
3503:     /* printf("row %d\n",i); */
3504:     nz    = ai[r[i]+1] - ai[r[i]];
3505:     ajtmp = aj + ai[r[i]];
3506:     v     = aa + ai[r[i]];
3507:     for (j=0; j<nz; j++) {
3508:       rtmp[ics[*ajtmp++]] = v[j];
3509:       /* printf(" (%d,%g),",ics[ajtmp[j]],rtmp[ics[ajtmp[j]]]); */
3510:     }
3511:     /* printf("\n"); */

3513:     /* numerical factorization */
3514:     bjtmp = bj + bi[i]; /* point to 1st entry of L(i,:) */
3515:     nzl   = bi[i+1] - bi[i]; /* num of entries in L(i,:) */
3516:     k     = 0;
3517:     while (k < nzl) {
3518:       row = *bjtmp++;
3519:       /* printf("  prow %d\n",row); */
3520:       pc         = rtmp + row;
3521:       pv         = b->a + bdiag[row]; /* 1./(diag of the pivot row) */
3522:       multiplier = (*pc) * (*pv);
3523:       *pc        = multiplier;
3524:       if (PetscAbsScalar(multiplier) > dt) {
3525:         pj = bj + bdiag[row+1] + 1;         /* point to 1st entry of U(row,:) */
3526:         pv = b->a + bdiag[row+1] + 1;
3527:         nz = bdiag[row] - bdiag[row+1] - 1;         /* num of entries in U(row,:), excluding diagonal */
3528:         for (j=0; j<nz; j++) rtmp[*pj++] -= multiplier * (*pv++);
3529:         PetscLogFlops(1+2*nz);
3530:       }
3531:       k++;
3532:     }

3534:     /* finished row so stick it into b->a */
3535:     /* L-part */
3536:     pv  = b->a + bi[i];
3537:     pj  = bj + bi[i];
3538:     nzl = bi[i+1] - bi[i];
3539:     for (j=0; j<nzl; j++) {
3540:       pv[j] = rtmp[pj[j]];
3541:       /* printf(" (%d,%g),",pj[j],pv[j]); */
3542:     }

3544:     /* diagonal: invert diagonal entries for simplier triangular solves */
3545:     if (rtmp[i] == 0.0) rtmp[i] = dt+shift;
3546:     b->a[bdiag[i]] = 1.0/rtmp[i];
3547:     /* printf(" (%d,%g),",i,b->a[bdiag[i]]); */

3549:     /* U-part */
3550:     pv  = b->a + bdiag[i+1] + 1;
3551:     pj  = bj + bdiag[i+1] + 1;
3552:     nzu = bdiag[i] - bdiag[i+1] - 1;
3553:     for (j=0; j<nzu; j++) {
3554:       pv[j] = rtmp[pj[j]];
3555:       /* printf(" (%d,%g),",pj[j],pv[j]); */
3556:     }
3557:     /* printf("\n"); */
3558:   }

3560:   PetscFree(rtmp);
3561:   ISRestoreIndices(isicol,&ic);
3562:   ISRestoreIndices(isrow,&r);

3564:   ISIdentity(isrow,&row_identity);
3565:   ISIdentity(isicol,&col_identity);
3566:   if (row_identity && col_identity) {
3567:     C->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
3568:   } else {
3569:     C->ops->solve = MatSolve_SeqAIJ;
3570:   }
3571:   C->ops->solveadd          = 0;
3572:   C->ops->solvetranspose    = 0;
3573:   C->ops->solvetransposeadd = 0;
3574:   C->ops->matsolve          = 0;
3575:   C->assembled              = PETSC_TRUE;
3576:   C->preallocated           = PETSC_TRUE;

3578:   PetscLogFlops(C->cmap->n);
3579:   return(0);
3580: }