xref: /petsc/src/mat/impls/aij/seq/aijfact.c (revision 207126cb90d7f2cb72e2179925bb8d6adc164b40) 
1 #define PETSCMAT_DLL
2 
3 #include "src/mat/impls/aij/seq/aij.h"
4 #include "src/inline/dot.h"
5 #include "src/inline/spops.h"
6 #include "petscbt.h"
7 #include "src/mat/utils/freespace.h"
8 
9 #undef __FUNCT__
10 #define __FUNCT__ "MatOrdering_Flow_SeqAIJ"
11 PetscErrorCode MatOrdering_Flow_SeqAIJ(Mat mat,const MatOrderingType type,IS *irow,IS *icol)
12 {
13   PetscFunctionBegin;
14 
15   SETERRQ(PETSC_ERR_SUP,"Code not written");
16 #if !defined(PETSC_USE_DEBUG)
17   PetscFunctionReturn(0);
18 #endif
19 }
20 
21 
22 #if !defined(PETSC_AVOID_GNUCOPYRIGHT_CODE)
23 EXTERN PetscErrorCode SPARSEKIT2dperm(PetscInt*,MatScalar*,PetscInt*,PetscInt*,MatScalar*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
24 EXTERN PetscErrorCode SPARSEKIT2ilutp(PetscInt*,MatScalar*,PetscInt*,PetscInt*,PetscInt*,PetscReal,PetscReal*,PetscInt*,MatScalar*,PetscInt*,PetscInt*,PetscInt*,MatScalar*,PetscInt*,PetscInt*,PetscErrorCode*);
25 EXTERN PetscErrorCode SPARSEKIT2msrcsr(PetscInt*,MatScalar*,PetscInt*,MatScalar*,PetscInt*,PetscInt*,MatScalar*,PetscInt*);
26 #endif
27 
28 #undef __FUNCT__
29 #define __FUNCT__ "MatILUDTFactor_SeqAIJ"
30   /* ------------------------------------------------------------
31 
32           This interface was contribed by Tony Caola
33 
34      This routine is an interface to the pivoting drop-tolerance
35      ILU routine written by Yousef Saad (saad@cs.umn.edu) as part of
36      SPARSEKIT2.
37 
38      The SPARSEKIT2 routines used here are covered by the GNU
39      copyright; see the file gnu in this directory.
40 
41      Thanks to Prof. Saad, Dr. Hysom, and Dr. Smith for their
42      help in getting this routine ironed out.
43 
44      The major drawback to this routine is that if info->fill is
45      not large enough it fails rather than allocating more space;
46      this can be fixed by hacking/improving the f2c version of
47      Yousef Saad's code.
48 
49      ------------------------------------------------------------
50 */
51 PetscErrorCode MatILUDTFactor_SeqAIJ(Mat A,IS isrow,IS iscol,MatFactorInfo *info,Mat *fact)
52 {
53 #if defined(PETSC_AVOID_GNUCOPYRIGHT_CODE)
54   PetscFunctionBegin;
55   SETERRQ(PETSC_ERR_SUP_SYS,"This distribution does not include GNU Copyright code\n\
56   You can obtain the drop tolerance routines by installing PETSc from\n\
57   www.mcs.anl.gov/petsc\n");
58 #else
59   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data,*b;
60   IS             iscolf,isicol,isirow;
61   PetscTruth     reorder;
62   PetscErrorCode ierr,sierr;
63   PetscInt       *c,*r,*ic,i,n = A->rmap.n;
64   PetscInt       *old_i = a->i,*old_j = a->j,*new_i,*old_i2 = 0,*old_j2 = 0,*new_j;
65   PetscInt       *ordcol,*iwk,*iperm,*jw;
66   PetscInt       jmax,lfill,job,*o_i,*o_j;
67   MatScalar      *old_a = a->a,*w,*new_a,*old_a2 = 0,*wk,*o_a;
68   PetscReal      af;
69 
70   PetscFunctionBegin;
71 
72   if (info->dt == PETSC_DEFAULT)      info->dt      = .005;
73   if (info->dtcount == PETSC_DEFAULT) info->dtcount = (PetscInt)(1.5*a->rmax);
74   if (info->dtcol == PETSC_DEFAULT)   info->dtcol   = .01;
75   if (info->fill == PETSC_DEFAULT)    info->fill    = ((double)(n*(info->dtcount+1)))/a->nz;
76   lfill   = (PetscInt)(info->dtcount/2.0);
77   jmax    = (PetscInt)(info->fill*a->nz);
78 
79 
80   /* ------------------------------------------------------------
81      If reorder=.TRUE., then the original matrix has to be
82      reordered to reflect the user selected ordering scheme, and
83      then de-reordered so it is in it's original format.
84      Because Saad's dperm() is NOT in place, we have to copy
85      the original matrix and allocate more storage. . .
86      ------------------------------------------------------------
87   */
88 
89   /* set reorder to true if either isrow or iscol is not identity */
90   ierr = ISIdentity(isrow,&reorder);CHKERRQ(ierr);
91   if (reorder) {ierr = ISIdentity(iscol,&reorder);CHKERRQ(ierr);}
92   reorder = PetscNot(reorder);
93 
94 
95   /* storage for ilu factor */
96   ierr = PetscMalloc((n+1)*sizeof(PetscInt),&new_i);CHKERRQ(ierr);
97   ierr = PetscMalloc(jmax*sizeof(PetscInt),&new_j);CHKERRQ(ierr);
98   ierr = PetscMalloc(jmax*sizeof(MatScalar),&new_a);CHKERRQ(ierr);
99   ierr = PetscMalloc(n*sizeof(PetscInt),&ordcol);CHKERRQ(ierr);
100 
101   /* ------------------------------------------------------------
102      Make sure that everything is Fortran formatted (1-Based)
103      ------------------------------------------------------------
104   */
105   for (i=old_i[0];i<old_i[n];i++) {
106     old_j[i]++;
107   }
108   for(i=0;i<n+1;i++) {
109     old_i[i]++;
110   };
111 
112 
113   if (reorder) {
114     ierr = ISGetIndices(iscol,&c);CHKERRQ(ierr);
115     ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
116     for(i=0;i<n;i++) {
117       r[i]  = r[i]+1;
118       c[i]  = c[i]+1;
119     }
120     ierr = PetscMalloc((n+1)*sizeof(PetscInt),&old_i2);CHKERRQ(ierr);
121     ierr = PetscMalloc((old_i[n]-old_i[0]+1)*sizeof(PetscInt),&old_j2);CHKERRQ(ierr);
122     ierr = PetscMalloc((old_i[n]-old_i[0]+1)*sizeof(MatScalar),&old_a2);CHKERRQ(ierr);
123     job  = 3; SPARSEKIT2dperm(&n,old_a,old_j,old_i,old_a2,old_j2,old_i2,r,c,&job);
124     for (i=0;i<n;i++) {
125       r[i]  = r[i]-1;
126       c[i]  = c[i]-1;
127     }
128     ierr = ISRestoreIndices(iscol,&c);CHKERRQ(ierr);
129     ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
130     o_a = old_a2;
131     o_j = old_j2;
132     o_i = old_i2;
133   } else {
134     o_a = old_a;
135     o_j = old_j;
136     o_i = old_i;
137   }
138 
139   /* ------------------------------------------------------------
140      Call Saad's ilutp() routine to generate the factorization
141      ------------------------------------------------------------
142   */
143 
144   ierr = PetscMalloc(2*n*sizeof(PetscInt),&iperm);CHKERRQ(ierr);
145   ierr = PetscMalloc(2*n*sizeof(PetscInt),&jw);CHKERRQ(ierr);
146   ierr = PetscMalloc(n*sizeof(PetscScalar),&w);CHKERRQ(ierr);
147 
148   SPARSEKIT2ilutp(&n,o_a,o_j,o_i,&lfill,(PetscReal)info->dt,&info->dtcol,&n,new_a,new_j,new_i,&jmax,w,jw,iperm,&sierr);
149   if (sierr) {
150     switch (sierr) {
151       case -3: SETERRQ2(PETSC_ERR_LIB,"ilutp(), matrix U overflows, need larger info->fill current fill %G space allocated %D",info->fill,jmax);
152       case -2: SETERRQ2(PETSC_ERR_LIB,"ilutp(), matrix L overflows, need larger info->fill current fill %G space allocated %D",info->fill,jmax);
153       case -5: SETERRQ(PETSC_ERR_LIB,"ilutp(), zero row encountered");
154       case -1: SETERRQ(PETSC_ERR_LIB,"ilutp(), input matrix may be wrong");
155       case -4: SETERRQ1(PETSC_ERR_LIB,"ilutp(), illegal info->fill value %D",jmax);
156       default: SETERRQ1(PETSC_ERR_LIB,"ilutp(), zero pivot detected on row %D",sierr);
157     }
158   }
159 
160   ierr = PetscFree(w);CHKERRQ(ierr);
161   ierr = PetscFree(jw);CHKERRQ(ierr);
162 
163   /* ------------------------------------------------------------
164      Saad's routine gives the result in Modified Sparse Row (msr)
165      Convert to Compressed Sparse Row format (csr)
166      ------------------------------------------------------------
167   */
168 
169   ierr = PetscMalloc(n*sizeof(PetscScalar),&wk);CHKERRQ(ierr);
170   ierr = PetscMalloc((n+1)*sizeof(PetscInt),&iwk);CHKERRQ(ierr);
171 
172   SPARSEKIT2msrcsr(&n,new_a,new_j,new_a,new_j,new_i,wk,iwk);
173 
174   ierr = PetscFree(iwk);CHKERRQ(ierr);
175   ierr = PetscFree(wk);CHKERRQ(ierr);
176 
177   if (reorder) {
178     ierr = PetscFree(old_a2);CHKERRQ(ierr);
179     ierr = PetscFree(old_j2);CHKERRQ(ierr);
180     ierr = PetscFree(old_i2);CHKERRQ(ierr);
181   } else {
182     /* fix permutation of old_j that the factorization introduced */
183     for (i=old_i[0]; i<old_i[n]; i++) {
184       old_j[i-1] = iperm[old_j[i-1]-1];
185     }
186   }
187 
188   /* get rid of the shift to indices starting at 1 */
189   for (i=0; i<n+1; i++) {
190     old_i[i]--;
191   }
192   for (i=old_i[0];i<old_i[n];i++) {
193     old_j[i]--;
194   }
195 
196   /* Make the factored matrix 0-based */
197   for (i=0; i<n+1; i++) {
198     new_i[i]--;
199   }
200   for (i=new_i[0];i<new_i[n];i++) {
201     new_j[i]--;
202   }
203 
204   /*-- due to the pivoting, we need to reorder iscol to correctly --*/
205   /*-- permute the right-hand-side and solution vectors           --*/
206   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
207   ierr = ISInvertPermutation(isrow,PETSC_DECIDE,&isirow);CHKERRQ(ierr);
208   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
209   for(i=0; i<n; i++) {
210     ordcol[i] = ic[iperm[i]-1];
211   };
212   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
213   ierr = ISDestroy(isicol);CHKERRQ(ierr);
214 
215   ierr = PetscFree(iperm);CHKERRQ(ierr);
216 
217   ierr = ISCreateGeneral(PETSC_COMM_SELF,n,ordcol,&iscolf);CHKERRQ(ierr);
218   ierr = PetscFree(ordcol);CHKERRQ(ierr);
219 
220   /*----- put together the new matrix -----*/
221 
222   ierr = MatCreate(((PetscObject)A)->comm,fact);CHKERRQ(ierr);
223   ierr = MatSetSizes(*fact,n,n,n,n);CHKERRQ(ierr);
224   ierr = MatSetType(*fact,((PetscObject)A)->type_name);CHKERRQ(ierr);
225   ierr = MatSeqAIJSetPreallocation_SeqAIJ(*fact,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr);
226   (*fact)->factor    = MAT_FACTOR_LU;
227   (*fact)->assembled = PETSC_TRUE;
228 
229   b = (Mat_SeqAIJ*)(*fact)->data;
230   b->free_a        = PETSC_TRUE;
231   b->free_ij       = PETSC_TRUE;
232   b->singlemalloc  = PETSC_FALSE;
233   b->a             = new_a;
234   b->j             = new_j;
235   b->i             = new_i;
236   b->ilen          = 0;
237   b->imax          = 0;
238   /*  I am not sure why these are the inverses of the row and column permutations; but the other way is NO GOOD */
239   b->row           = isirow;
240   b->col           = iscolf;
241   ierr = PetscMalloc((n+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
242   b->maxnz = b->nz = new_i[n];
243   ierr = MatMarkDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
244   (*fact)->info.factor_mallocs = 0;
245 
246   af = ((double)b->nz)/((double)a->nz) + .001;
247   ierr = PetscInfo2(A,"Fill ratio:given %G needed %G\n",info->fill,af);CHKERRQ(ierr);
248   ierr = PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);CHKERRQ(ierr);
249   ierr = PetscInfo1(A,"PCFactorSetFill(pc,%G);\n",af);CHKERRQ(ierr);
250   ierr = PetscInfo(A,"for best performance.\n");CHKERRQ(ierr);
251 
252   ierr = MatILUDTFactor_Inode(A,isrow,iscol,info,fact);CHKERRQ(ierr);
253 
254   PetscFunctionReturn(0);
255 #endif
256 }
257 
258 #undef __FUNCT__
259 #define __FUNCT__ "MatGetFactor_seqaij_petsc"
260 PetscErrorCode MatGetFactor_seqaij_petsc(Mat A,MatFactorType ftype,Mat *B)
261 {
262   PetscInt           n = A->rmap.n;
263   PetscErrorCode     ierr;
264 
265   PetscFunctionBegin;
266   ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr);
267   ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr);
268   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
269     ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr);
270   } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
271     ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr);
272     ierr = MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr);
273     (*B)->ops->iccfactorsymbolic      = MatICCFactorSymbolic_SeqAIJ;
274     (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqAIJ;
275     (*B)->ops->choleskyfactornumeric  = MatCholeskyFactorNumeric_SeqAIJ;
276   } else SETERRQ(PETSC_ERR_SUP,"Factor type not supported");
277   PetscFunctionReturn(0);
278 }
279 
280 #undef __FUNCT__
281 #define __FUNCT__ "MatLUFactorSymbolic_SeqAIJ"
282 PetscErrorCode MatLUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatFactorInfo *info,Mat *B)
283 {
284   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data,*b;
285   IS                 isicol;
286   PetscErrorCode     ierr;
287   PetscInt           *r,*ic,i,n=A->rmap.n,*ai=a->i,*aj=a->j;
288   PetscInt           *bi,*bj,*ajtmp;
289   PetscInt           *bdiag,row,nnz,nzi,reallocs=0,nzbd,*im;
290   PetscReal          f;
291   PetscInt           nlnk,*lnk,k,**bi_ptr;
292   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
293   PetscBT            lnkbt;
294 
295   PetscFunctionBegin;
296   if (A->rmap.N != A->cmap.N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
297   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
298   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
299   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
300 
301   /* get new row pointers */
302   ierr = PetscMalloc((n+1)*sizeof(PetscInt),&bi);CHKERRQ(ierr);
303   bi[0] = 0;
304 
305   /* bdiag is location of diagonal in factor */
306   ierr = PetscMalloc((n+1)*sizeof(PetscInt),&bdiag);CHKERRQ(ierr);
307   bdiag[0] = 0;
308 
309   /* linked list for storing column indices of the active row */
310   nlnk = n + 1;
311   ierr = PetscLLCreate(n,n,nlnk,lnk,lnkbt);CHKERRQ(ierr);
312 
313   ierr = PetscMalloc2(n+1,PetscInt**,&bi_ptr,n+1,PetscInt,&im);CHKERRQ(ierr);
314 
315   /* initial FreeSpace size is f*(ai[n]+1) */
316   f = info->fill;
317   ierr = PetscFreeSpaceGet((PetscInt)(f*(ai[n]+1)),&free_space);CHKERRQ(ierr);
318   current_space = free_space;
319 
320   for (i=0; i<n; i++) {
321     /* copy previous fill into linked list */
322     nzi = 0;
323     nnz = ai[r[i]+1] - ai[r[i]];
324     if (!nnz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
325     ajtmp = aj + ai[r[i]];
326     ierr = PetscLLAddPerm(nnz,ajtmp,ic,n,nlnk,lnk,lnkbt);CHKERRQ(ierr);
327     nzi += nlnk;
328 
329     /* add pivot rows into linked list */
330     row = lnk[n];
331     while (row < i) {
332       nzbd    = bdiag[row] - bi[row] + 1; /* num of entries in the row with column index <= row */
333       ajtmp   = bi_ptr[row] + nzbd; /* points to the entry next to the diagonal */
334       ierr = PetscLLAddSortedLU(ajtmp,row,nlnk,lnk,lnkbt,i,nzbd,im);CHKERRQ(ierr);
335       nzi += nlnk;
336       row  = lnk[row];
337     }
338     bi[i+1] = bi[i] + nzi;
339     im[i]   = nzi;
340 
341     /* mark bdiag */
342     nzbd = 0;
343     nnz  = nzi;
344     k    = lnk[n];
345     while (nnz-- && k < i){
346       nzbd++;
347       k = lnk[k];
348     }
349     bdiag[i] = bi[i] + nzbd;
350 
351     /* if free space is not available, make more free space */
352     if (current_space->local_remaining<nzi) {
353       nnz = (n - i)*nzi; /* estimated and max additional space needed */
354       ierr = PetscFreeSpaceGet(nnz,&current_space);CHKERRQ(ierr);
355       reallocs++;
356     }
357 
358     /* copy data into free space, then initialize lnk */
359     ierr = PetscLLClean(n,n,nzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
360     bi_ptr[i] = current_space->array;
361     current_space->array           += nzi;
362     current_space->local_used      += nzi;
363     current_space->local_remaining -= nzi;
364   }
365 #if defined(PETSC_USE_INFO)
366   if (ai[n] != 0) {
367     PetscReal af = ((PetscReal)bi[n])/((PetscReal)ai[n]);
368     ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,f,af);CHKERRQ(ierr);
369     ierr = PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);CHKERRQ(ierr);
370     ierr = PetscInfo1(A,"PCFactorSetFill(pc,%G);\n",af);CHKERRQ(ierr);
371     ierr = PetscInfo(A,"for best performance.\n");CHKERRQ(ierr);
372   } else {
373     ierr = PetscInfo(A,"Empty matrix\n");CHKERRQ(ierr);
374   }
375 #endif
376 
377   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
378   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
379 
380   /* destroy list of free space and other temporary array(s) */
381   ierr = PetscMalloc((bi[n]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr);
382   ierr = PetscFreeSpaceContiguous(&free_space,bj);CHKERRQ(ierr);
383   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
384   ierr = PetscFree2(bi_ptr,im);CHKERRQ(ierr);
385 
386   /* put together the new matrix */
387   ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr);
388   ierr = PetscLogObjectParent(*B,isicol);CHKERRQ(ierr);
389   b    = (Mat_SeqAIJ*)(*B)->data;
390   b->free_a       = PETSC_TRUE;
391   b->free_ij      = PETSC_TRUE;
392   b->singlemalloc = PETSC_FALSE;
393   ierr          = PetscMalloc((bi[n]+1)*sizeof(PetscScalar),&b->a);CHKERRQ(ierr);
394   b->j          = bj;
395   b->i          = bi;
396   b->diag       = bdiag;
397   b->ilen       = 0;
398   b->imax       = 0;
399   b->row        = isrow;
400   b->col        = iscol;
401   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
402   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
403   b->icol       = isicol;
404   ierr          = PetscMalloc((n+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
405 
406   /* In b structure:  Free imax, ilen, old a, old j.  Allocate solve_work, new a, new j */
407   ierr = PetscLogObjectMemory(*B,(bi[n]-n)*(sizeof(PetscInt)+sizeof(PetscScalar)));CHKERRQ(ierr);
408   b->maxnz = b->nz = bi[n] ;
409 
410   (*B)->factor                 =  MAT_FACTOR_LU;
411   (*B)->info.factor_mallocs    = reallocs;
412   (*B)->info.fill_ratio_given  = f;
413 
414   if (ai[n] != 0) {
415     (*B)->info.fill_ratio_needed = ((PetscReal)bi[n])/((PetscReal)ai[n]);
416   } else {
417     (*B)->info.fill_ratio_needed = 0.0;
418   }
419   ierr = MatLUFactorSymbolic_Inode(A,isrow,iscol,info,B);CHKERRQ(ierr);
420   (*B)->ops->lufactornumeric =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
421   PetscFunctionReturn(0);
422 }
423 
424 /*
425     Trouble in factorization, should we dump the original matrix?
426 */
427 #undef __FUNCT__
428 #define __FUNCT__ "MatFactorDumpMatrix"
429 PetscErrorCode MatFactorDumpMatrix(Mat A)
430 {
431   PetscErrorCode ierr;
432   PetscTruth     flg;
433 
434   PetscFunctionBegin;
435   ierr = PetscOptionsHasName(PETSC_NULL,"-mat_factor_dump_on_error",&flg);CHKERRQ(ierr);
436   if (flg) {
437     PetscViewer viewer;
438     char        filename[PETSC_MAX_PATH_LEN];
439 
440     ierr = PetscSNPrintf(filename,PETSC_MAX_PATH_LEN,"matrix_factor_error.%d",PetscGlobalRank);CHKERRQ(ierr);
441     ierr = PetscViewerBinaryOpen(((PetscObject)A)->comm,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
442     ierr = MatView(A,viewer);CHKERRQ(ierr);
443     ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
444   }
445   PetscFunctionReturn(0);
446 }
447 
448 /* ----------------------------------------------------------- */
449 #undef __FUNCT__
450 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJ"
451 PetscErrorCode MatLUFactorNumeric_SeqAIJ(Mat A,MatFactorInfo *info,Mat *B)
452 {
453   Mat            C=*B;
454   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ *)C->data;
455   IS             isrow = b->row,isicol = b->icol;
456   PetscErrorCode ierr;
457   PetscInt       *r,*ic,i,j,n=A->rmap.n,*bi=b->i,*bj=b->j;
458   PetscInt       *ajtmp,*bjtmp,nz,row,*ics;
459   PetscInt       *diag_offset = b->diag,diag,*pj;
460   PetscScalar    *rtmp,*pc,multiplier,*rtmps;
461   MatScalar      *v,*pv;
462   PetscScalar    d;
463   PetscReal      rs;
464   LUShift_Ctx    sctx;
465   PetscInt       newshift,*ddiag;
466 
467   PetscFunctionBegin;
468   ierr  = ISGetIndices(isrow,&r);CHKERRQ(ierr);
469   ierr  = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
470   ierr  = PetscMalloc((n+1)*sizeof(PetscScalar),&rtmp);CHKERRQ(ierr);
471   ierr  = PetscMemzero(rtmp,(n+1)*sizeof(PetscScalar));CHKERRQ(ierr);
472   rtmps = rtmp; ics = ic;
473 
474   sctx.shift_top  = 0;
475   sctx.nshift_max = 0;
476   sctx.shift_lo   = 0;
477   sctx.shift_hi   = 0;
478 
479   /* if both shift schemes are chosen by user, only use info->shiftpd */
480   if (info->shiftpd && info->shiftnz) info->shiftnz = 0.0;
481   if (info->shiftpd) { /* set sctx.shift_top=max{rs} */
482     PetscInt *aai = a->i;
483     ddiag          = a->diag;
484     sctx.shift_top = 0;
485     for (i=0; i<n; i++) {
486       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
487       d  = (a->a)[ddiag[i]];
488       rs = -PetscAbsScalar(d) - PetscRealPart(d);
489       v  = a->a+aai[i];
490       nz = aai[i+1] - aai[i];
491       for (j=0; j<nz; j++)
492 	rs += PetscAbsScalar(v[j]);
493       if (rs>sctx.shift_top) sctx.shift_top = rs;
494     }
495     if (sctx.shift_top < info->zeropivot) sctx.shift_top = info->zeropivot;
496     sctx.shift_top    *= 1.1;
497     sctx.nshift_max   = 5;
498     sctx.shift_lo     = 0.;
499     sctx.shift_hi     = 1.;
500   }
501 
502   sctx.shift_amount = 0;
503   sctx.nshift       = 0;
504   do {
505     sctx.lushift = PETSC_FALSE;
506     for (i=0; i<n; i++){
507       nz    = bi[i+1] - bi[i];
508       bjtmp = bj + bi[i];
509       for  (j=0; j<nz; j++) rtmps[bjtmp[j]] = 0.0;
510 
511       /* load in initial (unfactored row) */
512       nz    = a->i[r[i]+1] - a->i[r[i]];
513       ajtmp = a->j + a->i[r[i]];
514       v     = a->a + a->i[r[i]];
515       for (j=0; j<nz; j++) {
516         rtmp[ics[ajtmp[j]]] = v[j];
517       }
518       rtmp[ics[r[i]]] += sctx.shift_amount; /* shift the diagonal of the matrix */
519 
520       row = *bjtmp++;
521       while  (row < i) {
522         pc = rtmp + row;
523         if (*pc != 0.0) {
524           pv         = b->a + diag_offset[row];
525           pj         = b->j + diag_offset[row] + 1;
526           multiplier = *pc / *pv++;
527           *pc        = multiplier;
528           nz         = bi[row+1] - diag_offset[row] - 1;
529           for (j=0; j<nz; j++) rtmps[pj[j]] -= multiplier * pv[j];
530           ierr = PetscLogFlops(2*nz);CHKERRQ(ierr);
531         }
532         row = *bjtmp++;
533       }
534       /* finished row so stick it into b->a */
535       pv   = b->a + bi[i] ;
536       pj   = b->j + bi[i] ;
537       nz   = bi[i+1] - bi[i];
538       diag = diag_offset[i] - bi[i];
539       rs   = 0.0;
540       for (j=0; j<nz; j++) {
541         pv[j] = rtmps[pj[j]];
542         if (j != diag) rs += PetscAbsScalar(pv[j]);
543       }
544 
545       /* 9/13/02 Victor Eijkhout suggested scaling zeropivot by rs for matrices with funny scalings */
546       sctx.rs  = rs;
547       sctx.pv  = pv[diag];
548       ierr = MatLUCheckShift_inline(info,sctx,i,newshift);CHKERRQ(ierr);
549       if (newshift == 1) break;
550     }
551 
552     if (info->shiftpd && !sctx.lushift && info->shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
553       /*
554        * if no shift in this attempt & shifting & started shifting & can refine,
555        * then try lower shift
556        */
557       sctx.shift_hi        = info->shift_fraction;
558       info->shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
559       sctx.shift_amount    = info->shift_fraction * sctx.shift_top;
560       sctx.lushift         = PETSC_TRUE;
561       sctx.nshift++;
562     }
563   } while (sctx.lushift);
564 
565   /* invert diagonal entries for simplier triangular solves */
566   for (i=0; i<n; i++) {
567     b->a[diag_offset[i]] = 1.0/b->a[diag_offset[i]];
568   }
569   ierr = PetscFree(rtmp);CHKERRQ(ierr);
570   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
571   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
572   C->factor       = MAT_FACTOR_LU;
573   C->assembled    = PETSC_TRUE;
574   C->preallocated = PETSC_TRUE;
575   ierr = PetscLogFlops(C->cmap.n);CHKERRQ(ierr);
576   if (sctx.nshift){
577     if (info->shiftnz) {
578       ierr = PetscInfo2(A,"number of shift_nz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);CHKERRQ(ierr);
579     } else if (info->shiftpd) {
580       ierr = PetscInfo4(A,"number of shift_pd tries %D, shift_amount %G, diagonal shifted up by %e fraction top_value %e\n",sctx.nshift,sctx.shift_amount,info->shift_fraction,sctx.shift_top);CHKERRQ(ierr);
581     }
582   }
583   PetscFunctionReturn(0);
584 }
585 
586 /*
587    This routine implements inplace ILU(0) with row or/and column permutations.
588    Input:
589      A - original matrix
590    Output;
591      A - a->i (rowptr) is same as original rowptr, but factored i-the row is stored in rowperm[i]
592          a->j (col index) is permuted by the inverse of colperm, then sorted
593          a->a reordered accordingly with a->j
594          a->diag (ptr to diagonal elements) is updated.
595 */
596 #undef __FUNCT__
597 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJ_InplaceWithPerm"
598 PetscErrorCode MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(Mat A,MatFactorInfo *info,Mat *B)
599 {
600   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data;
601   IS             isrow = a->row,isicol = a->icol;
602   PetscErrorCode ierr;
603   PetscInt       *r,*ic,i,j,n=A->rmap.n,*ai=a->i,*aj=a->j;
604   PetscInt       *ajtmp,nz,row,*ics;
605   PetscInt       *diag = a->diag,nbdiag,*pj;
606   PetscScalar    *rtmp,*pc,multiplier,d;
607   MatScalar      *v,*pv;
608   PetscReal      rs;
609   LUShift_Ctx    sctx;
610   PetscInt       newshift;
611 
612   PetscFunctionBegin;
613   if (A != *B) SETERRQ(PETSC_ERR_ARG_INCOMP,"input and output matrix must have same address");
614   ierr  = ISGetIndices(isrow,&r);CHKERRQ(ierr);
615   ierr  = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
616   ierr  = PetscMalloc((n+1)*sizeof(PetscScalar),&rtmp);CHKERRQ(ierr);
617   ierr  = PetscMemzero(rtmp,(n+1)*sizeof(PetscScalar));CHKERRQ(ierr);
618   ics = ic;
619 
620   sctx.shift_top  = 0;
621   sctx.nshift_max = 0;
622   sctx.shift_lo   = 0;
623   sctx.shift_hi   = 0;
624 
625   /* if both shift schemes are chosen by user, only use info->shiftpd */
626   if (info->shiftpd && info->shiftnz) info->shiftnz = 0.0;
627   if (info->shiftpd) { /* set sctx.shift_top=max{rs} */
628     sctx.shift_top = 0;
629     for (i=0; i<n; i++) {
630       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
631       d  = (a->a)[diag[i]];
632       rs = -PetscAbsScalar(d) - PetscRealPart(d);
633       v  = a->a+ai[i];
634       nz = ai[i+1] - ai[i];
635       for (j=0; j<nz; j++)
636 	rs += PetscAbsScalar(v[j]);
637       if (rs>sctx.shift_top) sctx.shift_top = rs;
638     }
639     if (sctx.shift_top < info->zeropivot) sctx.shift_top = info->zeropivot;
640     sctx.shift_top    *= 1.1;
641     sctx.nshift_max   = 5;
642     sctx.shift_lo     = 0.;
643     sctx.shift_hi     = 1.;
644   }
645 
646   sctx.shift_amount = 0;
647   sctx.nshift       = 0;
648   do {
649     sctx.lushift = PETSC_FALSE;
650     for (i=0; i<n; i++){
651       /* load in initial unfactored row */
652       nz    = ai[r[i]+1] - ai[r[i]];
653       ajtmp = aj + ai[r[i]];
654       v     = a->a + ai[r[i]];
655       /* sort permuted ajtmp and values v accordingly */
656       for (j=0; j<nz; j++) ajtmp[j] = ics[ajtmp[j]];
657       ierr = PetscSortIntWithScalarArray(nz,ajtmp,v);CHKERRQ(ierr);
658 
659       diag[r[i]] = ai[r[i]];
660       for (j=0; j<nz; j++) {
661         rtmp[ajtmp[j]] = v[j];
662         if (ajtmp[j] < i) diag[r[i]]++; /* update a->diag */
663       }
664       rtmp[r[i]] += sctx.shift_amount; /* shift the diagonal of the matrix */
665 
666       row = *ajtmp++;
667       while  (row < i) {
668         pc = rtmp + row;
669         if (*pc != 0.0) {
670           pv         = a->a + diag[r[row]];
671           pj         = aj + diag[r[row]] + 1;
672 
673           multiplier = *pc / *pv++;
674           *pc        = multiplier;
675           nz         = ai[r[row]+1] - diag[r[row]] - 1;
676           for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
677           ierr = PetscLogFlops(2*nz);CHKERRQ(ierr);
678         }
679         row = *ajtmp++;
680       }
681       /* finished row so overwrite it onto a->a */
682       pv   = a->a + ai[r[i]] ;
683       pj   = aj + ai[r[i]] ;
684       nz   = ai[r[i]+1] - ai[r[i]];
685       nbdiag = diag[r[i]] - ai[r[i]]; /* num of entries before the diagonal */
686 
687       rs   = 0.0;
688       for (j=0; j<nz; j++) {
689         pv[j] = rtmp[pj[j]];
690         if (j != nbdiag) rs += PetscAbsScalar(pv[j]);
691       }
692 
693       /* 9/13/02 Victor Eijkhout suggested scaling zeropivot by rs for matrices with funny scalings */
694       sctx.rs  = rs;
695       sctx.pv  = pv[nbdiag];
696       ierr = MatLUCheckShift_inline(info,sctx,i,newshift);CHKERRQ(ierr);
697       if (newshift == 1) break;
698     }
699 
700     if (info->shiftpd && !sctx.lushift && info->shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
701       /*
702        * if no shift in this attempt & shifting & started shifting & can refine,
703        * then try lower shift
704        */
705       sctx.shift_hi        = info->shift_fraction;
706       info->shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
707       sctx.shift_amount    = info->shift_fraction * sctx.shift_top;
708       sctx.lushift         = PETSC_TRUE;
709       sctx.nshift++;
710     }
711   } while (sctx.lushift);
712 
713   /* invert diagonal entries for simplier triangular solves */
714   for (i=0; i<n; i++) {
715     a->a[diag[r[i]]] = 1.0/a->a[diag[r[i]]];
716   }
717 
718   ierr = PetscFree(rtmp);CHKERRQ(ierr);
719   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
720   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
721   A->factor     = MAT_FACTOR_LU;
722   A->ops->solve = MatSolve_SeqAIJ_InplaceWithPerm;
723   A->assembled = PETSC_TRUE;
724   A->preallocated = PETSC_TRUE;
725   ierr = PetscLogFlops(A->cmap.n);CHKERRQ(ierr);
726   if (sctx.nshift){
727     if (info->shiftnz) {
728       ierr = PetscInfo2(A,"number of shift_nz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);CHKERRQ(ierr);
729     } else if (info->shiftpd) {
730       ierr = PetscInfo4(A,"number of shift_pd tries %D, shift_amount %G, diagonal shifted up by %e fraction top_value %e\n",sctx.nshift,sctx.shift_amount,info->shift_fraction,sctx.shift_top);CHKERRQ(ierr);
731     }
732   }
733   PetscFunctionReturn(0);
734 }
735 
736 /* ----------------------------------------------------------- */
737 #undef __FUNCT__
738 #define __FUNCT__ "MatLUFactor_SeqAIJ"
739 PetscErrorCode MatLUFactor_SeqAIJ(Mat A,IS row,IS col,MatFactorInfo *info)
740 {
741   PetscErrorCode ierr;
742   Mat            C;
743 
744   PetscFunctionBegin;
745   ierr = MatGetFactor(A,"petsc",MAT_FACTOR_LU,&C);CHKERRQ(ierr);
746   ierr = MatLUFactorSymbolic(A,row,col,info,&C);CHKERRQ(ierr);
747   ierr = MatLUFactorNumeric(A,info,&C);CHKERRQ(ierr);
748   ierr = MatHeaderCopy(A,C);CHKERRQ(ierr);
749   ierr = PetscLogObjectParent(A,((Mat_SeqAIJ*)(A->data))->icol);CHKERRQ(ierr);
750   PetscFunctionReturn(0);
751 }
752 /* ----------------------------------------------------------- */
753 #undef __FUNCT__
754 #define __FUNCT__ "MatSolve_SeqAIJ"
755 PetscErrorCode MatSolve_SeqAIJ(Mat A,Vec bb,Vec xx)
756 {
757   Mat_SeqAIJ        *a = (Mat_SeqAIJ*)A->data;
758   IS                iscol = a->col,isrow = a->row;
759   PetscErrorCode    ierr;
760   PetscInt          *r,*c,i, n = A->rmap.n,*vi,*ai = a->i,*aj = a->j;
761   PetscInt          nz,*rout,*cout;
762   PetscScalar       *x,*tmp,*tmps,sum;
763   const PetscScalar *b;
764   const MatScalar   *aa = a->a,*v;
765 
766   PetscFunctionBegin;
767   if (!n) PetscFunctionReturn(0);
768 
769   ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);
770   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
771   tmp  = a->solve_work;
772 
773   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
774   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
775 
776   /* forward solve the lower triangular */
777   tmp[0] = b[*r++];
778   tmps   = tmp;
779   for (i=1; i<n; i++) {
780     v   = aa + ai[i] ;
781     vi  = aj + ai[i] ;
782     nz  = a->diag[i] - ai[i];
783     sum = b[*r++];
784     SPARSEDENSEMDOT(sum,tmps,v,vi,nz);
785     tmp[i] = sum;
786   }
787 
788   /* backward solve the upper triangular */
789   for (i=n-1; i>=0; i--){
790     v   = aa + a->diag[i] + 1;
791     vi  = aj + a->diag[i] + 1;
792     nz  = ai[i+1] - a->diag[i] - 1;
793     sum = tmp[i];
794     SPARSEDENSEMDOT(sum,tmps,v,vi,nz);
795     x[*c--] = tmp[i] = sum*aa[a->diag[i]];
796   }
797 
798   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
799   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
800   ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);
801   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
802   ierr = PetscLogFlops(2*a->nz - A->cmap.n);CHKERRQ(ierr);
803   PetscFunctionReturn(0);
804 }
805 
806 #undef __FUNCT__
807 #define __FUNCT__ "MatMatSolve_SeqAIJ"
808 PetscErrorCode MatMatSolve_SeqAIJ(Mat A,Mat B,Mat X)
809 {
810   Mat_SeqAIJ      *a = (Mat_SeqAIJ*)A->data;
811   IS              iscol = a->col,isrow = a->row;
812   PetscErrorCode  ierr;
813   PetscInt        *r,*c,i, n = A->rmap.n,*vi,*ai = a->i,*aj = a->j;
814   PetscInt        nz,*rout,*cout,neq;
815   PetscScalar     *x,*b,*tmp,*tmps,sum;
816   const MatScalar *aa = a->a,*v;
817 
818   PetscFunctionBegin;
819   if (!n) PetscFunctionReturn(0);
820 
821   ierr = MatGetArray(B,&b);CHKERRQ(ierr);
822   ierr = MatGetArray(X,&x);CHKERRQ(ierr);
823 
824   tmp  = a->solve_work;
825   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
826   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
827 
828   for (neq=0; neq<n; neq++){
829     /* forward solve the lower triangular */
830     tmp[0] = b[r[0]];
831     tmps   = tmp;
832     for (i=1; i<n; i++) {
833       v   = aa + ai[i] ;
834       vi  = aj + ai[i] ;
835       nz  = a->diag[i] - ai[i];
836       sum = b[r[i]];
837       SPARSEDENSEMDOT(sum,tmps,v,vi,nz);
838       tmp[i] = sum;
839     }
840     /* backward solve the upper triangular */
841     for (i=n-1; i>=0; i--){
842       v   = aa + a->diag[i] + 1;
843       vi  = aj + a->diag[i] + 1;
844       nz  = ai[i+1] - a->diag[i] - 1;
845       sum = tmp[i];
846       SPARSEDENSEMDOT(sum,tmps,v,vi,nz);
847       x[c[i]] = tmp[i] = sum*aa[a->diag[i]];
848     }
849 
850     b += n;
851     x += n;
852   }
853   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
854   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
855   ierr = MatRestoreArray(B,&b);CHKERRQ(ierr);
856   ierr = MatRestoreArray(X,&x);CHKERRQ(ierr);
857   ierr = PetscLogFlops(n*(2*a->nz - n));CHKERRQ(ierr);
858   PetscFunctionReturn(0);
859 }
860 
861 #undef __FUNCT__
862 #define __FUNCT__ "MatSolve_SeqAIJ_InplaceWithPerm"
863 PetscErrorCode MatSolve_SeqAIJ_InplaceWithPerm(Mat A,Vec bb,Vec xx)
864 {
865   Mat_SeqAIJ      *a = (Mat_SeqAIJ*)A->data;
866   IS              iscol = a->col,isrow = a->row;
867   PetscErrorCode  ierr;
868   PetscInt        *r,*c,i, n = A->rmap.n,*vi,*ai = a->i,*aj = a->j;
869   PetscInt        nz,*rout,*cout,row;
870   PetscScalar     *x,*b,*tmp,*tmps,sum;
871   const MatScalar *aa = a->a,*v;
872 
873   PetscFunctionBegin;
874   if (!n) PetscFunctionReturn(0);
875 
876   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
877   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
878   tmp  = a->solve_work;
879 
880   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
881   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
882 
883   /* forward solve the lower triangular */
884   tmp[0] = b[*r++];
885   tmps   = tmp;
886   for (row=1; row<n; row++) {
887     i   = rout[row]; /* permuted row */
888     v   = aa + ai[i] ;
889     vi  = aj + ai[i] ;
890     nz  = a->diag[i] - ai[i];
891     sum = b[*r++];
892     SPARSEDENSEMDOT(sum,tmps,v,vi,nz);
893     tmp[row] = sum;
894   }
895 
896   /* backward solve the upper triangular */
897   for (row=n-1; row>=0; row--){
898     i   = rout[row]; /* permuted row */
899     v   = aa + a->diag[i] + 1;
900     vi  = aj + a->diag[i] + 1;
901     nz  = ai[i+1] - a->diag[i] - 1;
902     sum = tmp[row];
903     SPARSEDENSEMDOT(sum,tmps,v,vi,nz);
904     x[*c--] = tmp[row] = sum*aa[a->diag[i]];
905   }
906 
907   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
908   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
909   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
910   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
911   ierr = PetscLogFlops(2*a->nz - A->cmap.n);CHKERRQ(ierr);
912   PetscFunctionReturn(0);
913 }
914 
915 /* ----------------------------------------------------------- */
916 #undef __FUNCT__
917 #define __FUNCT__ "MatSolve_SeqAIJ_NaturalOrdering"
918 PetscErrorCode MatSolve_SeqAIJ_NaturalOrdering(Mat A,Vec bb,Vec xx)
919 {
920   Mat_SeqAIJ        *a = (Mat_SeqAIJ*)A->data;
921   PetscErrorCode    ierr;
922   PetscInt          n = A->rmap.n;
923   const PetscInt    *ai = a->i,*aj = a->j,*adiag = a->diag,*vi;
924   PetscScalar       *x;
925   const PetscScalar *b;
926   const MatScalar   *aa = a->a;
927 #if !defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
928   PetscInt          adiag_i,i,nz,ai_i;
929   const MatScalar   *v;
930   PetscScalar       sum;
931 #endif
932 
933   PetscFunctionBegin;
934   if (!n) PetscFunctionReturn(0);
935 
936   ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);
937   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
938 
939 #if defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
940   fortransolveaij_(&n,x,ai,aj,adiag,aa,b);
941 #else
942   /* forward solve the lower triangular */
943   x[0] = b[0];
944   for (i=1; i<n; i++) {
945     ai_i = ai[i];
946     v    = aa + ai_i;
947     vi   = aj + ai_i;
948     nz   = adiag[i] - ai_i;
949     sum  = b[i];
950     while (nz--) sum -= *v++ * x[*vi++];
951     x[i] = sum;
952   }
953 
954   /* backward solve the upper triangular */
955   for (i=n-1; i>=0; i--){
956     adiag_i = adiag[i];
957     v       = aa + adiag_i + 1;
958     vi      = aj + adiag_i + 1;
959     nz      = ai[i+1] - adiag_i - 1;
960     sum     = x[i];
961     while (nz--) sum -= *v++ * x[*vi++];
962     x[i]    = sum*aa[adiag_i];
963   }
964 #endif
965   ierr = PetscLogFlops(2*a->nz - A->cmap.n);CHKERRQ(ierr);
966   ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);
967   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
968   PetscFunctionReturn(0);
969 }
970 
971 #undef __FUNCT__
972 #define __FUNCT__ "MatSolveAdd_SeqAIJ"
973 PetscErrorCode MatSolveAdd_SeqAIJ(Mat A,Vec bb,Vec yy,Vec xx)
974 {
975   Mat_SeqAIJ      *a = (Mat_SeqAIJ*)A->data;
976   IS              iscol = a->col,isrow = a->row;
977   PetscErrorCode  ierr;
978   PetscInt        *r,*c,i, n = A->rmap.n,*vi,*ai = a->i,*aj = a->j;
979   PetscInt        nz,*rout,*cout;
980   PetscScalar     *x,*b,*tmp,sum;
981   const MatScalar *aa = a->a,*v;
982 
983   PetscFunctionBegin;
984   if (yy != xx) {ierr = VecCopy(yy,xx);CHKERRQ(ierr);}
985 
986   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
987   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
988   tmp  = a->solve_work;
989 
990   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
991   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
992 
993   /* forward solve the lower triangular */
994   tmp[0] = b[*r++];
995   for (i=1; i<n; i++) {
996     v   = aa + ai[i] ;
997     vi  = aj + ai[i] ;
998     nz  = a->diag[i] - ai[i];
999     sum = b[*r++];
1000     while (nz--) sum -= *v++ * tmp[*vi++ ];
1001     tmp[i] = sum;
1002   }
1003 
1004   /* backward solve the upper triangular */
1005   for (i=n-1; i>=0; i--){
1006     v   = aa + a->diag[i] + 1;
1007     vi  = aj + a->diag[i] + 1;
1008     nz  = ai[i+1] - a->diag[i] - 1;
1009     sum = tmp[i];
1010     while (nz--) sum -= *v++ * tmp[*vi++ ];
1011     tmp[i] = sum*aa[a->diag[i]];
1012     x[*c--] += tmp[i];
1013   }
1014 
1015   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
1016   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
1017   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
1018   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
1019   ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr);
1020 
1021   PetscFunctionReturn(0);
1022 }
1023 /* -------------------------------------------------------------------*/
1024 #undef __FUNCT__
1025 #define __FUNCT__ "MatSolveTranspose_SeqAIJ"
1026 PetscErrorCode MatSolveTranspose_SeqAIJ(Mat A,Vec bb,Vec xx)
1027 {
1028   Mat_SeqAIJ      *a = (Mat_SeqAIJ*)A->data;
1029   IS              iscol = a->col,isrow = a->row;
1030   PetscErrorCode  ierr;
1031   PetscInt        *r,*c,i,n = A->rmap.n,*vi,*ai = a->i,*aj = a->j;
1032   PetscInt        nz,*rout,*cout,*diag = a->diag;
1033   PetscScalar     *x,*b,*tmp,s1;
1034   const MatScalar *aa = a->a,*v;
1035 
1036   PetscFunctionBegin;
1037   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
1038   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
1039   tmp  = a->solve_work;
1040 
1041   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
1042   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
1043 
1044   /* copy the b into temp work space according to permutation */
1045   for (i=0; i<n; i++) tmp[i] = b[c[i]];
1046 
1047   /* forward solve the U^T */
1048   for (i=0; i<n; i++) {
1049     v   = aa + diag[i] ;
1050     vi  = aj + diag[i] + 1;
1051     nz  = ai[i+1] - diag[i] - 1;
1052     s1  = tmp[i];
1053     s1 *= (*v++);  /* multiply by inverse of diagonal entry */
1054     while (nz--) {
1055       tmp[*vi++ ] -= (*v++)*s1;
1056     }
1057     tmp[i] = s1;
1058   }
1059 
1060   /* backward solve the L^T */
1061   for (i=n-1; i>=0; i--){
1062     v   = aa + diag[i] - 1 ;
1063     vi  = aj + diag[i] - 1 ;
1064     nz  = diag[i] - ai[i];
1065     s1  = tmp[i];
1066     while (nz--) {
1067       tmp[*vi-- ] -= (*v--)*s1;
1068     }
1069   }
1070 
1071   /* copy tmp into x according to permutation */
1072   for (i=0; i<n; i++) x[r[i]] = tmp[i];
1073 
1074   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
1075   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
1076   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
1077   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
1078 
1079   ierr = PetscLogFlops(2*a->nz-A->cmap.n);CHKERRQ(ierr);
1080   PetscFunctionReturn(0);
1081 }
1082 
1083 #undef __FUNCT__
1084 #define __FUNCT__ "MatSolveTransposeAdd_SeqAIJ"
1085 PetscErrorCode MatSolveTransposeAdd_SeqAIJ(Mat A,Vec bb,Vec zz,Vec xx)
1086 {
1087   Mat_SeqAIJ      *a = (Mat_SeqAIJ*)A->data;
1088   IS              iscol = a->col,isrow = a->row;
1089   PetscErrorCode  ierr;
1090   PetscInt        *r,*c,i,n = A->rmap.n,*vi,*ai = a->i,*aj = a->j;
1091   PetscInt        nz,*rout,*cout,*diag = a->diag;
1092   PetscScalar     *x,*b,*tmp;
1093   const MatScalar *aa = a->a,*v;
1094 
1095   PetscFunctionBegin;
1096   if (zz != xx) {ierr = VecCopy(zz,xx);CHKERRQ(ierr);}
1097 
1098   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
1099   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
1100   tmp = a->solve_work;
1101 
1102   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
1103   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
1104 
1105   /* copy the b into temp work space according to permutation */
1106   for (i=0; i<n; i++) tmp[i] = b[c[i]];
1107 
1108   /* forward solve the U^T */
1109   for (i=0; i<n; i++) {
1110     v   = aa + diag[i] ;
1111     vi  = aj + diag[i] + 1;
1112     nz  = ai[i+1] - diag[i] - 1;
1113     tmp[i] *= *v++;
1114     while (nz--) {
1115       tmp[*vi++ ] -= (*v++)*tmp[i];
1116     }
1117   }
1118 
1119   /* backward solve the L^T */
1120   for (i=n-1; i>=0; i--){
1121     v   = aa + diag[i] - 1 ;
1122     vi  = aj + diag[i] - 1 ;
1123     nz  = diag[i] - ai[i];
1124     while (nz--) {
1125       tmp[*vi-- ] -= (*v--)*tmp[i];
1126     }
1127   }
1128 
1129   /* copy tmp into x according to permutation */
1130   for (i=0; i<n; i++) x[r[i]] += tmp[i];
1131 
1132   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
1133   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
1134   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
1135   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
1136 
1137   ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr);
1138   PetscFunctionReturn(0);
1139 }
1140 /* ----------------------------------------------------------------*/
1141 EXTERN PetscErrorCode Mat_CheckInode(Mat,PetscTruth);
1142 EXTERN PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat,MatDuplicateOption,Mat*);
1143 
1144 #undef __FUNCT__
1145 #define __FUNCT__ "MatILUFactorSymbolic_SeqAIJ"
1146 PetscErrorCode MatILUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatFactorInfo *info,Mat *fact)
1147 {
1148   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data,*b;
1149   IS                 isicol;
1150   PetscErrorCode     ierr;
1151   PetscInt           *r,*ic,n=A->rmap.n,*ai=a->i,*aj=a->j,d;
1152   PetscInt           *bi,*cols,nnz,*cols_lvl;
1153   PetscInt           *bdiag,prow,fm,nzbd,len, reallocs=0,dcount=0;
1154   PetscInt           i,levels,diagonal_fill;
1155   PetscTruth         col_identity,row_identity;
1156   PetscReal          f;
1157   PetscInt           nlnk,*lnk,*lnk_lvl=PETSC_NULL;
1158   PetscBT            lnkbt;
1159   PetscInt           nzi,*bj,**bj_ptr,**bjlvl_ptr;
1160   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
1161   PetscFreeSpaceList free_space_lvl=PETSC_NULL,current_space_lvl=PETSC_NULL;
1162   PetscTruth         missing;
1163 
1164   PetscFunctionBegin;
1165   if (A->rmap.n != A->cmap.n) SETERRQ2(PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %D columns %D",A->rmap.n,A->cmap.n);
1166   f             = info->fill;
1167   levels        = (PetscInt)info->levels;
1168   diagonal_fill = (PetscInt)info->diagonal_fill;
1169   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
1170 
1171   /* special case that simply copies fill pattern */
1172   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
1173   ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
1174   if (!levels && row_identity && col_identity) {
1175     ierr = MatDuplicateNoCreate_SeqAIJ(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr);
1176     (*fact)->factor                 = MAT_FACTOR_LU;
1177     (*fact)->info.factor_mallocs    = 0;
1178     (*fact)->info.fill_ratio_given  = info->fill;
1179     (*fact)->info.fill_ratio_needed = 1.0;
1180     b               = (Mat_SeqAIJ*)(*fact)->data;
1181     ierr = MatMissingDiagonal(A,&missing,&d);CHKERRQ(ierr);
1182     if (missing) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
1183     b->row              = isrow;
1184     b->col              = iscol;
1185     b->icol             = isicol;
1186     ierr                = PetscMalloc(((*fact)->rmap.n+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
1187     (*fact)->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
1188     ierr                = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
1189     ierr                = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
1190     ierr = Mat_CheckInode(*fact,PETSC_FALSE);CHKERRQ(ierr);
1191     PetscFunctionReturn(0);
1192   }
1193 
1194   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
1195   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
1196 
1197   /* get new row pointers */
1198   ierr = PetscMalloc((n+1)*sizeof(PetscInt),&bi);CHKERRQ(ierr);
1199   bi[0] = 0;
1200   /* bdiag is location of diagonal in factor */
1201   ierr = PetscMalloc((n+1)*sizeof(PetscInt),&bdiag);CHKERRQ(ierr);
1202   bdiag[0]  = 0;
1203 
1204   ierr = PetscMalloc((2*n+1)*sizeof(PetscInt**),&bj_ptr);CHKERRQ(ierr);
1205   bjlvl_ptr = (PetscInt**)(bj_ptr + n);
1206 
1207   /* create a linked list for storing column indices of the active row */
1208   nlnk = n + 1;
1209   ierr = PetscIncompleteLLCreate(n,n,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr);
1210 
1211   /* initial FreeSpace size is f*(ai[n]+1) */
1212   ierr = PetscFreeSpaceGet((PetscInt)(f*(ai[n]+1)),&free_space);CHKERRQ(ierr);
1213   current_space = free_space;
1214   ierr = PetscFreeSpaceGet((PetscInt)(f*(ai[n]+1)),&free_space_lvl);CHKERRQ(ierr);
1215   current_space_lvl = free_space_lvl;
1216 
1217   for (i=0; i<n; i++) {
1218     nzi = 0;
1219     /* copy current row into linked list */
1220     nnz  = ai[r[i]+1] - ai[r[i]];
1221     if (!nnz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
1222     cols = aj + ai[r[i]];
1223     lnk[i] = -1; /* marker to indicate if diagonal exists */
1224     ierr = PetscIncompleteLLInit(nnz,cols,n,ic,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr);
1225     nzi += nlnk;
1226 
1227     /* make sure diagonal entry is included */
1228     if (diagonal_fill && lnk[i] == -1) {
1229       fm = n;
1230       while (lnk[fm] < i) fm = lnk[fm];
1231       lnk[i]     = lnk[fm]; /* insert diagonal into linked list */
1232       lnk[fm]    = i;
1233       lnk_lvl[i] = 0;
1234       nzi++; dcount++;
1235     }
1236 
1237     /* add pivot rows into the active row */
1238     nzbd = 0;
1239     prow = lnk[n];
1240     while (prow < i) {
1241       nnz      = bdiag[prow];
1242       cols     = bj_ptr[prow] + nnz + 1;
1243       cols_lvl = bjlvl_ptr[prow] + nnz + 1;
1244       nnz      = bi[prow+1] - bi[prow] - nnz - 1;
1245       ierr = PetscILULLAddSorted(nnz,cols,levels,cols_lvl,prow,nlnk,lnk,lnk_lvl,lnkbt,prow);CHKERRQ(ierr);
1246       nzi += nlnk;
1247       prow = lnk[prow];
1248       nzbd++;
1249     }
1250     bdiag[i] = nzbd;
1251     bi[i+1]  = bi[i] + nzi;
1252 
1253     /* if free space is not available, make more free space */
1254     if (current_space->local_remaining<nzi) {
1255       nnz = nzi*(n - i); /* estimated and max additional space needed */
1256       ierr = PetscFreeSpaceGet(nnz,&current_space);CHKERRQ(ierr);
1257       ierr = PetscFreeSpaceGet(nnz,&current_space_lvl);CHKERRQ(ierr);
1258       reallocs++;
1259     }
1260 
1261     /* copy data into free_space and free_space_lvl, then initialize lnk */
1262     ierr = PetscIncompleteLLClean(n,n,nzi,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);CHKERRQ(ierr);
1263     bj_ptr[i]    = current_space->array;
1264     bjlvl_ptr[i] = current_space_lvl->array;
1265 
1266     /* make sure the active row i has diagonal entry */
1267     if (*(bj_ptr[i]+bdiag[i]) != i) {
1268       SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,"Row %D has missing diagonal in factored matrix\n\
1269     try running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill",i);
1270     }
1271 
1272     current_space->array           += nzi;
1273     current_space->local_used      += nzi;
1274     current_space->local_remaining -= nzi;
1275     current_space_lvl->array           += nzi;
1276     current_space_lvl->local_used      += nzi;
1277     current_space_lvl->local_remaining -= nzi;
1278   }
1279 
1280   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
1281   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
1282 
1283   /* destroy list of free space and other temporary arrays */
1284   ierr = PetscMalloc((bi[n]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr);
1285   ierr = PetscFreeSpaceContiguous(&free_space,bj);CHKERRQ(ierr);
1286   ierr = PetscIncompleteLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
1287   ierr = PetscFreeSpaceDestroy(free_space_lvl);CHKERRQ(ierr);
1288   ierr = PetscFree(bj_ptr);CHKERRQ(ierr);
1289 
1290 #if defined(PETSC_USE_INFO)
1291   {
1292     PetscReal af = ((PetscReal)bi[n])/((PetscReal)ai[n]);
1293     ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,f,af);CHKERRQ(ierr);
1294     ierr = PetscInfo1(A,"Run with -[sub_]pc_factor_fill %G or use \n",af);CHKERRQ(ierr);
1295     ierr = PetscInfo1(A,"PCFactorSetFill([sub]pc,%G);\n",af);CHKERRQ(ierr);
1296     ierr = PetscInfo(A,"for best performance.\n");CHKERRQ(ierr);
1297     if (diagonal_fill) {
1298       ierr = PetscInfo1(A,"Detected and replaced %D missing diagonals",dcount);CHKERRQ(ierr);
1299     }
1300   }
1301 #endif
1302 
1303   /* put together the new matrix */
1304   ierr = PetscLogObjectParent(*fact,isicol);CHKERRQ(ierr);
1305   b = (Mat_SeqAIJ*)(*fact)->data;
1306   b->free_a       = PETSC_TRUE;
1307   b->free_ij      = PETSC_TRUE;
1308   b->singlemalloc = PETSC_FALSE;
1309   len = (bi[n] )*sizeof(PetscScalar);
1310   ierr = PetscMalloc(len+1,&b->a);CHKERRQ(ierr);
1311   b->j          = bj;
1312   b->i          = bi;
1313   for (i=0; i<n; i++) bdiag[i] += bi[i];
1314   b->diag       = bdiag;
1315   b->ilen       = 0;
1316   b->imax       = 0;
1317   b->row        = isrow;
1318   b->col        = iscol;
1319   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
1320   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
1321   b->icol       = isicol;
1322   ierr = PetscMalloc((n+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
1323   /* In b structure:  Free imax, ilen, old a, old j.
1324      Allocate bdiag, solve_work, new a, new j */
1325   ierr = PetscLogObjectMemory(*fact,(bi[n]-n) * (sizeof(PetscInt)+sizeof(PetscScalar)));CHKERRQ(ierr);
1326   b->maxnz             = b->nz = bi[n] ;
1327   (*fact)->factor = MAT_FACTOR_LU;
1328   (*fact)->info.factor_mallocs    = reallocs;
1329   (*fact)->info.fill_ratio_given  = f;
1330   (*fact)->info.fill_ratio_needed = ((PetscReal)bi[n])/((PetscReal)ai[n]);
1331 
1332   ierr = MatILUFactorSymbolic_Inode(A,isrow,iscol,info,fact);CHKERRQ(ierr);
1333   (*fact)->ops->lufactornumeric =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
1334 
1335   PetscFunctionReturn(0);
1336 }
1337 
1338 #include "src/mat/impls/sbaij/seq/sbaij.h"
1339 #undef __FUNCT__
1340 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqAIJ"
1341 PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ(Mat A,MatFactorInfo *info,Mat *B)
1342 {
1343   Mat            C = *B;
1344   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data;
1345   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
1346   IS             ip=b->row,iip = b->icol;
1347   PetscErrorCode ierr;
1348   PetscInt       *rip,*riip,i,j,mbs=A->rmap.n,*bi=b->i,*bj=b->j,*bcol;
1349   PetscInt       *ai=a->i,*aj=a->j;
1350   PetscInt       k,jmin,jmax,*jl,*il,col,nexti,ili,nz;
1351   MatScalar      *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
1352   PetscReal      zeropivot,rs,shiftnz;
1353   PetscReal      shiftpd;
1354   ChShift_Ctx    sctx;
1355   PetscInt       newshift;
1356 
1357   PetscFunctionBegin;
1358 
1359   shiftnz   = info->shiftnz;
1360   shiftpd   = info->shiftpd;
1361   zeropivot = info->zeropivot;
1362 
1363   ierr  = ISGetIndices(ip,&rip);CHKERRQ(ierr);
1364   ierr  = ISGetIndices(iip,&riip);CHKERRQ(ierr);
1365 
1366   /* initialization */
1367   nz   = (2*mbs+1)*sizeof(PetscInt)+mbs*sizeof(MatScalar);
1368   ierr = PetscMalloc(nz,&il);CHKERRQ(ierr);
1369   jl   = il + mbs;
1370   rtmp = (MatScalar*)(jl + mbs);
1371 
1372   sctx.shift_amount = 0;
1373   sctx.nshift       = 0;
1374   do {
1375     sctx.chshift = PETSC_FALSE;
1376     for (i=0; i<mbs; i++) {
1377       rtmp[i] = 0.0; jl[i] = mbs; il[0] = 0;
1378     }
1379 
1380     for (k = 0; k<mbs; k++){
1381       bval = ba + bi[k];
1382       /* initialize k-th row by the perm[k]-th row of A */
1383       jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
1384       for (j = jmin; j < jmax; j++){
1385         col = riip[aj[j]];
1386         if (col >= k){ /* only take upper triangular entry */
1387           rtmp[col] = aa[j];
1388           *bval++  = 0.0; /* for in-place factorization */
1389         }
1390       }
1391       /* shift the diagonal of the matrix */
1392       if (sctx.nshift) rtmp[k] += sctx.shift_amount;
1393 
1394       /* modify k-th row by adding in those rows i with U(i,k)!=0 */
1395       dk = rtmp[k];
1396       i = jl[k]; /* first row to be added to k_th row  */
1397 
1398       while (i < k){
1399         nexti = jl[i]; /* next row to be added to k_th row */
1400 
1401         /* compute multiplier, update diag(k) and U(i,k) */
1402         ili = il[i];  /* index of first nonzero element in U(i,k:bms-1) */
1403         uikdi = - ba[ili]*ba[bi[i]];  /* diagonal(k) */
1404         dk += uikdi*ba[ili];
1405         ba[ili] = uikdi; /* -U(i,k) */
1406 
1407         /* add multiple of row i to k-th row */
1408         jmin = ili + 1; jmax = bi[i+1];
1409         if (jmin < jmax){
1410           for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
1411           /* update il and jl for row i */
1412           il[i] = jmin;
1413           j = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
1414         }
1415         i = nexti;
1416       }
1417 
1418       /* shift the diagonals when zero pivot is detected */
1419       /* compute rs=sum of abs(off-diagonal) */
1420       rs   = 0.0;
1421       jmin = bi[k]+1;
1422       nz   = bi[k+1] - jmin;
1423       bcol = bj + jmin;
1424       while (nz--){
1425         rs += PetscAbsScalar(rtmp[*bcol]);
1426         bcol++;
1427       }
1428 
1429       sctx.rs = rs;
1430       sctx.pv = dk;
1431       ierr = MatCholeskyCheckShift_inline(info,sctx,k,newshift);CHKERRQ(ierr);
1432 
1433       if (newshift == 1) {
1434         if (!sctx.shift_amount) {
1435           sctx.shift_amount = 1e-5;
1436         }
1437         break;
1438       }
1439 
1440       /* copy data into U(k,:) */
1441       ba[bi[k]] = 1.0/dk; /* U(k,k) */
1442       jmin = bi[k]+1; jmax = bi[k+1];
1443       if (jmin < jmax) {
1444         for (j=jmin; j<jmax; j++){
1445           col = bj[j]; ba[j] = rtmp[col]; rtmp[col] = 0.0;
1446         }
1447         /* add the k-th row into il and jl */
1448         il[k] = jmin;
1449         i = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
1450       }
1451     }
1452   } while (sctx.chshift);
1453   ierr = PetscFree(il);CHKERRQ(ierr);
1454 
1455   ierr = ISRestoreIndices(ip,&rip);CHKERRQ(ierr);
1456   ierr = ISRestoreIndices(iip,&riip);CHKERRQ(ierr);
1457   C->factor       = MAT_FACTOR_CHOLESKY;
1458   C->assembled    = PETSC_TRUE;
1459   C->preallocated = PETSC_TRUE;
1460   ierr = PetscLogFlops(C->rmap.n);CHKERRQ(ierr);
1461   if (sctx.nshift){
1462     if (shiftnz) {
1463       ierr = PetscInfo2(A,"number of shiftnz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);CHKERRQ(ierr);
1464     } else if (shiftpd) {
1465       ierr = PetscInfo2(A,"number of shiftpd tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);CHKERRQ(ierr);
1466     }
1467   }
1468   PetscFunctionReturn(0);
1469 }
1470 
1471 #undef __FUNCT__
1472 #define __FUNCT__ "MatICCFactorSymbolic_SeqAIJ"
1473 PetscErrorCode MatICCFactorSymbolic_SeqAIJ(Mat A,IS perm,MatFactorInfo *info,Mat *fact)
1474 {
1475   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data;
1476   Mat_SeqSBAIJ       *b;
1477   PetscErrorCode     ierr;
1478   PetscTruth         perm_identity,missing;
1479   PetscInt           reallocs=0,*rip,*riip,i,*ai=a->i,*aj=a->j,am=A->rmap.n,*ui;
1480   PetscInt           jmin,jmax,nzk,k,j,*jl,prow,*il,nextprow;
1481   PetscInt           nlnk,*lnk,*lnk_lvl=PETSC_NULL,d;
1482   PetscInt           ncols,ncols_upper,*cols,*ajtmp,*uj,**uj_ptr,**uj_lvl_ptr;
1483   PetscReal          fill=info->fill,levels=info->levels;
1484   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
1485   PetscFreeSpaceList free_space_lvl=PETSC_NULL,current_space_lvl=PETSC_NULL;
1486   PetscBT            lnkbt;
1487   IS                 iperm;
1488 
1489   PetscFunctionBegin;
1490   if (A->rmap.n != A->cmap.n) SETERRQ2(PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %D columns %D",A->rmap.n,A->cmap.n);
1491   ierr = MatMissingDiagonal(A,&missing,&d);CHKERRQ(ierr);
1492   if (missing) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
1493   ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr);
1494   ierr = ISInvertPermutation(perm,PETSC_DECIDE,&iperm);CHKERRQ(ierr);
1495 
1496   ierr = PetscMalloc((am+1)*sizeof(PetscInt),&ui);CHKERRQ(ierr);
1497   ui[0] = 0;
1498 
1499   /* ICC(0) without matrix ordering: simply copies fill pattern */
1500   if (!levels && perm_identity) {
1501 
1502     for (i=0; i<am; i++) {
1503       ui[i+1] = ui[i] + ai[i+1] - a->diag[i];
1504     }
1505     ierr = PetscMalloc((ui[am]+1)*sizeof(PetscInt),&uj);CHKERRQ(ierr);
1506     cols = uj;
1507     for (i=0; i<am; i++) {
1508       aj    = a->j + a->diag[i];
1509       ncols = ui[i+1] - ui[i];
1510       for (j=0; j<ncols; j++) *cols++ = *aj++;
1511     }
1512   } else { /* case: levels>0 || (levels=0 && !perm_identity) */
1513     ierr = ISGetIndices(iperm,&riip);CHKERRQ(ierr);
1514     ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr);
1515 
1516     /* initialization */
1517     ierr  = PetscMalloc((am+1)*sizeof(PetscInt),&ajtmp);CHKERRQ(ierr);
1518 
1519     /* jl: linked list for storing indices of the pivot rows
1520        il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
1521     ierr = PetscMalloc((2*am+1)*sizeof(PetscInt)+2*am*sizeof(PetscInt**),&jl);CHKERRQ(ierr);
1522     il         = jl + am;
1523     uj_ptr     = (PetscInt**)(il + am);
1524     uj_lvl_ptr = (PetscInt**)(uj_ptr + am);
1525     for (i=0; i<am; i++){
1526       jl[i] = am; il[i] = 0;
1527     }
1528 
1529     /* create and initialize a linked list for storing column indices of the active row k */
1530     nlnk = am + 1;
1531     ierr = PetscIncompleteLLCreate(am,am,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr);
1532 
1533     /* initial FreeSpace size is fill*(ai[am]+1) */
1534     ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+1)),&free_space);CHKERRQ(ierr);
1535     current_space = free_space;
1536     ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+1)),&free_space_lvl);CHKERRQ(ierr);
1537     current_space_lvl = free_space_lvl;
1538 
1539     for (k=0; k<am; k++){  /* for each active row k */
1540       /* initialize lnk by the column indices of row rip[k] of A */
1541       nzk   = 0;
1542       ncols = ai[rip[k]+1] - ai[rip[k]];
1543       if (!ncols) SETERRQ(PETSC_ERR_MAT_CH_ZRPVT,"Empty row in matrix");
1544       ncols_upper = 0;
1545       for (j=0; j<ncols; j++){
1546         i = *(aj + ai[rip[k]] + j); /* unpermuted column index */
1547         if (riip[i] >= k){ /* only take upper triangular entry */
1548           ajtmp[ncols_upper] = i;
1549           ncols_upper++;
1550         }
1551       }
1552       ierr = PetscIncompleteLLInit(ncols_upper,ajtmp,am,riip,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr);
1553       nzk += nlnk;
1554 
1555       /* update lnk by computing fill-in for each pivot row to be merged in */
1556       prow = jl[k]; /* 1st pivot row */
1557 
1558       while (prow < k){
1559         nextprow = jl[prow];
1560 
1561         /* merge prow into k-th row */
1562         jmin = il[prow] + 1;  /* index of the 2nd nzero entry in U(prow,k:am-1) */
1563         jmax = ui[prow+1];
1564         ncols = jmax-jmin;
1565         i     = jmin - ui[prow];
1566         cols  = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
1567         uj    = uj_lvl_ptr[prow] + i; /* levels of cols */
1568         j     = *(uj - 1);
1569         ierr = PetscICCLLAddSorted(ncols,cols,levels,uj,am,nlnk,lnk,lnk_lvl,lnkbt,j);CHKERRQ(ierr);
1570         nzk += nlnk;
1571 
1572         /* update il and jl for prow */
1573         if (jmin < jmax){
1574           il[prow] = jmin;
1575           j = *cols; jl[prow] = jl[j]; jl[j] = prow;
1576         }
1577         prow = nextprow;
1578       }
1579 
1580       /* if free space is not available, make more free space */
1581       if (current_space->local_remaining<nzk) {
1582         i = am - k + 1; /* num of unfactored rows */
1583         i = PetscMin(i*nzk, i*(i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1584         ierr = PetscFreeSpaceGet(i,&current_space);CHKERRQ(ierr);
1585         ierr = PetscFreeSpaceGet(i,&current_space_lvl);CHKERRQ(ierr);
1586         reallocs++;
1587       }
1588 
1589       /* copy data into free_space and free_space_lvl, then initialize lnk */
1590       if (nzk == 0) SETERRQ1(PETSC_ERR_ARG_WRONG,"Empty row %D in ICC matrix factor",k);
1591       ierr = PetscIncompleteLLClean(am,am,nzk,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);CHKERRQ(ierr);
1592 
1593       /* add the k-th row into il and jl */
1594       if (nzk > 1){
1595         i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
1596         jl[k] = jl[i]; jl[i] = k;
1597         il[k] = ui[k] + 1;
1598       }
1599       uj_ptr[k]     = current_space->array;
1600       uj_lvl_ptr[k] = current_space_lvl->array;
1601 
1602       current_space->array           += nzk;
1603       current_space->local_used      += nzk;
1604       current_space->local_remaining -= nzk;
1605 
1606       current_space_lvl->array           += nzk;
1607       current_space_lvl->local_used      += nzk;
1608       current_space_lvl->local_remaining -= nzk;
1609 
1610       ui[k+1] = ui[k] + nzk;
1611     }
1612 
1613 #if defined(PETSC_USE_INFO)
1614     if (ai[am] != 0) {
1615       PetscReal af = (PetscReal)ui[am]/((PetscReal)ai[am]);
1616       ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,fill,af);CHKERRQ(ierr);
1617       ierr = PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);CHKERRQ(ierr);
1618       ierr = PetscInfo1(A,"PCFactorSetFill(pc,%G) for best performance.\n",af);CHKERRQ(ierr);
1619     } else {
1620       ierr = PetscInfo(A,"Empty matrix.\n");CHKERRQ(ierr);
1621     }
1622 #endif
1623 
1624     ierr = ISRestoreIndices(perm,&rip);CHKERRQ(ierr);
1625     ierr = ISRestoreIndices(iperm,&riip);CHKERRQ(ierr);
1626     ierr = PetscFree(jl);CHKERRQ(ierr);
1627     ierr = PetscFree(ajtmp);CHKERRQ(ierr);
1628 
1629     /* destroy list of free space and other temporary array(s) */
1630     ierr = PetscMalloc((ui[am]+1)*sizeof(PetscInt),&uj);CHKERRQ(ierr);
1631     ierr = PetscFreeSpaceContiguous(&free_space,uj);CHKERRQ(ierr);
1632     ierr = PetscIncompleteLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
1633     ierr = PetscFreeSpaceDestroy(free_space_lvl);CHKERRQ(ierr);
1634 
1635   } /* end of case: levels>0 || (levels=0 && !perm_identity) */
1636 
1637   /* put together the new matrix in MATSEQSBAIJ format */
1638 
1639   b    = (Mat_SeqSBAIJ*)(*fact)->data;
1640   b->singlemalloc = PETSC_FALSE;
1641   ierr = PetscMalloc((ui[am]+1)*sizeof(MatScalar),&b->a);CHKERRQ(ierr);
1642   b->j    = uj;
1643   b->i    = ui;
1644   b->diag = 0;
1645   b->ilen = 0;
1646   b->imax = 0;
1647   b->row  = perm;
1648   b->col  = perm;
1649   ierr    = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr);
1650   ierr    = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr);
1651   b->icol = iperm;
1652   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
1653   ierr    = PetscMalloc((am+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
1654   ierr = PetscLogObjectMemory((*fact),(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr);
1655   b->maxnz   = b->nz = ui[am];
1656   b->free_a  = PETSC_TRUE;
1657   b->free_ij = PETSC_TRUE;
1658 
1659   (*fact)->factor                 = MAT_FACTOR_CHOLESKY;
1660   (*fact)->info.factor_mallocs    = reallocs;
1661   (*fact)->info.fill_ratio_given  = fill;
1662   if (ai[am] != 0) {
1663     (*fact)->info.fill_ratio_needed = ((PetscReal)ui[am])/((PetscReal)ai[am]);
1664   } else {
1665     (*fact)->info.fill_ratio_needed = 0.0;
1666   }
1667   (*fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ;
1668   if (perm_identity){
1669     (*fact)->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1670     (*fact)->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1671   }
1672   PetscFunctionReturn(0);
1673 }
1674 
1675 #undef __FUNCT__
1676 #define __FUNCT__ "MatCholeskyFactorSymbolic_SeqAIJ"
1677 PetscErrorCode MatCholeskyFactorSymbolic_SeqAIJ(Mat A,IS perm,MatFactorInfo *info,Mat *fact)
1678 {
1679   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data;
1680   Mat_SeqSBAIJ       *b;
1681   PetscErrorCode     ierr;
1682   PetscTruth         perm_identity;
1683   PetscReal          fill = info->fill;
1684   PetscInt           *rip,*riip,i,am=A->rmap.n,*ai=a->i,*aj=a->j,reallocs=0,prow;
1685   PetscInt           *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
1686   PetscInt           nlnk,*lnk,ncols,ncols_upper,*cols,*uj,**ui_ptr,*uj_ptr;
1687   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
1688   PetscBT            lnkbt;
1689   IS                 iperm;
1690 
1691   PetscFunctionBegin;
1692   if (A->rmap.n != A->cmap.n) SETERRQ2(PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %D columns %D",A->rmap.n,A->cmap.n);
1693   /* check whether perm is the identity mapping */
1694   ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr);
1695   ierr = ISInvertPermutation(perm,PETSC_DECIDE,&iperm);CHKERRQ(ierr);
1696   ierr = ISGetIndices(iperm,&riip);CHKERRQ(ierr);
1697   ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr);
1698 
1699   /* initialization */
1700   ierr  = PetscMalloc((am+1)*sizeof(PetscInt),&ui);CHKERRQ(ierr);
1701   ui[0] = 0;
1702 
1703   /* jl: linked list for storing indices of the pivot rows
1704      il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
1705   ierr = PetscMalloc((3*am+1)*sizeof(PetscInt)+am*sizeof(PetscInt**),&jl);CHKERRQ(ierr);
1706   il     = jl + am;
1707   cols   = il + am;
1708   ui_ptr = (PetscInt**)(cols + am);
1709   for (i=0; i<am; i++){
1710     jl[i] = am; il[i] = 0;
1711   }
1712 
1713   /* create and initialize a linked list for storing column indices of the active row k */
1714   nlnk = am + 1;
1715   ierr = PetscLLCreate(am,am,nlnk,lnk,lnkbt);CHKERRQ(ierr);
1716 
1717   /* initial FreeSpace size is fill*(ai[am]+1) */
1718   ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+1)),&free_space);CHKERRQ(ierr);
1719   current_space = free_space;
1720 
1721   for (k=0; k<am; k++){  /* for each active row k */
1722     /* initialize lnk by the column indices of row rip[k] of A */
1723     nzk   = 0;
1724     ncols = ai[rip[k]+1] - ai[rip[k]];
1725     if (!ncols) SETERRQ(PETSC_ERR_MAT_CH_ZRPVT,"Empty row in matrix");
1726     ncols_upper = 0;
1727     for (j=0; j<ncols; j++){
1728       i = riip[*(aj + ai[rip[k]] + j)];
1729       if (i >= k){ /* only take upper triangular entry */
1730         cols[ncols_upper] = i;
1731         ncols_upper++;
1732       }
1733     }
1734     ierr = PetscLLAdd(ncols_upper,cols,am,nlnk,lnk,lnkbt);CHKERRQ(ierr);
1735     nzk += nlnk;
1736 
1737     /* update lnk by computing fill-in for each pivot row to be merged in */
1738     prow = jl[k]; /* 1st pivot row */
1739 
1740     while (prow < k){
1741       nextprow = jl[prow];
1742       /* merge prow into k-th row */
1743       jmin = il[prow] + 1;  /* index of the 2nd nzero entry in U(prow,k:am-1) */
1744       jmax = ui[prow+1];
1745       ncols = jmax-jmin;
1746       uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:am-1) */
1747       ierr = PetscLLAddSorted(ncols,uj_ptr,am,nlnk,lnk,lnkbt);CHKERRQ(ierr);
1748       nzk += nlnk;
1749 
1750       /* update il and jl for prow */
1751       if (jmin < jmax){
1752         il[prow] = jmin;
1753         j = *uj_ptr; jl[prow] = jl[j]; jl[j] = prow;
1754       }
1755       prow = nextprow;
1756     }
1757 
1758     /* if free space is not available, make more free space */
1759     if (current_space->local_remaining<nzk) {
1760       i = am - k + 1; /* num of unfactored rows */
1761       i = PetscMin(i*nzk, i*(i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1762       ierr = PetscFreeSpaceGet(i,&current_space);CHKERRQ(ierr);
1763       reallocs++;
1764     }
1765 
1766     /* copy data into free space, then initialize lnk */
1767     ierr = PetscLLClean(am,am,nzk,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
1768 
1769     /* add the k-th row into il and jl */
1770     if (nzk-1 > 0){
1771       i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
1772       jl[k] = jl[i]; jl[i] = k;
1773       il[k] = ui[k] + 1;
1774     }
1775     ui_ptr[k] = current_space->array;
1776     current_space->array           += nzk;
1777     current_space->local_used      += nzk;
1778     current_space->local_remaining -= nzk;
1779 
1780     ui[k+1] = ui[k] + nzk;
1781   }
1782 
1783 #if defined(PETSC_USE_INFO)
1784   if (ai[am] != 0) {
1785     PetscReal af = (PetscReal)(ui[am])/((PetscReal)ai[am]);
1786     ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,fill,af);CHKERRQ(ierr);
1787     ierr = PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);CHKERRQ(ierr);
1788     ierr = PetscInfo1(A,"PCFactorSetFill(pc,%G) for best performance.\n",af);CHKERRQ(ierr);
1789   } else {
1790      ierr = PetscInfo(A,"Empty matrix.\n");CHKERRQ(ierr);
1791   }
1792 #endif
1793 
1794   ierr = ISRestoreIndices(perm,&rip);CHKERRQ(ierr);
1795   ierr = ISRestoreIndices(iperm,&riip);CHKERRQ(ierr);
1796   ierr = PetscFree(jl);CHKERRQ(ierr);
1797 
1798   /* destroy list of free space and other temporary array(s) */
1799   ierr = PetscMalloc((ui[am]+1)*sizeof(PetscInt),&uj);CHKERRQ(ierr);
1800   ierr = PetscFreeSpaceContiguous(&free_space,uj);CHKERRQ(ierr);
1801   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
1802 
1803   /* put together the new matrix in MATSEQSBAIJ format */
1804 
1805   b = (Mat_SeqSBAIJ*)(*fact)->data;
1806   b->singlemalloc = PETSC_FALSE;
1807   b->free_a       = PETSC_TRUE;
1808   b->free_ij      = PETSC_TRUE;
1809   ierr = PetscMalloc((ui[am]+1)*sizeof(MatScalar),&b->a);CHKERRQ(ierr);
1810   b->j    = uj;
1811   b->i    = ui;
1812   b->diag = 0;
1813   b->ilen = 0;
1814   b->imax = 0;
1815   b->row  = perm;
1816   b->col  = perm;
1817   ierr    = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr);
1818   ierr    = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr);
1819   b->icol = iperm;
1820   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
1821   ierr    = PetscMalloc((am+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
1822   ierr    = PetscLogObjectMemory(*fact,(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr);
1823   b->maxnz = b->nz = ui[am];
1824 
1825   (*fact)->factor                 = MAT_FACTOR_CHOLESKY;
1826   (*fact)->info.factor_mallocs    = reallocs;
1827   (*fact)->info.fill_ratio_given  = fill;
1828   if (ai[am] != 0) {
1829     (*fact)->info.fill_ratio_needed = ((PetscReal)ui[am])/((PetscReal)ai[am]);
1830   } else {
1831     (*fact)->info.fill_ratio_needed = 0.0;
1832   }
1833   (*fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ;
1834   if (perm_identity){
1835     (*fact)->ops->solve           = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1836     (*fact)->ops->solvetranspose  = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1837     (*fact)->ops->forwardsolve    = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
1838     (*fact)->ops->backwardsolve   = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;
1839   } else {
1840     (*fact)->ops->solve           = MatSolve_SeqSBAIJ_1;
1841     (*fact)->ops->solvetranspose  = MatSolve_SeqSBAIJ_1;
1842     (*fact)->ops->forwardsolve    = MatForwardSolve_SeqSBAIJ_1;
1843     (*fact)->ops->backwardsolve   = MatBackwardSolve_SeqSBAIJ_1;
1844   }
1845   PetscFunctionReturn(0);
1846 }
1847