Actual source code: aijfact.c

petsc-main 2021-04-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: static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type)
 93: {
 95:   *type = MATSOLVERPETSC;
 96:   return(0);
 97: }

 99: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_petsc(Mat A,MatFactorType ftype,Mat *B)
100: {
101:   PetscInt       n = A->rmap->n;

105: #if defined(PETSC_USE_COMPLEX)
106:   if (A->hermitian && !A->symmetric && (ftype == MAT_FACTOR_CHOLESKY||ftype == MAT_FACTOR_ICC)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian CHOLESKY or ICC Factor is not supported");
107: #endif
108:   MatCreate(PetscObjectComm((PetscObject)A),B);
109:   MatSetSizes(*B,n,n,n,n);
110:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU || ftype == MAT_FACTOR_ILUDT) {
111:     MatSetType(*B,MATSEQAIJ);

113:     (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqAIJ;
114:     (*B)->ops->lufactorsymbolic  = MatLUFactorSymbolic_SeqAIJ;

116:     MatSetBlockSizesFromMats(*B,A,A);
117:     PetscStrallocpy(MATORDERINGND,(char**)&(*B)->preferredordering[MAT_FACTOR_LU]);
118:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ILU]);
119:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ILUDT]);
120:   } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
121:     MatSetType(*B,MATSEQSBAIJ);
122:     MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,NULL);

124:     (*B)->ops->iccfactorsymbolic      = MatICCFactorSymbolic_SeqAIJ;
125:     (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqAIJ;
126:     PetscStrallocpy(MATORDERINGND,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]);
127:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC]);
128:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
129:   (*B)->factortype = ftype;

131:   PetscFree((*B)->solvertype);
132:   PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);
133:   (*B)->canuseordering = PETSC_TRUE;
134:   PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc);
135:   return(0);
136: }

138: PetscErrorCode MatLUFactorSymbolic_SeqAIJ_inplace(Mat B,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
139: {
140:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data,*b;
141:   IS                 isicol;
142:   PetscErrorCode     ierr;
143:   const PetscInt     *r,*ic;
144:   PetscInt           i,n=A->rmap->n,*ai=a->i,*aj=a->j;
145:   PetscInt           *bi,*bj,*ajtmp;
146:   PetscInt           *bdiag,row,nnz,nzi,reallocs=0,nzbd,*im;
147:   PetscReal          f;
148:   PetscInt           nlnk,*lnk,k,**bi_ptr;
149:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
150:   PetscBT            lnkbt;
151:   PetscBool          missing;

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

158:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
159:   ISGetIndices(isrow,&r);
160:   ISGetIndices(isicol,&ic);

162:   /* get new row pointers */
163:   PetscMalloc1(n+1,&bi);
164:   bi[0] = 0;

166:   /* bdiag is location of diagonal in factor */
167:   PetscMalloc1(n+1,&bdiag);
168:   bdiag[0] = 0;

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

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

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

182:   for (i=0; i<n; i++) {
183:     /* copy previous fill into linked list */
184:     nzi = 0;
185:     nnz = ai[r[i]+1] - ai[r[i]];
186:     ajtmp = aj + ai[r[i]];
187:     PetscLLAddPerm(nnz,ajtmp,ic,n,nlnk,lnk,lnkbt);
188:     nzi  += nlnk;

190:     /* add pivot rows into linked list */
191:     row = lnk[n];
192:     while (row < i) {
193:       nzbd  = bdiag[row] - bi[row] + 1;   /* num of entries in the row with column index <= row */
194:       ajtmp = bi_ptr[row] + nzbd;   /* points to the entry next to the diagonal */
195:       PetscLLAddSortedLU(ajtmp,row,nlnk,lnk,lnkbt,i,nzbd,im);
196:       nzi  += nlnk;
197:       row   = lnk[row];
198:     }
199:     bi[i+1] = bi[i] + nzi;
200:     im[i]   = nzi;

202:     /* mark bdiag */
203:     nzbd = 0;
204:     nnz  = nzi;
205:     k    = lnk[n];
206:     while (nnz-- && k < i) {
207:       nzbd++;
208:       k = lnk[k];
209:     }
210:     bdiag[i] = bi[i] + nzbd;

212:     /* if free space is not available, make more free space */
213:     if (current_space->local_remaining<nzi) {
214:       nnz  = PetscIntMultTruncate(n - i,nzi); /* estimated and max additional space needed */
215:       PetscFreeSpaceGet(nnz,&current_space);
216:       reallocs++;
217:     }

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

222:     bi_ptr[i]                       = current_space->array;
223:     current_space->array           += nzi;
224:     current_space->local_used      += nzi;
225:     current_space->local_remaining -= nzi;
226:   }
227: #if defined(PETSC_USE_INFO)
228:   if (ai[n] != 0) {
229:     PetscReal af = ((PetscReal)bi[n])/((PetscReal)ai[n]);
230:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);
231:     PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
232:     PetscInfo1(A,"PCFactorSetFill(pc,%g);\n",(double)af);
233:     PetscInfo(A,"for best performance.\n");
234:   } else {
235:     PetscInfo(A,"Empty matrix\n");
236:   }
237: #endif

239:   ISRestoreIndices(isrow,&r);
240:   ISRestoreIndices(isicol,&ic);

242:   /* destroy list of free space and other temporary array(s) */
243:   PetscMalloc1(bi[n]+1,&bj);
244:   PetscFreeSpaceContiguous(&free_space,bj);
245:   PetscLLDestroy(lnk,lnkbt);
246:   PetscFree2(bi_ptr,im);

248:   /* put together the new matrix */
249:   MatSeqAIJSetPreallocation_SeqAIJ(B,MAT_SKIP_ALLOCATION,NULL);
250:   PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);
251:   b    = (Mat_SeqAIJ*)(B)->data;

253:   b->free_a       = PETSC_TRUE;
254:   b->free_ij      = PETSC_TRUE;
255:   b->singlemalloc = PETSC_FALSE;

257:   PetscMalloc1(bi[n]+1,&b->a);
258:   b->j    = bj;
259:   b->i    = bi;
260:   b->diag = bdiag;
261:   b->ilen = NULL;
262:   b->imax = NULL;
263:   b->row  = isrow;
264:   b->col  = iscol;
265:   PetscObjectReference((PetscObject)isrow);
266:   PetscObjectReference((PetscObject)iscol);
267:   b->icol = isicol;
268:   PetscMalloc1(n+1,&b->solve_work);

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

274:   (B)->factortype            = MAT_FACTOR_LU;
275:   (B)->info.factor_mallocs   = reallocs;
276:   (B)->info.fill_ratio_given = f;

278:   if (ai[n]) {
279:     (B)->info.fill_ratio_needed = ((PetscReal)bi[n])/((PetscReal)ai[n]);
280:   } else {
281:     (B)->info.fill_ratio_needed = 0.0;
282:   }
283:   (B)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_inplace;
284:   if (a->inode.size) {
285:     (B)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode_inplace;
286:   }
287:   return(0);
288: }

