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