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