290: PetscErrorCode MatLUFactorSymbolic_SeqAIJ(Mat B,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
291: {
292:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data,*b;
293:   IS                 isicol;
294:   PetscErrorCode     ierr;
295:   const PetscInt     *r,*ic,*ai=a->i,*aj=a->j,*ajtmp;
296:   PetscInt           i,n=A->rmap->n;
297:   PetscInt           *bi,*bj;
298:   PetscInt           *bdiag,row,nnz,nzi,reallocs=0,nzbd,*im;
299:   PetscReal          f;
300:   PetscInt           nlnk,*lnk,k,**bi_ptr;
301:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
302:   PetscBT            lnkbt;
303:   PetscBool          missing;

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

310:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
311:   ISGetIndices(isrow,&r);
312:   ISGetIndices(isicol,&ic);

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

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

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

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

331:   for (i=0; i<n; i++) {
332:     /* copy previous fill into linked list */
333:     nzi = 0;
334:     nnz = ai[r[i]+1] - ai[r[i]];
335:     ajtmp = aj + ai[r[i]];
336:     PetscLLAddPerm(nnz,ajtmp,ic,n,nlnk,lnk,lnkbt);
337:     nzi  += nlnk;

339:     /* add pivot rows into linked list */
340:     row = lnk[n];
341:     while (row < i) {
342:       nzbd  = bdiag[row] + 1; /* num of entries in the row with column index <= row */
343:       ajtmp = bi_ptr[row] + nzbd; /* points to the entry next to the diagonal */
344:       PetscLLAddSortedLU(ajtmp,row,nlnk,lnk,lnkbt,i,nzbd,im);
345:       nzi  += nlnk;
346:       row   = lnk[row];
347:     }
348:     bi[i+1] = bi[i] + nzi;
349:     im[i]   = nzi;

351:     /* mark bdiag */
352:     nzbd = 0;
353:     nnz  = nzi;
354:     k    = lnk[n];
355:     while (nnz-- && k < i) {
356:       nzbd++;
357:       k = lnk[k];
358:     }
359:     bdiag[i] = nzbd; /* note: bdiag[i] = nnzL as input for PetscFreeSpaceContiguous_LU() */

361:     /* if free space is not available, make more free space */
362:     if (current_space->local_remaining<nzi) {
363:       /* estimated additional space needed */
364:       nnz  = PetscIntMultTruncate(2,PetscIntMultTruncate(n-1,nzi));
365:       PetscFreeSpaceGet(nnz,&current_space);
366:       reallocs++;
367:     }

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

372:     bi_ptr[i]                       = current_space->array;
373:     current_space->array           += nzi;
374:     current_space->local_used      += nzi;
375:     current_space->local_remaining -= nzi;
376:   }

378:   ISRestoreIndices(isrow,&r);
379:   ISRestoreIndices(isicol,&ic);

381:   /*   copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
382:   PetscMalloc1(bi[n]+1,&bj);
383:   PetscFreeSpaceContiguous_LU(&free_space,bj,n,bi,bdiag);
384:   PetscLLDestroy(lnk,lnkbt);
385:   PetscFree2(bi_ptr,im);

387:   /* put together the new matrix */
388:   MatSeqAIJSetPreallocation_SeqAIJ(B,MAT_SKIP_ALLOCATION,NULL);
389:   PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);
390:   b    = (Mat_SeqAIJ*)(B)->data;

392:   b->free_a       = PETSC_TRUE;
393:   b->free_ij      = PETSC_TRUE;
394:   b->singlemalloc = PETSC_FALSE;

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

398:   b->j    = bj;
399:   b->i    = bi;
400:   b->diag = bdiag;
401:   b->ilen = NULL;
402:   b->imax = NULL;
403:   b->row  = isrow;
404:   b->col  = iscol;
405:   PetscObjectReference((PetscObject)isrow);
406:   PetscObjectReference((PetscObject)iscol);
407:   b->icol = isicol;
408:   PetscMalloc1(n+1,&b->solve_work);

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

414:   B->factortype            = MAT_FACTOR_LU;
415:   B->info.factor_mallocs   = reallocs;
416:   B->info.fill_ratio_given = f;

418:   if (ai[n]) {
419:     B->info.fill_ratio_needed = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]);
420:   } else {
421:     B->info.fill_ratio_needed = 0.0;
422:   }
423: #if defined(PETSC_USE_INFO)
424:   if (ai[n] != 0) {
425:     PetscReal af = B->info.fill_ratio_needed;
426:     PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);
427:     PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
428:     PetscInfo1(A,"PCFactorSetFill(pc,%g);\n",(double)af);
429:     PetscInfo(A,"for best performance.\n");
430:   } else {
431:     PetscInfo(A,"Empty matrix\n");
432:   }
433: #endif
434:   B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ;
435:   if (a->inode.size) {
436:     B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
437:   }
438:   MatSeqAIJCheckInode_FactorLU(B);
439:   return(0);
440: }

442: /*
443:     Trouble in factorization, should we dump the original matrix?
444: */
445: PetscErrorCode MatFactorDumpMatrix(Mat A)
446: {
448:   PetscBool      flg = PETSC_FALSE;

451:   PetscOptionsGetBool(((PetscObject)A)->options,NULL,"-mat_factor_dump_on_error",&flg,NULL);
452:   if (flg) {
453:     PetscViewer viewer;
454:     char        filename[PETSC_MAX_PATH_LEN];

456:     PetscSNPrintf(filename,PETSC_MAX_PATH_LEN,"matrix_factor_error.%d",PetscGlobalRank);
457:     PetscViewerBinaryOpen(PetscObjectComm((PetscObject)A),filename,FILE_MODE_WRITE,&viewer);
458:     MatView(A,viewer);
459:     PetscViewerDestroy(&viewer);
460:   }
461:   return(0);
462: }

464: PetscErrorCode MatLUFactorNumeric_SeqAIJ(Mat B,Mat A,const MatFactorInfo *info)
465: {
466:   Mat             C     =B;
467:   Mat_SeqAIJ      *a    =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)C->data;
468:   IS              isrow = b->row,isicol = b->icol;
469:   PetscErrorCode  ierr;
470:   const PetscInt  *r,*ic,*ics;
471:   const PetscInt  n=A->rmap->n,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bdiag=b->diag;
472:   PetscInt        i,j,k,nz,nzL,row,*pj;
473:   const PetscInt  *ajtmp,*bjtmp;
474:   MatScalar       *rtmp,*pc,multiplier,*pv;
475:   const MatScalar *aa=a->a,*v;
476:   PetscBool       row_identity,col_identity;
477:   FactorShiftCtx  sctx;
478:   const PetscInt  *ddiag;
479:   PetscReal       rs;
480:   MatScalar       d;

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

486:   if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
487:     ddiag          = a->diag;
488:     sctx.shift_top = info->zeropivot;
489:     for (i=0; i<n; i++) {
490:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
491:       d  = (aa)[ddiag[i]];
492:       rs = -PetscAbsScalar(d) - PetscRealPart(d);
493:       v  = aa+ai[i];
494:       nz = ai[i+1] - ai[i];
495:       for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
496:       if (rs>sctx.shift_top) sctx.shift_top = rs;
497:     }
498:     sctx.shift_top *= 1.1;
499:     sctx.nshift_max = 5;
500:     sctx.shift_lo   = 0.;
501:     sctx.shift_hi   = 1.;
502:   }

504:   ISGetIndices(isrow,&r);
505:   ISGetIndices(isicol,&ic);
506:   PetscMalloc1(n+1,&rtmp);
507:   ics  = ic;

509:   do {
510:     sctx.newshift = PETSC_FALSE;
511:     for (i=0; i<n; i++) {
512:       /* zero rtmp */
513:       /* L part */
514:       nz    = bi[i+1] - bi[i];
515:       bjtmp = bj + bi[i];
516:       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

518:       /* U part */
519:       nz    = bdiag[i]-bdiag[i+1];
520:       bjtmp = bj + bdiag[i+1]+1;
521:       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

523:       /* load in initial (unfactored row) */
524:       nz    = ai[r[i]+1] - ai[r[i]];
525:       ajtmp = aj + ai[r[i]];
526:       v     = aa + ai[r[i]];
527:       for (j=0; j<nz; j++) {
528:         rtmp[ics[ajtmp[j]]] = v[j];
529:       }
530:       /* ZeropivotApply() */
531:       rtmp[i] += sctx.shift_amount;  /* shift the diagonal of the matrix */

533:       /* elimination */
534:       bjtmp = bj + bi[i];
535:       row   = *bjtmp++;
536:       nzL   = bi[i+1] - bi[i];
537:       for (k=0; k < nzL; k++) {
538:         pc = rtmp + row;
539:         if (*pc != 0.0) {
540:           pv         = b->a + bdiag[row];
541:           multiplier = *pc * (*pv);
542:           *pc        = multiplier;

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

548:           for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
549:           PetscLogFlops(1+2.0*nz);
550:         }
551:         row = *bjtmp++;
552:       }

554:       /* finished row so stick it into b->a */
555:       rs = 0.0;
556:       /* L part */
557:       pv = b->a + bi[i];
558:       pj = b->j + bi[i];
559:       nz = bi[i+1] - bi[i];
560:       for (j=0; j<nz; j++) {
561:         pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
562:       }

564:       /* U part */
565:       pv = b->a + bdiag[i+1]+1;
566:       pj = b->j + bdiag[i+1]+1;
567:       nz = bdiag[i] - bdiag[i+1]-1;
568:       for (j=0; j<nz; j++) {
569:         pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
570:       }

572:       sctx.rs = rs;
573:       sctx.pv = rtmp[i];
574:       MatPivotCheck(B,A,info,&sctx,i);
575:       if (sctx.newshift) break; /* break for-loop */
576:       rtmp[i] = sctx.pv; /* sctx.pv might be updated in the case of MAT_SHIFT_INBLOCKS */

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

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

584:     /* MatPivotRefine() */
585:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
586:       /*
587:        * if no shift in this attempt & shifting & started shifting & can refine,
588:        * then try lower shift
589:        */
590:       sctx.shift_hi       = sctx.shift_fraction;
591:       sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
592:       sctx.shift_amount   = sctx.shift_fraction * sctx.shift_top;
593:       sctx.newshift       = PETSC_TRUE;
594:       sctx.nshift++;
595:     }
596:   } while (sctx.newshift);

598:   PetscFree(rtmp);
599:   ISRestoreIndices(isicol,&ic);
600:   ISRestoreIndices(isrow,&r);

602:   ISIdentity(isrow,&row_identity);
603:   ISIdentity(isicol,&col_identity);
604:   if (b->inode.size) {
605:     C->ops->solve = MatSolve_SeqAIJ_Inode;
606:   } else if (row_identity && col_identity) {
607:     C->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
608:   } else {
609:     C->ops->solve = MatSolve_SeqAIJ;
610:   }
611:   C->ops->solveadd          = MatSolveAdd_SeqAIJ;
612:   C->ops->solvetranspose    = MatSolveTranspose_SeqAIJ;
613:   C->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ;
614:   C->ops->matsolve          = MatMatSolve_SeqAIJ;
615:   C->assembled              = PETSC_TRUE;
616:   C->preallocated           = PETSC_TRUE;

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

620:   /* MatShiftView(A,info,&sctx) */
621:   if (sctx.nshift) {
622:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
623:       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);
624:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
625:       PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
626:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
627:       PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %g\n",sctx.nshift,(double)info->shiftamount);
628:     }
629:   }
630:   return(0);
631: }

633: PetscErrorCode MatLUFactorNumeric_SeqAIJ_inplace(Mat B,Mat A,const MatFactorInfo *info)
634: {
635:   Mat             C     =B;
636:   Mat_SeqAIJ      *a    =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)C->data;
637:   IS              isrow = b->row,isicol = b->icol;
638:   PetscErrorCode  ierr;
639:   const PetscInt  *r,*ic,*ics;
640:   PetscInt        nz,row,i,j,n=A->rmap->n,diag;
641:   const PetscInt  *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
642:   const PetscInt  *ajtmp,*bjtmp,*diag_offset = b->diag,*pj;
643:   MatScalar       *pv,*rtmp,*pc,multiplier,d;
644:   const MatScalar *v,*aa=a->a;
645:   PetscReal       rs=0.0;
646:   FactorShiftCtx  sctx;
647:   const PetscInt  *ddiag;
648:   PetscBool       row_identity, col_identity;

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

654:   if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
655:     ddiag          = a->diag;
656:     sctx.shift_top = info->zeropivot;
657:     for (i=0; i<n; i++) {
658:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
659:       d  = (aa)[ddiag[i]];
660:       rs = -PetscAbsScalar(d) - PetscRealPart(d);
661:       v  = aa+ai[i];
662:       nz = ai[i+1] - ai[i];
663:       for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
664:       if (rs>sctx.shift_top) sctx.shift_top = rs;
665:     }
666:     sctx.shift_top *= 1.1;
667:     sctx.nshift_max = 5;
668:     sctx.shift_lo   = 0.;
669:     sctx.shift_hi   = 1.;
670:   }

672:   ISGetIndices(isrow,&r);
673:   ISGetIndices(isicol,&ic);
674:   PetscMalloc1(n+1,&rtmp);
675:   ics  = ic;

677:   do {
678:     sctx.newshift = PETSC_FALSE;
679:     for (i=0; i<n; i++) {
680:       nz    = bi[i+1] - bi[i];
681:       bjtmp = bj + bi[i];
682:       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

684:       /* load in initial (unfactored row) */
685:       nz    = ai[r[i]+1] - ai[r[i]];
686:       ajtmp = aj + ai[r[i]];
687:       v     = aa + ai[r[i]];
688:       for (j=0; j<nz; j++) {
689:         rtmp[ics[ajtmp[j]]] = v[j];
690:       }
691:       rtmp[ics[r[i]]] += sctx.shift_amount; /* shift the diagonal of the matrix */

693:       row = *bjtmp++;
694:       while  (row < i) {
695:         pc = rtmp + row;
696:         if (*pc != 0.0) {
697:           pv         = b->a + diag_offset[row];
698:           pj         = b->j + diag_offset[row] + 1;
699:           multiplier = *pc / *pv++;
700:           *pc        = multiplier;
701:           nz         = bi[row+1] - diag_offset[row] - 1;
702:           for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
703:           PetscLogFlops(1+2.0*nz);
704:         }
705:         row = *bjtmp++;
706:       }
707:       /* finished row so stick it into b->a */
708:       pv   = b->a + bi[i];
709:       pj   = b->j + bi[i];
710:       nz   = bi[i+1] - bi[i];
711:       diag = diag_offset[i] - bi[i];
712:       rs   = 0.0;
713:       for (j=0; j<nz; j++) {
714:         pv[j] = rtmp[pj[j]];
715:         rs   += PetscAbsScalar(pv[j]);
716:       }
717:       rs -= PetscAbsScalar(pv[diag]);

719:       sctx.rs = rs;
720:       sctx.pv = pv[diag];
721:       MatPivotCheck(B,A,info,&sctx,i);
722:       if (sctx.newshift) break;
723:       pv[diag] = sctx.pv;
724:     }

726:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
727:       /*
728:        * if no shift in this attempt & shifting & started shifting & can refine,
729:        * then try lower shift
730:        */
731:       sctx.shift_hi       = sctx.shift_fraction;
732:       sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
733:       sctx.shift_amount   = sctx.shift_fraction * sctx.shift_top;
734:       sctx.newshift       = PETSC_TRUE;
735:       sctx.nshift++;
736:     }
737:   } while (sctx.newshift);

739:   /* invert diagonal entries for simplier triangular solves */
740:   for (i=0; i<n; i++) {
741:     b->a[diag_offset[i]] = 1.0/b->a[diag_offset[i]];
742:   }
743:   PetscFree(rtmp);
744:   ISRestoreIndices(isicol,&ic);
745:   ISRestoreIndices(isrow,&r);

747:   ISIdentity(isrow,&row_identity);
748:   ISIdentity(isicol,&col_identity);
749:   if (row_identity && col_identity) {
750:     C->ops->solve = MatSolve_SeqAIJ_NaturalOrdering_inplace;
751:   } else {
752:     C->ops->solve = MatSolve_SeqAIJ_inplace;
753:   }
754:   C->ops->solveadd          = MatSolveAdd_SeqAIJ_inplace;
755:   C->ops->solvetranspose    = MatSolveTranspose_SeqAIJ_inplace;
756:   C->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ_inplace;
757:   C->ops->matsolve          = MatMatSolve_SeqAIJ_inplace;

759:   C->assembled    = PETSC_TRUE;
760:   C->preallocated = PETSC_TRUE;

762:   PetscLogFlops(C->cmap->n);
763:   if (sctx.nshift) {
764:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
765:       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);
766:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
767:       PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
768:     }
769:   }
770:   (C)->ops->solve          = MatSolve_SeqAIJ_inplace;
771:   (C)->ops->solvetranspose = MatSolveTranspose_SeqAIJ_inplace;

773:   MatSeqAIJCheckInode(C);
774:   return(0);
775: }

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

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

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

808:   if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
809:     const PetscInt *ddiag = a->diag;
810:     sctx.shift_top = info->zeropivot;
811:     for (i=0; i<n; i++) {
812:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
813:       d    = (aa)[ddiag[i]];
814:       rs   = -PetscAbsScalar(d) - PetscRealPart(d);
815:       vtmp = aa+ai[i];
816:       nz   = ai[i+1] - ai[i];
817:       for (j=0; j<nz; j++) rs += PetscAbsScalar(vtmp[j]);
818:       if (rs>sctx.shift_top) sctx.shift_top = rs;
819:     }
820:     sctx.shift_top *= 1.1;
821:     sctx.nshift_max = 5;
822:     sctx.shift_lo   = 0.;
823:     sctx.shift_hi   = 1.;
824:   }

826:   ISGetIndices(isrow,&r);
827:   ISGetIndices(isicol,&ic);
828:   PetscMalloc1(n+1,&rtmp);
829:   PetscArrayzero(rtmp,n+1);
830:   ics  = ic;

832: #if defined(MV)
833:   sctx.shift_top      = 0.;
834:   sctx.nshift_max     = 0;
835:   sctx.shift_lo       = 0.;
836:   sctx.shift_hi       = 0.;
837:   sctx.shift_fraction = 0.;

839:   if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
840:     sctx.shift_top = 0.;
841:     for (i=0; i<n; i++) {
842:       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
843:       d  = (a->a)[diag[i]];
844:       rs = -PetscAbsScalar(d) - PetscRealPart(d);
845:       v  = a->a+ai[i];
846:       nz = ai[i+1] - ai[i];
847:       for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
848:       if (rs>sctx.shift_top) sctx.shift_top = rs;
849:     }
850:     if (sctx.shift_top < info->zeropivot) sctx.shift_top = info->zeropivot;
851:     sctx.shift_top *= 1.1;
852:     sctx.nshift_max = 5;
853:     sctx.shift_lo   = 0.;
854:     sctx.shift_hi   = 1.;
855:   }

857:   sctx.shift_amount = 0.;
858:   sctx.nshift       = 0;
859: #endif

861:   do {
862:     sctx.newshift = PETSC_FALSE;
863:     for (i=0; i<n; i++) {
864:       /* load in initial unfactored row */
865:       nz    = ai[r[i]+1] - ai[r[i]];
866:       ajtmp = aj + ai[r[i]];
867:       v     = a->a + ai[r[i]];
868:       /* sort permuted ajtmp and values v accordingly */
869:       for (j=0; j<nz; j++) ajtmp[j] = ics[ajtmp[j]];
870:       PetscSortIntWithScalarArray(nz,ajtmp,v);

872:       diag[r[i]] = ai[r[i]];
873:       for (j=0; j<nz; j++) {
874:         rtmp[ajtmp[j]] = v[j];
875:         if (ajtmp[j] < i) diag[r[i]]++; /* update a->diag */
876:       }
877:       rtmp[r[i]] += sctx.shift_amount; /* shift the diagonal of the matrix */

879:       row = *ajtmp++;
880:       while  (row < i) {
881:         pc = rtmp + row;
882:         if (*pc != 0.0) {
883:           pv = a->a + diag[r[row]];
884:           pj = aj + diag[r[row]] + 1;

886:           multiplier = *pc / *pv++;
887:           *pc        = multiplier;
888:           nz         = ai[r[row]+1] - diag[r[row]] - 1;
889:           for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
890:           PetscLogFlops(1+2.0*nz);
891:         }
892:         row = *ajtmp++;
893:       }
894:       /* finished row so overwrite it onto a->a */
895:       pv     = a->a + ai[r[i]];
896:       pj     = aj + ai[r[i]];
897:       nz     = ai[r[i]+1] - ai[r[i]];
898:       nbdiag = diag[r[i]] - ai[r[i]]; /* num of entries before the diagonal */

900:       rs = 0.0;
901:       for (j=0; j<nz; j++) {
902:         pv[j] = rtmp[pj[j]];
903:         if (j != nbdiag) rs += PetscAbsScalar(pv[j]);
904:       }

906:       sctx.rs = rs;
907:       sctx.pv = pv[nbdiag];
908:       MatPivotCheck(B,A,info,&sctx,i);
909:       if (sctx.newshift) break;
910:       pv[nbdiag] = sctx.pv;
911:     }

913:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
914:       /*
915:        * if no shift in this attempt & shifting & started shifting & can refine,
916:        * then try lower shift
917:        */
918:       sctx.shift_hi       = sctx.shift_fraction;
919:       sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
920:       sctx.shift_amount   = sctx.shift_fraction * sctx.shift_top;
921:       sctx.newshift       = PETSC_TRUE;
922:       sctx.nshift++;
923:     }
924:   } while (sctx.newshift);

926:   /* invert diagonal entries for simplier triangular solves */
927:   for (i=0; i<n; i++) {
928:     a->a[diag[r[i]]] = 1.0/a->a[diag[r[i]]];
929:   }

931:   PetscFree(rtmp);
932:   ISRestoreIndices(isicol,&ic);
933:   ISRestoreIndices(isrow,&r);

935:   A->ops->solve             = MatSolve_SeqAIJ_InplaceWithPerm;
936:   A->ops->solveadd          = MatSolveAdd_SeqAIJ_inplace;
937:   A->ops->solvetranspose    = MatSolveTranspose_SeqAIJ_inplace;
938:   A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ_inplace;

940:   A->assembled    = PETSC_TRUE;
941:   A->preallocated = PETSC_TRUE;

943:   PetscLogFlops(A->cmap->n);
944:   if (sctx.nshift) {
945:     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
946:       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);
947:     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
948:       PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
949:     }
950:   }
951:   return(0);
952: }

954: /* ----------------------------------------------------------- */
955: PetscErrorCode MatLUFactor_SeqAIJ(Mat A,IS row,IS col,const MatFactorInfo *info)
956: {
958:   Mat            C;

961:   MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&C);
962:   MatLUFactorSymbolic(C,A,row,col,info);
963:   MatLUFactorNumeric(C,A,info);

965:   A->ops->solve          = C->ops->solve;
966:   A->ops->solvetranspose = C->ops->solvetranspose;

968:   MatHeaderMerge(A,&C);
969:   PetscLogObjectParent((PetscObject)A,(PetscObject)((Mat_SeqAIJ*)(A->data))->icol);
970:   return(0);
971: }
972: /* ----------------------------------------------------------- */


975: PetscErrorCode MatSolve_SeqAIJ_inplace(Mat A,Vec bb,Vec xx)
976: {
977:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
978:   IS                iscol = a->col,isrow = a->row;
979:   PetscErrorCode    ierr;
980:   PetscInt          i, n = A->rmap->n,*vi,*ai = a->i,*aj = a->j;
981:   PetscInt          nz;
982:   const PetscInt    *rout,*cout,*r,*c;
983:   PetscScalar       *x,*tmp,*tmps,sum;
984:   const PetscScalar *b;
985:   const MatScalar   *aa = a->a,*v;

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

990:   VecGetArrayRead(bb,&b);
991:   VecGetArrayWrite(xx,&x);
992:   tmp  = a->solve_work;

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

997:   /* forward solve the lower triangular */
998:   tmp[0] = b[*r++];
999:   tmps   = tmp;
1000:   for (i=1; i<n; i++) {
1001:     v   = aa + ai[i];
1002:     vi  = aj + ai[i];
1003:     nz  = a->diag[i] - ai[i];
1004:     sum = b[*r++];
1005:     PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
1006:     tmp[i] = sum;
1007:   }

1009:   /* backward solve the upper triangular */
1010:   for (i=n-1; i>=0; i--) {
1011:     v   = aa + a->diag[i] + 1;
1012:     vi  = aj + a->diag[i] + 1;
1013:     nz  = ai[i+1] - a->diag[i] - 1;
1014:     sum = tmp[i];
1015:     PetscSparseDenseMinusDot(sum,tmps,v,vi,nz);
1016:     x[*c--] = tmp[i] = sum*aa[a->diag[i]];
1017:   }

1019:   ISRestoreIndices(isrow,&rout);
1020:   ISRestoreIndices(iscol,&cout);
1021:   VecRestoreArrayRead(bb,&b);
1022:   VecRestoreArrayWrite(xx,&x);
1023:   PetscLogFlops(2.0*a->nz - A->cmap->n);
1024:   return(0);
1025: }

1027: PetscErrorCode MatMatSolve_SeqAIJ_inplace(Mat A,Mat B,Mat X)
1028: {
1029:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1030:   IS                iscol = a->col,isrow = a->row;
1031:   PetscErrorCode    ierr;
1032:   PetscInt          i, n = A->rmap->n,*vi,*ai = a->i,*aj = a->j;
1033:   PetscInt          nz,neq,ldb,ldx;
1034:   const PetscInt    *rout,*cout,*r,*c;
1035:   PetscScalar       *x,*tmp = a->solve_work,*tmps,sum;
1036:   const PetscScalar *b,*aa = a->a,*v;
1037:   PetscBool         isdense;

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

1085: PetscErrorCode MatMatSolve_SeqAIJ(Mat A,Mat B,Mat X)
1086: {
1087:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ*)A->data;
1088:   IS                iscol = a->col,isrow = a->row;
1089:   PetscErrorCode    ierr;
1090:   PetscInt          i, n = A->rmap->n,*vi,*ai = a->i,*aj = a->j,*adiag = a->diag;
1091:   PetscInt          nz,neq,ldb,ldx;
1092:   const PetscInt    *rout,*cout,*r,*c;
1093:   PetscScalar       *x,*tmp = a->solve_work,sum;
1094:   const PetscScalar *b,*aa = a->a,*v;
1095:   PetscBool         isdense;

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

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

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

1158:   VecGetArrayRead(bb,&b);
1159:   VecGetArrayWrite(xx,&x);
1160:   tmp  = a->solve_work;

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

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

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

1189:   ISRestoreIndices(isrow,&rout);
1190:   ISRestoreIndices(iscol,&cout);
1191:   VecRestoreArrayRead(bb,&b);
1192:   VecRestoreArrayWrite(xx,&x);
1193:   PetscLogFlops(2.0*a->nz - A->cmap->n);
1194:   return(0);
1195: }

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

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

1218:   VecGetArrayRead(bb,&b);
1219:   VecGetArrayWrite(xx,&x);

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

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

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

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

1268:   VecGetArrayRead(bb,&b);
1269:   VecGetArray(xx,&x);
1270:   tmp  = a->solve_work;

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

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

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

1297:   ISRestoreIndices(isrow,&rout);
1298:   ISRestoreIndices(iscol,&cout);
1299:   VecRestoreArrayRead(bb,&b);
1300:   VecRestoreArray(xx,&x);
1301:   PetscLogFlops(2.0*a->nz);
1302:   return(0);
1303: }

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

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

1320:   VecGetArrayRead(bb,&b);
1321:   VecGetArray(xx,&x);
1322:   tmp  = a->solve_work;

1324:   ISGetIndices(isrow,&rout); r = rout;
1325:   ISGetIndices(iscol,&cout); c = cout;

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

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

1352:   ISRestoreIndices(isrow,&rout);
1353:   ISRestoreIndices(iscol,&cout);
1354:   VecRestoreArrayRead(bb,&b);
1355:   VecRestoreArray(xx,&x);
1356:   PetscLogFlops(2.0*a->nz);
1357:   return(0);
1358: }

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

1373:   VecGetArrayRead(bb,&b);
1374:   VecGetArrayWrite(xx,&x);
1375:   tmp  = a->solve_work;

1377:   ISGetIndices(isrow,&rout); r = rout;
1378:   ISGetIndices(iscol,&cout); c = cout;

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

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

1394:   /* backward solve the L^T */
1395:   for (i=n-1; i>=0; i--) {
1396:     v  = aa + diag[i] - 1;
1397:     vi = aj + diag[i] - 1;
1398:     nz = diag[i] - ai[i];
1399:     s1 = tmp[i];
1400:     for (j=0; j>-nz; j--) tmp[vi[j]] -= s1*v[j];
1401:   }

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

1406:   ISRestoreIndices(isrow,&rout);
1407:   ISRestoreIndices(iscol,&cout);
1408:   VecRestoreArrayRead(bb,&b);
1409:   VecRestoreArrayWrite(xx,&x);

1411:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1412:   return(0);
1413: }

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

1428:   VecGetArrayRead(bb,&b);
1429:   VecGetArrayWrite(xx,&x);
1430:   tmp  = a->solve_work;

1432:   ISGetIndices(isrow,&rout); r = rout;
1433:   ISGetIndices(iscol,&cout); c = cout;

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

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

1449:   /* backward solve the L^T */
1450:   for (i=n-1; i>=0; i--) {
1451:     v  = aa + ai[i];
1452:     vi = aj + ai[i];
1453:     nz = ai[i+1] - ai[i];
1454:     s1 = tmp[i];
1455:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1456:   }

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

1461:   ISRestoreIndices(isrow,&rout);
1462:   ISRestoreIndices(iscol,&cout);
1463:   VecRestoreArrayRead(bb,&b);
1464:   VecRestoreArrayWrite(xx,&x);

1466:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1467:   return(0);
1468: }

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

1483:   if (zz != xx) {VecCopy(zz,xx);}
1484:   VecGetArrayRead(bb,&b);
1485:   VecGetArray(xx,&x);
1486:   tmp  = a->solve_work;

1488:   ISGetIndices(isrow,&rout); r = rout;
1489:   ISGetIndices(iscol,&cout); c = cout;

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

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

1505:   /* backward solve the L^T */
1506:   for (i=n-1; i>=0; i--) {
1507:     v  = aa + diag[i] - 1;
1508:     vi = aj + diag[i] - 1;
1509:     nz = diag[i] - ai[i];
1510:     s1 = tmp[i];
1511:     for (j=0; j>-nz; j--) tmp[vi[j]] -= s1*v[j];
1512:   }

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

1517:   ISRestoreIndices(isrow,&rout);
1518:   ISRestoreIndices(iscol,&cout);
1519:   VecRestoreArrayRead(bb,&b);
1520:   VecRestoreArray(xx,&x);

1522:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1523:   return(0);
1524: }

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

1539:   if (zz != xx) {VecCopy(zz,xx);}
1540:   VecGetArrayRead(bb,&b);
1541:   VecGetArray(xx,&x);
1542:   tmp  = a->solve_work;

1544:   ISGetIndices(isrow,&rout); r = rout;
1545:   ISGetIndices(iscol,&cout); c = cout;

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

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


1562:   /* backward solve the L^T */
1563:   for (i=n-1; i>=0; i--) {
1564:     v  = aa + ai[i];
1565:     vi = aj + ai[i];
1566:     nz = ai[i+1] - ai[i];
1567:     s1 = tmp[i];
1568:     for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j];
1569:   }

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

1574:   ISRestoreIndices(isrow,&rout);
1575:   ISRestoreIndices(iscol,&cout);
1576:   VecRestoreArrayRead(bb,&b);
1577:   VecRestoreArray(xx,&x);

1579:   PetscLogFlops(2.0*a->nz-A->cmap->n);
1580:   return(0);
1581: }

1583: /* ----------------------------------------------------------------*/

1585: /*
1586:    ilu() under revised new data structure.
1587:    Factored arrays bj and ba are stored as
1588:      L(0,:), L(1,:), ...,L(n-1,:),  U(n-1,:),...,U(i,:),U(i-1,:),...,U(0,:)

1590:    bi=fact->i is an array of size n+1, in which
1591:    bi+
1592:      bi[i]:  points to 1st entry of L(i,:),i=0,...,n-1
1593:      bi[n]:  points to L(n-1,n-1)+1

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

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

1611:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
1612:   MatDuplicateNoCreate_SeqAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_FALSE);
1613:   b    = (Mat_SeqAIJ*)(fact)->data;

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

1619:   b->singlemalloc = PETSC_TRUE;
1620:   if (!b->diag) {
1621:     PetscMalloc1(n+1,&b->diag);
1622:     PetscLogObjectMemory((PetscObject)fact,(n+1)*sizeof(PetscInt));
1623:   }
1624:   bdiag = b->diag;

1626:   if (n > 0) {
1627:     PetscArrayzero(b->a,ai[n]);
1628:   }

1630:   /* set bi and bj with new data structure */
1631:   bi = b->i;
1632:   bj = b->j;

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

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

1661:   fact->factortype             = MAT_FACTOR_ILU;
1662:   fact->info.factor_mallocs    = 0;
1663:   fact->info.fill_ratio_given  = info->fill;
1664:   fact->info.fill_ratio_needed = 1.0;
1665:   fact->ops->lufactornumeric   = MatLUFactorNumeric_SeqAIJ;
1666:   MatSeqAIJCheckInode_FactorLU(fact);

1668:   b       = (Mat_SeqAIJ*)(fact)->data;
1669:   b->row  = isrow;
1670:   b->col  = iscol;
1671:   b->icol = isicol;
1672:   PetscMalloc1(fact->rmap->n+1,&b->solve_work);
1673:   PetscObjectReference((PetscObject)isrow);
1674:   PetscObjectReference((PetscObject)iscol);
1675:   return(0);
1676: }

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

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

1701:   levels = (PetscInt)info->levels;
1702:   ISIdentity(isrow,&row_identity);
1703:   ISIdentity(iscol,&col_identity);
1704:   if (!levels && row_identity && col_identity) {
1705:     /* special case: ilu(0) with natural ordering */
1706:     MatILUFactorSymbolic_SeqAIJ_ilu0(fact,A,isrow,iscol,info);
1707:     if (a->inode.size) {
1708:       fact->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode;
1709:     }
1710:     return(0);
1711:   }

1713:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
1714:   ISGetIndices(isrow,&r);
1715:   ISGetIndices(isicol,&ic);

1717:   /* get new row and diagonal pointers, must be allocated separately because they will be given to the Mat_SeqAIJ and freed separately */
1718:   PetscMalloc1(n+1,&bi);
1719:   PetscMalloc1(n+1,&bdiag);
1720:   bi[0] = bdiag[0] = 0;
1721:   PetscMalloc2(n,&bj_ptr,n,&bjlvl_ptr);

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

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

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

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

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

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

1785:     current_space->array               += nzi;
1786:     current_space->local_used          += nzi;
1787:     current_space->local_remaining     -= nzi;
1788:     current_space_lvl->array           += nzi;
1789:     current_space_lvl->local_used      += nzi;
1790:     current_space_lvl->local_remaining -= nzi;
1791:   }

1793:   ISRestoreIndices(isrow,&r);
1794:   ISRestoreIndices(isicol,&ic);
1795:   /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
1796:   PetscMalloc1(bi[n]+1,&bj);
1797:   PetscFreeSpaceContiguous_LU(&free_space,bj,n,bi,bdiag);

1799:   PetscIncompleteLLDestroy(lnk,lnkbt);
1800:   PetscFreeSpaceDestroy(free_space_lvl);
1801:   PetscFree2(bj_ptr,bjlvl_ptr);

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

1820:   b->free_a       = PETSC_TRUE;
1821:   b->free_ij      = PETSC_TRUE;
1822:   b->singlemalloc = PETSC_FALSE;

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

1826:   b->j    = bj;
1827:   b->i    = bi;
1828:   b->diag = bdiag;
1829:   b->ilen = NULL;
1830:   b->imax = NULL;
1831:   b->row  = isrow;
1832:   b->col  = iscol;
1833:   PetscObjectReference((PetscObject)isrow);
1834:   PetscObjectReference((PetscObject)iscol);
1835:   b->icol = isicol;

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

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

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

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

1878:   f             = info->fill;
1879:   levels        = (PetscInt)info->levels;
1880:   diagonal_fill = (PetscInt)info->diagonal_fill;

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

1884:   ISIdentity(isrow,&row_identity);
1885:   ISIdentity(iscol,&col_identity);
1886:   if (!levels && row_identity && col_identity) { /* special case: ilu(0) with natural ordering */
1887:     MatDuplicateNoCreate_SeqAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_TRUE);

1889:     (fact)->ops->lufactornumeric =  MatLUFactorNumeric_SeqAIJ_inplace;
1890:     if (a->inode.size) {
1891:       (fact)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode_inplace;
1892:     }
1893:     fact->factortype               = MAT_FACTOR_ILU;
1894:     (fact)->info.factor_mallocs    = 0;
1895:     (fact)->info.fill_ratio_given  = info->fill;
1896:     (fact)->info.fill_ratio_needed = 1.0;

1898:     b    = (Mat_SeqAIJ*)(fact)->data;
1899:     b->row  = isrow;
1900:     b->col  = iscol;
1901:     b->icol = isicol;
1902:     PetscMalloc1((fact)->rmap->n+1,&b->solve_work);
1903:     PetscObjectReference((PetscObject)isrow);
1904:     PetscObjectReference((PetscObject)iscol);
1905:     return(0);
1906:   }

1908:   ISGetIndices(isrow,&r);
1909:   ISGetIndices(isicol,&ic);

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

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

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

1922:   /* initial FreeSpace size is f*(ai[n]+1) */
1923:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
1924:   current_space     = free_space;
1925:   PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space_lvl);
1926:   current_space_lvl = free_space_lvl;

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

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

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

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

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

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

1980:     current_space->array               += nzi;
1981:     current_space->local_used          += nzi;
1982:     current_space->local_remaining     -= nzi;
1983:     current_space_lvl->array           += nzi;
1984:     current_space_lvl->local_used      += nzi;
1985:     current_space_lvl->local_remaining -= nzi;
1986:   }

1988:   ISRestoreIndices(isrow,&r);
1989:   ISRestoreIndices(isicol,&ic);

1991:   /* destroy list of free space and other temporary arrays */
1992:   PetscMalloc1(bi[n]+1,&bj);
1993:   PetscFreeSpaceContiguous(&free_space,bj); /* copy free_space -> bj */
1994:   PetscIncompleteLLDestroy(lnk,lnkbt);
1995:   PetscFreeSpaceDestroy(free_space_lvl);
1996:   PetscFree2(bj_ptr,bjlvl_ptr);

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

2011:   /* put together the new matrix */
2012:   MatSeqAIJSetPreallocation_SeqAIJ(fact,MAT_SKIP_ALLOCATION,NULL);
2013:   PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol);
2014:   b    = (Mat_SeqAIJ*)(fact)->data;

2016:   b->free_a       = PETSC_TRUE;
2017:   b->free_ij      = PETSC_TRUE;
2018:   b->singlemalloc = PETSC_FALSE;

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

2038:   (fact)->info.factor_mallocs    = reallocs;
2039:   (fact)->info.fill_ratio_given  = f;
2040:   (fact)->info.fill_ratio_needed = ((PetscReal)bi[n])/((PetscReal)ai[n]);
2041:   (fact)->ops->lufactornumeric   =  MatLUFactorNumeric_SeqAIJ_inplace;
2042:   if (a->inode.size) {
2043:     (fact)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Inode_inplace;
2044:   }
2045:   return(0);
2046: }

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

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

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

2086:   ISGetIndices(ip,&rip);
2087:   ISGetIndices(iip,&riip);

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

2095:   do {
2096:     sctx.newshift = PETSC_FALSE;

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

2101:     for (k = 0; k<mbs; k++) {
2102:       /* zero rtmp */
2103:       nz    = bi[k+1] - bi[k];
2104:       bjtmp = bj + bi[k];
2105:       for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

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

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

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

2127:         /* compute multiplier, update diag(k) and U(i,k) */
2128:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
2129:         uikdi   = -ba[ili]*ba[bdiag[i]]; /* diagonal(k) */
2130:         dk     += uikdi*ba[ili]; /* update diag[k] */
2131:         ba[ili] = uikdi; /* -U(i,k) */

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

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

2156:       /* MatPivotCheck() */
2157:       sctx.rs = rs;
2158:       sctx.pv = dk;
2159:       MatPivotCheck(B,A,info,&sctx,i);
2160:       if (sctx.newshift) break;
2161:       dk = sctx.pv;

2163:       ba[bdiag[k]] = 1.0/dk; /* U(k,k) */
2164:     }
2165:   } while (sctx.newshift);

2167:   PetscFree3(rtmp,il,c2r);
2168:   ISRestoreIndices(ip,&rip);
2169:   ISRestoreIndices(iip,&riip);

2171:   ISIdentity(ip,&perm_identity);
2172:   if (perm_identity) {
2173:     B->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering;
2174:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
2175:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
2176:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;
2177:   } else {
2178:     B->ops->solve          = MatSolve_SeqSBAIJ_1;
2179:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1;
2180:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1;
2181:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1;
2182:   }

2184:   C->assembled    = PETSC_TRUE;
2185:   C->preallocated = PETSC_TRUE;

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

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

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

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

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

2240:   ISGetIndices(ip,&rip);
2241:   ISGetIndices(iip,&riip);

2243:   /* initialization */
2244:   PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&jl);

2246:   do {
2247:     sctx.newshift = PETSC_FALSE;

2249:     for (i=0; i<mbs; i++) jl[i] = mbs;
2250:     il[0] = 0;

2252:     for (k = 0; k<mbs; k++) {
2253:       /* zero rtmp */
2254:       nz    = bi[k+1] - bi[k];
2255:       bjtmp = bj + bi[k];
2256:       for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;

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

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

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

2278:         /* compute multiplier, update diag(k) and U(i,k) */
2279:         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
2280:         uikdi   = -ba[ili]*ba[bi[i]]; /* diagonal(k) */
2281:         dk     += uikdi*ba[ili];
2282:         ba[ili] = uikdi; /* -U(i,k) */

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

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

2305:       sctx.rs = rs;
2306:       sctx.pv = dk;
2307:       MatPivotCheck(B,A,info,&sctx,k);
2308:       if (sctx.newshift) break;
2309:       dk = sctx.pv;

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

2325:   PetscFree3(rtmp,il,jl);
2326:   ISRestoreIndices(ip,&rip);
2327:   ISRestoreIndices(iip,&riip);

2329:   ISIdentity(ip,&perm_identity);
2330:   if (perm_identity) {
2331:     B->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2332:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2333:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2334:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
2335:   } else {
2336:     B->ops->solve          = MatSolve_SeqSBAIJ_1_inplace;
2337:     B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
2338:     B->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_inplace;
2339:     B->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_inplace;
2340:   }

2342:   C->assembled    = PETSC_TRUE;
2343:   C->preallocated = PETSC_TRUE;

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

2356: /*
2357:    icc() under revised new data structure.
2358:    Factored arrays bj and ba are stored as
2359:      U(0,:),...,U(i,:),U(n-1,:)

2361:    ui=fact->i is an array of size n+1, in which
2362:    ui+
2363:      ui[i]:  points to 1st entry of U(i,:),i=0,...,n-1
2364:      ui[n]:  points to U(n-1,n-1)+1

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

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

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

2391:   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);
2392:   MatMissingDiagonal(A,&missing,&d);
2393:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
2394:   ISIdentity(perm,&perm_identity);
2395:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);

2397:   PetscMalloc1(am+1,&ui);
2398:   PetscMalloc1(am+1,&udiag);
2399:   ui[0] = 0;

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

2420:     /* initialization */
2421:     PetscMalloc1(am+1,&ajtmp);

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

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

2434:     /* initial FreeSpace size is fill*(ai[am]+am)/2 */
2435:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,(ai[am]+am)/2),&free_space);
2436:     current_space     = free_space;
2437:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,(ai[am]+am)/2),&free_space_lvl);
2438:     current_space_lvl = free_space_lvl;

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

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

2459:       while (prow < k) {
2460:         nextprow = jl[prow];

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

2473:         /* update il and jl for prow */
2474:         if (jmin < jmax) {
2475:           il[prow] = jmin;
2476:           j        = *cols; jl[prow] = jl[j]; jl[j] = prow;
2477:         }
2478:         prow = nextprow;
2479:       }

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

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

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

2503:       current_space->array           += nzk;
2504:       current_space->local_used      += nzk;
2505:       current_space->local_remaining -= nzk;

2507:       current_space_lvl->array           += nzk;
2508:       current_space_lvl->local_used      += nzk;
2509:       current_space_lvl->local_remaining -= nzk;

2511:       ui[k+1] = ui[k] + nzk;
2512:     }

2514:     ISRestoreIndices(perm,&rip);
2515:     ISRestoreIndices(iperm,&riip);
2516:     PetscFree4(uj_ptr,uj_lvl_ptr,jl,il);
2517:     PetscFree(ajtmp);

2519:     /* copy free_space into uj and free free_space; set ui, uj, udiag in new datastructure; */
2520:     PetscMalloc1(ui[am]+1,&uj);
2521:     PetscFreeSpaceContiguous_Cholesky(&free_space,uj,am,ui,udiag); /* store matrix factor  */
2522:     PetscIncompleteLLDestroy(lnk,lnkbt);
2523:     PetscFreeSpaceDestroy(free_space_lvl);

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

2527:   /* put together the new matrix in MATSEQSBAIJ format */
2528:   b               = (Mat_SeqSBAIJ*)(fact)->data;
2529:   b->singlemalloc = PETSC_FALSE;

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

2533:   b->j             = uj;
2534:   b->i             = ui;
2535:   b->diag          = udiag;
2536:   b->free_diag     = PETSC_TRUE;
2537:   b->ilen          = NULL;
2538:   b->imax          = NULL;
2539:   b->row           = perm;
2540:   b->col           = perm;
2541:   PetscObjectReference((PetscObject)perm);
2542:   PetscObjectReference((PetscObject)perm);
2543:   b->icol          = iperm;
2544:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */

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

2549:   b->maxnz   = b->nz = ui[am];
2550:   b->free_a  = PETSC_TRUE;
2551:   b->free_ij = PETSC_TRUE;

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

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

2593:   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);
2594:   MatMissingDiagonal(A,&missing,&d);
2595:   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
2596:   ISIdentity(perm,&perm_identity);
2597:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);

2599:   PetscMalloc1(am+1,&ui);
2600:   PetscMalloc1(am+1,&udiag);
2601:   ui[0] = 0;

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

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

2621:     /* initialization */
2622:     PetscMalloc1(am+1,&ajtmp);

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

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

2635:     /* initial FreeSpace size is fill*(ai[am]+1) */
2636:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,ai[am]+1),&free_space);
2637:     current_space     = free_space;
2638:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,ai[am]+1),&free_space_lvl);
2639:     current_space_lvl = free_space_lvl;

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

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

2660:       while (prow < k) {
2661:         nextprow = jl[prow];

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

2674:         /* update il and jl for prow */
2675:         if (jmin < jmax) {
2676:           il[prow] = jmin;
2677:           j        = *cols; jl[prow] = jl[j]; jl[j] = prow;
2678:         }
2679:         prow = nextprow;
2680:       }

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

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

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

2704:       current_space->array           += nzk;
2705:       current_space->local_used      += nzk;
2706:       current_space->local_remaining -= nzk;

2708:       current_space_lvl->array           += nzk;
2709:       current_space_lvl->local_used      += nzk;
2710:       current_space_lvl->local_remaining -= nzk;

2712:       ui[k+1] = ui[k] + nzk;
2713:     }

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

2726:     ISRestoreIndices(perm,&rip);
2727:     ISRestoreIndices(iperm,&riip);
2728:     PetscFree4(uj_ptr,uj_lvl_ptr,jl,il);
2729:     PetscFree(ajtmp);

2731:     /* destroy list of free space and other temporary array(s) */
2732:     PetscMalloc1(ui[am]+1,&uj);
2733:     PetscFreeSpaceContiguous(&free_space,uj);
2734:     PetscIncompleteLLDestroy(lnk,lnkbt);
2735:     PetscFreeSpaceDestroy(free_space_lvl);

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

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

2741:   b               = (Mat_SeqSBAIJ*)fact->data;
2742:   b->singlemalloc = PETSC_FALSE;

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

2746:   b->j         = uj;
2747:   b->i         = ui;
2748:   b->diag      = udiag;
2749:   b->free_diag = PETSC_TRUE;
2750:   b->ilen      = NULL;
2751:   b->imax      = NULL;
2752:   b->row       = perm;
2753:   b->col       = perm;

2755:   PetscObjectReference((PetscObject)perm);
2756:   PetscObjectReference((PetscObject)perm);

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

2766:   fact->info.factor_mallocs   = reallocs;
2767:   fact->info.fill_ratio_given = fill;
2768:   if (ai[am] != 0) {
2769:     fact->info.fill_ratio_needed = ((PetscReal)ui[am])/((PetscReal)ai[am]);
2770:   } else {
2771:     fact->info.fill_ratio_needed = 0.0;
2772:   }
2773:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_inplace;
2774:   return(0);
2775: }

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

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

2797:   /* check whether perm is the identity mapping */
2798:   ISIdentity(perm,&perm_identity);
2799:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);
2800:   ISGetIndices(iperm,&riip);
2801:   ISGetIndices(perm,&rip);

2803:   /* initialization */
2804:   PetscMalloc1(am+1,&ui);
2805:   PetscMalloc1(am+1,&udiag);
2806:   ui[0] = 0;

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

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

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

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

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

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

2852:       /* update il and jl for prow */
2853:       if (jmin < jmax) {
2854:         il[prow] = jmin;
2855:         j        = *uj_ptr;
2856:         jl[prow] = jl[j];
2857:         jl[j]    = prow;
2858:       }
2859:       prow = nextprow;
2860:     }

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

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

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

2881:     current_space->array           += nzk;
2882:     current_space->local_used      += nzk;
2883:     current_space->local_remaining -= nzk;

2885:     ui[k+1] = ui[k] + nzk;
2886:   }

2888:   ISRestoreIndices(perm,&rip);
2889:   ISRestoreIndices(iperm,&riip);
2890:   PetscFree4(ui_ptr,jl,il,cols);

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

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

2899:   b               = (Mat_SeqSBAIJ*)fact->data;
2900:   b->singlemalloc = PETSC_FALSE;
2901:   b->free_a       = PETSC_TRUE;
2902:   b->free_ij      = PETSC_TRUE;

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

2906:   b->j         = uj;
2907:   b->i         = ui;
2908:   b->diag      = udiag;
2909:   b->free_diag = PETSC_TRUE;
2910:   b->ilen      = NULL;
2911:   b->imax      = NULL;
2912:   b->row       = perm;
2913:   b->col       = perm;

2915:   PetscObjectReference((PetscObject)perm);
2916:   PetscObjectReference((PetscObject)perm);

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

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

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

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

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

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

2968:   /* check whether perm is the identity mapping */
2969:   ISIdentity(perm,&perm_identity);
2970:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);
2971:   ISGetIndices(iperm,&riip);
2972:   ISGetIndices(perm,&rip);

2974:   /* initialization */
2975:   PetscMalloc1(am+1,&ui);
2976:   ui[0] = 0;

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

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

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

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

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

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

3022:       /* update il and jl for prow */
3023:       if (jmin < jmax) {
3024:         il[prow] = jmin;
3025:         j        = *uj_ptr; jl[prow] = jl[j]; jl[j] = prow;
3026:       }
3027:       prow = nextprow;
3028:     }

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

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

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

3049:     current_space->array           += nzk;
3050:     current_space->local_used      += nzk;
3051:     current_space->local_remaining -= nzk;

3053:     ui[k+1] = ui[k] + nzk;
3054:   }

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

3067:   ISRestoreIndices(perm,&rip);
3068:   ISRestoreIndices(iperm,&riip);
3069:   PetscFree4(ui_ptr,jl,il,cols);

3071:   /* destroy list of free space and other temporary array(s) */
3072:   PetscMalloc1(ui[am]+1,&uj);
3073:   PetscFreeSpaceContiguous(&free_space,uj);
3074:   PetscLLDestroy(lnk,lnkbt);

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

3078:   b               = (Mat_SeqSBAIJ*)fact->data;
3079:   b->singlemalloc = PETSC_FALSE;
3080:   b->free_a       = PETSC_TRUE;
3081:   b->free_ij      = PETSC_TRUE;

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

3085:   b->j    = uj;
3086:   b->i    = ui;
3087:   b->diag = NULL;
3088:   b->ilen = NULL;
3089:   b->imax = NULL;
3090:   b->row  = perm;
3091:   b->col  = perm;

3093:   PetscObjectReference((PetscObject)perm);
3094:   PetscObjectReference((PetscObject)perm);

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

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

3103:   fact->info.factor_mallocs   = reallocs;
3104:   fact->info.fill_ratio_given = fill;
3105:   if (ai[am] != 0) {
3106:     fact->info.fill_ratio_needed = ((PetscReal)ui[am])/((PetscReal)ai[am]);
3107:   } else {
3108:     fact->info.fill_ratio_needed = 0.0;
3109:   }
3110:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_inplace;
3111:   return(0);
3112: }

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

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

3128:   VecGetArrayRead(bb,&b);
3129:   VecGetArrayWrite(xx,&x);

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

3144:   /* backward solve the upper triangular */
3145:   for (i=n-1; i>=0; i--) {
3146:     v   = aa + adiag[i+1] + 1;
3147:     vi  = aj + adiag[i+1] + 1;
3148:     nz  = adiag[i] - adiag[i+1]-1;
3149:     sum = x[i];
3150:     PetscSparseDenseMinusDot(sum,x,v,vi,nz);
3151:     x[i] = sum*v[nz]; /* x[i]=aa[adiag[i]]*sum; v++; */
3152:   }

3154:   PetscLogFlops(2.0*a->nz - A->cmap->n);
3155:   VecRestoreArrayRead(bb,&b);
3156:   VecRestoreArrayWrite(xx,&x);
3157:   return(0);
3158: }

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

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

3174:   VecGetArrayRead(bb,&b);
3175:   VecGetArrayWrite(xx,&x);
3176:   tmp  = a->solve_work;

3178:   ISGetIndices(isrow,&rout); r = rout;
3179:   ISGetIndices(iscol,&cout); c = cout;

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

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

3203:   ISRestoreIndices(isrow,&rout);
3204:   ISRestoreIndices(iscol,&cout);
3205:   VecRestoreArrayRead(bb,&b);
3206:   VecRestoreArrayWrite(xx,&x);
3207:   PetscLogFlops(2.0*a->nz - A->cmap->n);
3208:   return(0);
3209: }

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

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

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

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

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

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

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

3256:   PetscMalloc1(nnz_max+1,&bj);
3257:   PetscMalloc1(nnz_max+1,&ba);

3259:   /* put together the new matrix */
3260:   MatSeqAIJSetPreallocation_SeqAIJ(B,MAT_SKIP_ALLOCATION,NULL);
3261:   PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);
3262:   b    = (Mat_SeqAIJ*)B->data;

3264:   b->free_a       = PETSC_TRUE;
3265:   b->free_ij      = PETSC_TRUE;
3266:   b->singlemalloc = PETSC_FALSE;

3268:   b->a    = ba;
3269:   b->j    = bj;
3270:   b->i    = bi;
3271:   b->diag = bdiag;
3272:   b->ilen = NULL;
3273:   b->imax = NULL;
3274:   b->row  = isrow;
3275:   b->col  = iscol;
3276:   PetscObjectReference((PetscObject)isrow);
3277:   PetscObjectReference((PetscObject)iscol);
3278:   b->icol = isicol;

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

3284:   B->factortype            = MAT_FACTOR_ILUDT;
3285:   B->info.factor_mallocs   = 0;
3286:   B->info.fill_ratio_given = ((PetscReal)nnz_max)/((PetscReal)ai[n]);
3287:   /* ------- end of symbolic factorization ---------*/

3289:   ISGetIndices(isrow,&r);
3290:   ISGetIndices(isicol,&ic);
3291:   ics  = ic;

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

3297:   /* im: used by PetscLLAddSortedLU(); jtmp: working array for column indices of active row */
3298:   PetscMalloc2(n,&im,n,&jtmp);
3299:   /* rtmp, vtmp: working arrays for sparse and contiguous row entries of active row */
3300:   PetscMalloc2(n,&rtmp,n,&vtmp);
3301:   PetscArrayzero(rtmp,n);

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

3316:     /* load in initial (unfactored row) */
3317:     aatmp = a->a + ai[r[i]];
3318:     for (j=0; j<nzi; j++) {
3319:       rtmp[ics[*ajtmp++]] = *aatmp++;
3320:     }

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

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

3335:     /* numerical factorization */
3336:     bjtmp = jtmp;
3337:     row   = *bjtmp++; /* 1st pivot row */
3338:     while (row < i) {
3339:       pc         = rtmp + row;
3340:       pv         = ba + bdiag[row]; /* 1./(diag of the pivot row) */
3341:       multiplier = (*pc) * (*pv);
3342:       *pc        = multiplier;
3343:       if (PetscAbsScalar(*pc) > dt) { /* apply tolerance dropping rule */
3344:         pj = bj + bdiag[row+1] + 1;         /* point to 1st entry of U(row,:) */
3345:         pv = ba + bdiag[row+1] + 1;
3346:         nz = bdiag[row] - bdiag[row+1] - 1;         /* num of entries in U(row,:), excluding diagonal */
3347:         for (j=0; j<nz; j++) rtmp[*pj++] -= multiplier * (*pv++);
3348:         PetscLogFlops(1+2.0*nz);
3349:       }
3350:       row = *bjtmp++;
3351:     }

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

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

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

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

3406:     bjtmp = bj + bdiag[i+1]+1;
3407:     batmp = ba + bdiag[i+1]+1;
3408:     for (k=0; k<ncut; k++) {
3409:       bjtmp[k] = jtmp[nzi_bl+1+k];
3410:       batmp[k] = vtmp[nzi_bl+1+k];
3411:     }

3413:     im[i] = nzi;   /* used by PetscLLAddSortedLU() */
3414:   } /* for (i=0; i<n; i++) */
3415:   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]);

3417:   ISRestoreIndices(isrow,&r);
3418:   ISRestoreIndices(isicol,&ic);

3420:   PetscLLDestroy(lnk,lnkbt);
3421:   PetscFree2(im,jtmp);
3422:   PetscFree2(rtmp,vtmp);
3423:   PetscFree(bdiag_rev);

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

3428:   ISIdentity(isrow,&row_identity);
3429:   ISIdentity(isicol,&icol_identity);
3430:   if (row_identity && icol_identity) {
3431:     B->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
3432:   } else {
3433:     B->ops->solve = MatSolve_SeqAIJ;
3434:   }

3436:   B->ops->solveadd          = NULL;
3437:   B->ops->solvetranspose    = NULL;
3438:   B->ops->solvetransposeadd = NULL;
3439:   B->ops->matsolve          = NULL;
3440:   B->assembled              = PETSC_TRUE;
3441:   B->preallocated           = PETSC_TRUE;
3442:   return(0);
3443: }

3445: /* a wraper of MatILUDTFactor_SeqAIJ() */
3446: /*
3447:     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
3448: */

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

3455:   MatILUDTFactor_SeqAIJ(A,row,col,info,&fact);
3456:   return(0);
3457: }

3459: /*
3460:    same as MatLUFactorNumeric_SeqAIJ(), except using contiguous array matrix factors
3461:    - intend to replace existing MatLUFactorNumeric_SeqAIJ()
3462: */
3463: /*
3464:     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
3465: */

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

3481:   ISGetIndices(isrow,&r);
3482:   ISGetIndices(isicol,&ic);
3483:   PetscMalloc1(n+1,&rtmp);
3484:   ics  = ic;

3486:   for (i=0; i<n; i++) {
3487:     /* initialize rtmp array */
3488:     nzl   = bi[i+1] - bi[i];       /* num of nozeros in L(i,:) */
3489:     bjtmp = bj + bi[i];
3490:     for  (j=0; j<nzl; j++) rtmp[*bjtmp++] = 0.0;
3491:     rtmp[i] = 0.0;
3492:     nzu     = bdiag[i] - bdiag[i+1]; /* num of nozeros in U(i,:) */
3493:     bjtmp   = bj + bdiag[i+1] + 1;
3494:     for  (j=0; j<nzu; j++) rtmp[*bjtmp++] = 0.0;

3496:     /* load in initial unfactored row of A */
3497:     nz    = ai[r[i]+1] - ai[r[i]];
3498:     ajtmp = aj + ai[r[i]];
3499:     v     = aa + ai[r[i]];
3500:     for (j=0; j<nz; j++) {
3501:       rtmp[ics[*ajtmp++]] = v[j];
3502:     }

3504:     /* numerical factorization */
3505:     bjtmp = bj + bi[i]; /* point to 1st entry of L(i,:) */
3506:     nzl   = bi[i+1] - bi[i]; /* num of entries in L(i,:) */
3507:     k     = 0;
3508:     while (k < nzl) {
3509:       row        = *bjtmp++;
3510:       pc         = rtmp + row;
3511:       pv         = b->a + bdiag[row]; /* 1./(diag of the pivot row) */
3512:       multiplier = (*pc) * (*pv);
3513:       *pc        = multiplier;
3514:       if (PetscAbsScalar(multiplier) > dt) {
3515:         pj = bj + bdiag[row+1] + 1;         /* point to 1st entry of U(row,:) */
3516:         pv = b->a + bdiag[row+1] + 1;
3517:         nz = bdiag[row] - bdiag[row+1] - 1;         /* num of entries in U(row,:), excluding diagonal */
3518:         for (j=0; j<nz; j++) rtmp[*pj++] -= multiplier * (*pv++);
3519:         PetscLogFlops(1+2.0*nz);
3520:       }
3521:       k++;
3522:     }

3524:     /* finished row so stick it into b->a */
3525:     /* L-part */
3526:     pv  = b->a + bi[i];
3527:     pj  = bj + bi[i];
3528:     nzl = bi[i+1] - bi[i];
3529:     for (j=0; j<nzl; j++) {
3530:       pv[j] = rtmp[pj[j]];
3531:     }

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

3537:     /* U-part */
3538:     pv  = b->a + bdiag[i+1] + 1;
3539:     pj  = bj + bdiag[i+1] + 1;
3540:     nzu = bdiag[i] - bdiag[i+1] - 1;
3541:     for (j=0; j<nzu; j++) {
3542:       pv[j] = rtmp[pj[j]];
3543:     }
3544:   }

3546:   PetscFree(rtmp);
3547:   ISRestoreIndices(isicol,&ic);
3548:   ISRestoreIndices(isrow,&r);

3550:   ISIdentity(isrow,&row_identity);
3551:   ISIdentity(isicol,&col_identity);
3552:   if (row_identity && col_identity) {
3553:     C->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
3554:   } else {
3555:     C->ops->solve = MatSolve_SeqAIJ;
3556:   }
3557:   C->ops->solveadd          = NULL;
3558:   C->ops->solvetranspose    = NULL;
3559:   C->ops->solvetransposeadd = NULL;
3560:   C->ops->matsolve          = NULL;
3561:   C->assembled              = PETSC_TRUE;
3562:   C->preallocated           = PETSC_TRUE;

3564:   PetscLogFlops(C->cmap->n);
3565:   return(0);
3566: }