xref: /petsc/src/mat/impls/aij/seq/aijfact.c (revision 94b7f48cc472a54ea2ce57edf1fe19e8a254237c) 
1 /*$Id: aijfact.c,v 1.167 2001/09/11 16:32:26 bsmith Exp $*/
2 
3 #include "src/mat/impls/aij/seq/aij.h"
4 #include "src/vec/vecimpl.h"
5 #include "src/inline/dot.h"
6 #include "src/inline/spops.h"
7 
8 #undef __FUNCT__
9 #define __FUNCT__ "MatOrdering_Flow_SeqAIJ"
10 int MatOrdering_Flow_SeqAIJ(Mat mat,MatOrderingType type,IS *irow,IS *icol)
11 {
12   PetscFunctionBegin;
13 
14   SETERRQ(PETSC_ERR_SUP,"Code not written");
15 #if !defined(PETSC_USE_DEBUG)
16   PetscFunctionReturn(0);
17 #endif
18 }
19 
20 
21 EXTERN int MatMarkDiagonal_SeqAIJ(Mat);
22 EXTERN int Mat_AIJ_CheckInode(Mat,PetscTruth);
23 
24 EXTERN int SPARSEKIT2dperm(int*,PetscScalar*,int*,int*,PetscScalar*,int*,int*,int*,int*,int*);
25 EXTERN int SPARSEKIT2ilutp(int*,PetscScalar*,int*,int*,int*,PetscReal,PetscReal*,int*,PetscScalar*,int*,int*,int*,PetscScalar*,int*,int*,int*);
26 EXTERN int SPARSEKIT2msrcsr(int*,PetscScalar*,int*,PetscScalar*,int*,int*,PetscScalar*,int*);
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 int MatILUDTFactor_SeqAIJ(Mat A,MatFactorInfo *info,IS isrow,IS iscol,Mat *fact)
52 {
53 #if defined(PETSC_AVOID_GNUCOPYRIGHT_CODE)
54   PetscFunctionBegin;
55   SETERRQ(1,"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   int          *c,*r,*ic,ierr,i,n = A->m;
63   int          *old_i = a->i,*old_j = a->j,*new_i,*old_i2 = 0,*old_j2 = 0,*new_j;
64   int          *ordcol,*iwk,*iperm,*jw;
65   int          jmax,lfill,job,*o_i,*o_j;
66   PetscScalar  *old_a = a->a,*w,*new_a,*old_a2 = 0,*wk,*o_a;
67   PetscReal    permtol,af;
68 
69   PetscFunctionBegin;
70 
71   if (info->dt == PETSC_DEFAULT)      info->dt      = .005;
72   if (info->dtcount == PETSC_DEFAULT) info->dtcount = (int)(1.5*a->rmax);
73   if (info->dtcol == PETSC_DEFAULT)   info->dtcol   = .01;
74   if (info->fill == PETSC_DEFAULT)    info->fill    = ((double)(n*(info->dtcount+1)))/a->nz;
75   lfill   = (int)(info->dtcount/2.0);
76   jmax    = (int)(info->fill*a->nz);
77   permtol = info->dtcol;
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(int),&new_i);CHKERRQ(ierr);
97   ierr = PetscMalloc(jmax*sizeof(int),&new_j);CHKERRQ(ierr);
98   ierr = PetscMalloc(jmax*sizeof(PetscScalar),&new_a);CHKERRQ(ierr);
99   ierr = PetscMalloc(n*sizeof(int),&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(int),&old_i2);CHKERRQ(ierr);
121     ierr = PetscMalloc((old_i[n]-old_i[0]+1)*sizeof(int),&old_j2);CHKERRQ(ierr);
122     ierr = PetscMalloc((old_i[n]-old_i[0]+1)*sizeof(PetscScalar),&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(int),&iperm);CHKERRQ(ierr);
145   ierr = PetscMalloc(2*n*sizeof(int),&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,&permtol,&n,new_a,new_j,new_i,&jmax,w,jw,iperm,&ierr);
149   if (ierr) {
150     switch (ierr) {
151       case -3: SETERRQ2(1,"ilutp(), matrix U overflows, need larger info->fill current fill %g space allocated %d",info->fill,jmax);
152       case -2: SETERRQ2(1,"ilutp(), matrix L overflows, need larger info->fill current fill %g space allocated %d",info->fill,jmax);
153       case -5: SETERRQ(1,"ilutp(), zero row encountered");
154       case -1: SETERRQ(1,"ilutp(), input matrix may be wrong");
155       case -4: SETERRQ1(1,"ilutp(), illegal info->fill value %d",jmax);
156       default: SETERRQ1(1,"ilutp(), zero pivot detected on row %d",ierr);
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(int),&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 = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,fact);CHKERRQ(ierr);
223   (*fact)->factor    = FACTOR_LU;
224   (*fact)->assembled = PETSC_TRUE;
225 
226   b = (Mat_SeqAIJ*)(*fact)->data;
227   ierr = PetscFree(b->imax);CHKERRQ(ierr);
228   b->sorted        = PETSC_FALSE;
229   b->singlemalloc  = PETSC_FALSE;
230   /* the next line frees the default space generated by the MatCreate() */
231   ierr             = PetscFree(b->a);CHKERRQ(ierr);
232   ierr             = PetscFree(b->ilen);CHKERRQ(ierr);
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   ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);
247 
248   /* check out for identical nodes. If found, use inode functions */
249   ierr = Mat_AIJ_CheckInode(*fact,PETSC_FALSE);CHKERRQ(ierr);
250 
251   af = ((double)b->nz)/((double)a->nz) + .001;
252   PetscLogInfo(A,"MatILUDTFactor_SeqAIJ:Fill ratio:given %g needed %g\n",info->fill,af);
253   PetscLogInfo(A,"MatILUDTFactor_SeqAIJ:Run with -pc_ilu_fill %g or use \n",af);
254   PetscLogInfo(A,"MatILUDTFactor_SeqAIJ:PCILUSetFill(pc,%g);\n",af);
255   PetscLogInfo(A,"MatILUDTFactor_SeqAIJ:for best performance.\n");
256 
257   PetscFunctionReturn(0);
258 #endif
259 }
260 
261 /*
262     Factorization code for AIJ format.
263 */
264 #undef __FUNCT__
265 #define __FUNCT__ "MatLUFactorSymbolic_SeqAIJ"
266 int MatLUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatFactorInfo *info,Mat *B)
267 {
268   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b;
269   IS         isicol;
270   int        *r,*ic,ierr,i,n = A->m,*ai = a->i,*aj = a->j;
271   int        *ainew,*ajnew,jmax,*fill,*ajtmp,nz;
272   int        *idnew,idx,row,m,fm,nnz,nzi,realloc = 0,nzbd,*im;
273   PetscReal  f;
274 
275   PetscFunctionBegin;
276   if (A->M != A->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
277   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
278   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
279   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
280 
281   /* get new row pointers */
282   ierr = PetscMalloc((n+1)*sizeof(int),&ainew);CHKERRQ(ierr);
283   ainew[0] = 0;
284   /* don't know how many column pointers are needed so estimate */
285   f = info->fill;
286   jmax  = (int)(f*ai[n]+1);
287   ierr = PetscMalloc((jmax)*sizeof(int),&ajnew);CHKERRQ(ierr);
288   /* fill is a linked list of nonzeros in active row */
289   ierr = PetscMalloc((2*n+1)*sizeof(int),&fill);CHKERRQ(ierr);
290   im   = fill + n + 1;
291   /* idnew is location of diagonal in factor */
292   ierr = PetscMalloc((n+1)*sizeof(int),&idnew);CHKERRQ(ierr);
293   idnew[0] = 0;
294 
295   for (i=0; i<n; i++) {
296     /* first copy previous fill into linked list */
297     nnz     = nz    = ai[r[i]+1] - ai[r[i]];
298     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
299     ajtmp   = aj + ai[r[i]];
300     fill[n] = n;
301     while (nz--) {
302       fm  = n;
303       idx = ic[*ajtmp++];
304       do {
305         m  = fm;
306         fm = fill[m];
307       } while (fm < idx);
308       fill[m]   = idx;
309       fill[idx] = fm;
310     }
311     row = fill[n];
312     while (row < i) {
313       ajtmp = ajnew + idnew[row] + 1;
314       nzbd  = 1 + idnew[row] - ainew[row];
315       nz    = im[row] - nzbd;
316       fm    = row;
317       while (nz-- > 0) {
318         idx = *ajtmp++ ;
319         nzbd++;
320         if (idx == i) im[row] = nzbd;
321         do {
322           m  = fm;
323           fm = fill[m];
324         } while (fm < idx);
325         if (fm != idx) {
326           fill[m]   = idx;
327           fill[idx] = fm;
328           fm        = idx;
329           nnz++;
330         }
331       }
332       row = fill[row];
333     }
334     /* copy new filled row into permanent storage */
335     ainew[i+1] = ainew[i] + nnz;
336     if (ainew[i+1] > jmax) {
337 
338       /* estimate how much additional space we will need */
339       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
340       /* just double the memory each time */
341       int maxadd = jmax;
342       /* maxadd = (int)((f*(ai[n]+(!shift))*(n-i+5))/n); */
343       if (maxadd < nnz) maxadd = (n-i)*(nnz+1);
344       jmax += maxadd;
345 
346       /* allocate a longer ajnew */
347       ierr = PetscMalloc(jmax*sizeof(int),&ajtmp);CHKERRQ(ierr);
348       ierr  = PetscMemcpy(ajtmp,ajnew,(ainew[i])*sizeof(int));CHKERRQ(ierr);
349       ierr  = PetscFree(ajnew);CHKERRQ(ierr);
350       ajnew = ajtmp;
351       realloc++; /* count how many times we realloc */
352     }
353     ajtmp = ajnew + ainew[i];
354     fm    = fill[n];
355     nzi   = 0;
356     im[i] = nnz;
357     while (nnz--) {
358       if (fm < i) nzi++;
359       *ajtmp++ = fm ;
360       fm       = fill[fm];
361     }
362     idnew[i] = ainew[i] + nzi;
363   }
364   if (ai[n] != 0) {
365     PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
366     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af);
367     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Run with -pc_lu_fill %g or use \n",af);
368     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:PCLUSetFill(pc,%g);\n",af);
369     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:for best performance.\n");
370   } else {
371     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ: Empty matrix\n");
372   }
373 
374   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
375   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
376 
377   ierr = PetscFree(fill);CHKERRQ(ierr);
378 
379   /* put together the new matrix */
380   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,B);CHKERRQ(ierr);
381   PetscLogObjectParent(*B,isicol);
382   b = (Mat_SeqAIJ*)(*B)->data;
383   ierr = PetscFree(b->imax);CHKERRQ(ierr);
384   b->singlemalloc = PETSC_FALSE;
385   /* the next line frees the default space generated by the Create() */
386   ierr = PetscFree(b->a);CHKERRQ(ierr);
387   ierr = PetscFree(b->ilen);CHKERRQ(ierr);
388   ierr          = PetscMalloc((ainew[n]+1)*sizeof(PetscScalar),&b->a);CHKERRQ(ierr);
389   b->j          = ajnew;
390   b->i          = ainew;
391   b->diag       = idnew;
392   b->ilen       = 0;
393   b->imax       = 0;
394   b->row        = isrow;
395   b->col        = iscol;
396   b->lu_damping        = info->damping;
397   b->lu_zeropivot      = info->zeropivot;
398   b->lu_shift          = info->shift;
399   b->lu_shift_fraction = info->shift_fraction;
400   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
401   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
402   b->icol       = isicol;
403   ierr          = PetscMalloc((n+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
404   /* In b structure:  Free imax, ilen, old a, old j.
405      Allocate idnew, solve_work, new a, new j */
406   PetscLogObjectMemory(*B,(ainew[n]-n)*(sizeof(int)+sizeof(PetscScalar)));
407   b->maxnz = b->nz = ainew[n] ;
408 
409   (*B)->factor                 =  FACTOR_LU;
410   (*B)->info.factor_mallocs    = realloc;
411   (*B)->info.fill_ratio_given  = f;
412   ierr = Mat_AIJ_CheckInode(*B,PETSC_FALSE);CHKERRQ(ierr);
413   (*B)->ops->lufactornumeric   =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
414 
415   if (ai[n] != 0) {
416     (*B)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
417   } else {
418     (*B)->info.fill_ratio_needed = 0.0;
419   }
420   PetscFunctionReturn(0);
421 }
422 /* ----------------------------------------------------------- */
423 EXTERN int Mat_AIJ_CheckInode(Mat,PetscTruth);
424 
425 #undef __FUNCT__
426 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJ"
427 int MatLUFactorNumeric_SeqAIJ(Mat A,Mat *B)
428 {
429   Mat          C = *B;
430   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ *)C->data;
431   IS           isrow = b->row,isicol = b->icol;
432   int          *r,*ic,ierr,i,j,n = A->m,*ai = b->i,*aj = b->j;
433   int          *ajtmpold,*ajtmp,nz,row,*ics;
434   int          *diag_offset = b->diag,diag,*pj,ndamp = 0, nshift=0;
435   PetscScalar  *rtmp,*v,*pc,multiplier,*pv,*rtmps;
436   PetscReal    damping = b->lu_damping, zeropivot = b->lu_zeropivot,rs,d;
437   PetscReal    row_shift,shift_fraction,shift_amount,shift_lo=0., shift_hi=1., shift_top=0.;
438   PetscTruth   damp,lushift;
439 
440   PetscFunctionBegin;
441   ierr  = ISGetIndices(isrow,&r);CHKERRQ(ierr);
442   ierr  = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
443   ierr  = PetscMalloc((n+1)*sizeof(PetscScalar),&rtmp);CHKERRQ(ierr);
444   ierr  = PetscMemzero(rtmp,(n+1)*sizeof(PetscScalar));CHKERRQ(ierr);
445   rtmps = rtmp; ics = ic;
446 
447   if (!a->diag) {
448     ierr = MatMarkDiagonal_SeqAIJ(A); CHKERRQ(ierr);
449   }
450 
451   if (b->lu_shift) { /* set max shift */
452     int *aai = a->i,*ddiag = a->diag;
453     shift_top = 0;
454     for (i=0; i<n; i++) {
455       d = PetscAbsScalar((a->a)[ddiag[i]]);
456       /* calculate amt of shift needed for this row */
457       if (d<=0) {
458         row_shift = 0;
459       } else {
460         row_shift = -2*d;
461       }
462       v = a->a+aai[i];
463       for (j=0; j<aai[i+1]-aai[i]; j++)
464 	row_shift += PetscAbsScalar(v[j]);
465       if (row_shift>shift_top) shift_top = row_shift;
466     }
467   }
468 
469   shift_fraction = 0; shift_amount = 0;
470   do {
471     damp    = PETSC_FALSE;
472     lushift = PETSC_FALSE;
473     for (i=0; i<n; i++) {
474       nz    = ai[i+1] - ai[i];
475       ajtmp = aj + ai[i];
476       for  (j=0; j<nz; j++) rtmps[ajtmp[j]] = 0.0;
477 
478       /* load in initial (unfactored row) */
479       nz       = a->i[r[i]+1] - a->i[r[i]];
480       ajtmpold = a->j + a->i[r[i]];
481       v        = a->a + a->i[r[i]];
482       for (j=0; j<nz; j++) {
483         rtmp[ics[ajtmpold[j]]] = v[j];
484       }
485       rtmp[ics[r[i]]] += damping + shift_amount; /* damp the diagonal of the matrix */
486 
487       row = *ajtmp++ ;
488       while  (row < i) {
489         pc = rtmp + row;
490         if (*pc != 0.0) {
491           pv         = b->a + diag_offset[row] ;
492           pj         = b->j + diag_offset[row] + 1;
493           multiplier = *pc / *pv++;
494           *pc        = multiplier;
495           nz         = ai[row+1] - diag_offset[row] - 1;
496           for (j=0; j<nz; j++) rtmps[pj[j]] -= multiplier * pv[j];
497           PetscLogFlops(2*nz);
498         }
499         row = *ajtmp++;
500       }
501       /* finished row so stick it into b->a */
502       pv   = b->a + ai[i] ;
503       pj   = b->j + ai[i] ;
504       nz   = ai[i+1] - ai[i];
505       diag = diag_offset[i] - ai[i];
506       /* 9/13/02 Victor Eijkhout suggested scaling zeropivot by rs for matrices with funny scalings */
507       rs   = 0.0;
508       for (j=0; j<nz; j++) {
509         pv[j] = rtmps[pj[j]];
510         if (j != diag) rs += PetscAbsScalar(pv[j]);
511       }
512 #define MAX_NSHIFT 5
513       if (PetscRealPart(pv[diag]) < zeropivot*rs && b->lu_shift) {
514 	if (nshift>MAX_NSHIFT) {
515 	  SETERRQ(1,"Unable to determine shift to enforce positive definite preconditioner");
516 	} else if (nshift==MAX_NSHIFT) {
517 	  shift_fraction = shift_hi;
518 	  lushift      = PETSC_FALSE;
519 	} else {
520 	  shift_lo = shift_fraction; shift_fraction = (shift_hi+shift_lo)/2.;
521 	  lushift      = PETSC_TRUE;
522 	}
523 	shift_amount = shift_fraction * shift_top;
524         nshift++;
525         break;
526       }
527       if (PetscAbsScalar(pv[diag]) < zeropivot*rs) {
528         if (damping) {
529           if (ndamp) damping *= 2.0;
530           damp = PETSC_TRUE;
531           ndamp++;
532           break;
533         } else {
534           SETERRQ4(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %d value %g tolerance %g * rs %g",i,PetscAbsScalar(pv[diag]),zeropivot,rs);
535         }
536       }
537     }
538     if (!lushift && b->lu_shift && shift_fraction>0 && nshift<MAX_NSHIFT) {
539       /*
540        * if no shift in this attempt & shifting & started shifting & can refine,
541        * then try lower shift
542        */
543       shift_hi       = shift_fraction;
544       shift_fraction = (shift_hi+shift_lo)/2.;
545       shift_amount   = shift_fraction * shift_top;
546       lushift        = PETSC_TRUE;
547       nshift++;
548     }
549   } while (damp || lushift);
550 
551   /* invert diagonal entries for simplier triangular solves */
552   for (i=0; i<n; i++) {
553     b->a[diag_offset[i]] = 1.0/b->a[diag_offset[i]];
554   }
555 
556   ierr = PetscFree(rtmp);CHKERRQ(ierr);
557   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
558   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
559   C->factor = FACTOR_LU;
560   (*B)->ops->lufactornumeric   =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
561   C->assembled = PETSC_TRUE;
562   PetscLogFlops(C->n);
563   if (ndamp) {
564     PetscLogInfo(0,"MatLUFactorNumerical_SeqAIJ: number of damping tries %d damping value %g\n",ndamp,damping);
565   }
566   if (nshift) {
567     b->lu_shift_fraction = shift_fraction;
568     PetscLogInfo(0,"MatLUFactorNumerical_SeqAIJ: diagonal shifted up by %e fraction top_value %e number shifts %d\n",shift_fraction,shift_top,nshift);
569   }
570   PetscFunctionReturn(0);
571 }
572 
573 #undef __FUNCT__
574 #define __FUNCT__ "MatUsePETSc_SeqAIJ"
575 int MatUsePETSc_SeqAIJ(Mat A)
576 {
577   PetscFunctionBegin;
578   A->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ;
579   A->ops->lufactornumeric  = MatLUFactorNumeric_SeqAIJ;
580   PetscFunctionReturn(0);
581 }
582 
583 
584 /* ----------------------------------------------------------- */
585 #undef __FUNCT__
586 #define __FUNCT__ "MatLUFactor_SeqAIJ"
587 int MatLUFactor_SeqAIJ(Mat A,IS row,IS col,MatFactorInfo *info)
588 {
589   int ierr;
590   Mat C;
591 
592   PetscFunctionBegin;
593   ierr = MatLUFactorSymbolic(A,row,col,info,&C);CHKERRQ(ierr);
594   ierr = MatLUFactorNumeric(A,&C);CHKERRQ(ierr);
595   ierr = MatHeaderCopy(A,C);CHKERRQ(ierr);
596   PetscLogObjectParent(A,((Mat_SeqAIJ*)(A->data))->icol);
597   PetscFunctionReturn(0);
598 }
599 /* ----------------------------------------------------------- */
600 #undef __FUNCT__
601 #define __FUNCT__ "MatSolve_SeqAIJ"
602 int MatSolve_SeqAIJ(Mat A,Vec bb,Vec xx)
603 {
604   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
605   IS           iscol = a->col,isrow = a->row;
606   int          *r,*c,ierr,i, n = A->m,*vi,*ai = a->i,*aj = a->j;
607   int          nz,*rout,*cout;
608   PetscScalar  *x,*b,*tmp,*tmps,*aa = a->a,sum,*v;
609 
610   PetscFunctionBegin;
611   if (!n) PetscFunctionReturn(0);
612 
613   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
614   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
615   tmp  = a->solve_work;
616 
617   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
618   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
619 
620   /* forward solve the lower triangular */
621   tmp[0] = b[*r++];
622   tmps   = tmp;
623   for (i=1; i<n; i++) {
624     v   = aa + ai[i] ;
625     vi  = aj + ai[i] ;
626     nz  = a->diag[i] - ai[i];
627     sum = b[*r++];
628     SPARSEDENSEMDOT(sum,tmps,v,vi,nz);
629     tmp[i] = sum;
630   }
631 
632   /* backward solve the upper triangular */
633   for (i=n-1; i>=0; i--){
634     v   = aa + a->diag[i] + 1;
635     vi  = aj + a->diag[i] + 1;
636     nz  = ai[i+1] - a->diag[i] - 1;
637     sum = tmp[i];
638     SPARSEDENSEMDOT(sum,tmps,v,vi,nz);
639     x[*c--] = tmp[i] = sum*aa[a->diag[i]];
640   }
641 
642   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
643   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
644   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
645   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
646   PetscLogFlops(2*a->nz - A->n);
647   PetscFunctionReturn(0);
648 }
649 
650 /* ----------------------------------------------------------- */
651 #undef __FUNCT__
652 #define __FUNCT__ "MatSolve_SeqAIJ_NaturalOrdering"
653 int MatSolve_SeqAIJ_NaturalOrdering(Mat A,Vec bb,Vec xx)
654 {
655   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
656   int          n = A->m,*ai = a->i,*aj = a->j,*adiag = a->diag,ierr;
657   PetscScalar  *x,*b,*aa = a->a;
658 #if !defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
659   int          adiag_i,i,*vi,nz,ai_i;
660   PetscScalar  *v,sum;
661 #endif
662 
663   PetscFunctionBegin;
664   if (!n) PetscFunctionReturn(0);
665 
666   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
667   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
668 
669 #if defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
670   fortransolveaij_(&n,x,ai,aj,adiag,aa,b);
671 #else
672   /* forward solve the lower triangular */
673   x[0] = b[0];
674   for (i=1; i<n; i++) {
675     ai_i = ai[i];
676     v    = aa + ai_i;
677     vi   = aj + ai_i;
678     nz   = adiag[i] - ai_i;
679     sum  = b[i];
680     while (nz--) sum -= *v++ * x[*vi++];
681     x[i] = sum;
682   }
683 
684   /* backward solve the upper triangular */
685   for (i=n-1; i>=0; i--){
686     adiag_i = adiag[i];
687     v       = aa + adiag_i + 1;
688     vi      = aj + adiag_i + 1;
689     nz      = ai[i+1] - adiag_i - 1;
690     sum     = x[i];
691     while (nz--) sum -= *v++ * x[*vi++];
692     x[i]    = sum*aa[adiag_i];
693   }
694 #endif
695   PetscLogFlops(2*a->nz - A->n);
696   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
697   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
698   PetscFunctionReturn(0);
699 }
700 
701 #undef __FUNCT__
702 #define __FUNCT__ "MatSolveAdd_SeqAIJ"
703 int MatSolveAdd_SeqAIJ(Mat A,Vec bb,Vec yy,Vec xx)
704 {
705   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
706   IS           iscol = a->col,isrow = a->row;
707   int          *r,*c,ierr,i, n = A->m,*vi,*ai = a->i,*aj = a->j;
708   int          nz,*rout,*cout;
709   PetscScalar  *x,*b,*tmp,*aa = a->a,sum,*v;
710 
711   PetscFunctionBegin;
712   if (yy != xx) {ierr = VecCopy(yy,xx);CHKERRQ(ierr);}
713 
714   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
715   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
716   tmp  = a->solve_work;
717 
718   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
719   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
720 
721   /* forward solve the lower triangular */
722   tmp[0] = b[*r++];
723   for (i=1; i<n; i++) {
724     v   = aa + ai[i] ;
725     vi  = aj + ai[i] ;
726     nz  = a->diag[i] - ai[i];
727     sum = b[*r++];
728     while (nz--) sum -= *v++ * tmp[*vi++ ];
729     tmp[i] = sum;
730   }
731 
732   /* backward solve the upper triangular */
733   for (i=n-1; i>=0; i--){
734     v   = aa + a->diag[i] + 1;
735     vi  = aj + a->diag[i] + 1;
736     nz  = ai[i+1] - a->diag[i] - 1;
737     sum = tmp[i];
738     while (nz--) sum -= *v++ * tmp[*vi++ ];
739     tmp[i] = sum*aa[a->diag[i]];
740     x[*c--] += tmp[i];
741   }
742 
743   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
744   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
745   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
746   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
747   PetscLogFlops(2*a->nz);
748 
749   PetscFunctionReturn(0);
750 }
751 /* -------------------------------------------------------------------*/
752 #undef __FUNCT__
753 #define __FUNCT__ "MatSolveTranspose_SeqAIJ"
754 int MatSolveTranspose_SeqAIJ(Mat A,Vec bb,Vec xx)
755 {
756   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
757   IS           iscol = a->col,isrow = a->row;
758   int          *r,*c,ierr,i,n = A->m,*vi,*ai = a->i,*aj = a->j;
759   int          nz,*rout,*cout,*diag = a->diag;
760   PetscScalar  *x,*b,*tmp,*aa = a->a,*v,s1;
761 
762   PetscFunctionBegin;
763   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
764   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
765   tmp  = a->solve_work;
766 
767   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
768   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
769 
770   /* copy the b into temp work space according to permutation */
771   for (i=0; i<n; i++) tmp[i] = b[c[i]];
772 
773   /* forward solve the U^T */
774   for (i=0; i<n; i++) {
775     v   = aa + diag[i] ;
776     vi  = aj + diag[i] + 1;
777     nz  = ai[i+1] - diag[i] - 1;
778     s1  = tmp[i];
779     s1 *= (*v++);  /* multiply by inverse of diagonal entry */
780     while (nz--) {
781       tmp[*vi++ ] -= (*v++)*s1;
782     }
783     tmp[i] = s1;
784   }
785 
786   /* backward solve the L^T */
787   for (i=n-1; i>=0; i--){
788     v   = aa + diag[i] - 1 ;
789     vi  = aj + diag[i] - 1 ;
790     nz  = diag[i] - ai[i];
791     s1  = tmp[i];
792     while (nz--) {
793       tmp[*vi-- ] -= (*v--)*s1;
794     }
795   }
796 
797   /* copy tmp into x according to permutation */
798   for (i=0; i<n; i++) x[r[i]] = tmp[i];
799 
800   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
801   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
802   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
803   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
804 
805   PetscLogFlops(2*a->nz-A->n);
806   PetscFunctionReturn(0);
807 }
808 
809 #undef __FUNCT__
810 #define __FUNCT__ "MatSolveTransposeAdd_SeqAIJ"
811 int MatSolveTransposeAdd_SeqAIJ(Mat A,Vec bb,Vec zz,Vec xx)
812 {
813   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
814   IS           iscol = a->col,isrow = a->row;
815   int          *r,*c,ierr,i,n = A->m,*vi,*ai = a->i,*aj = a->j;
816   int          nz,*rout,*cout,*diag = a->diag;
817   PetscScalar  *x,*b,*tmp,*aa = a->a,*v;
818 
819   PetscFunctionBegin;
820   if (zz != xx) VecCopy(zz,xx);
821 
822   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
823   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
824   tmp = a->solve_work;
825 
826   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
827   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
828 
829   /* copy the b into temp work space according to permutation */
830   for (i=0; i<n; i++) tmp[i] = b[c[i]];
831 
832   /* forward solve the U^T */
833   for (i=0; i<n; i++) {
834     v   = aa + diag[i] ;
835     vi  = aj + diag[i] + 1;
836     nz  = ai[i+1] - diag[i] - 1;
837     tmp[i] *= *v++;
838     while (nz--) {
839       tmp[*vi++ ] -= (*v++)*tmp[i];
840     }
841   }
842 
843   /* backward solve the L^T */
844   for (i=n-1; i>=0; i--){
845     v   = aa + diag[i] - 1 ;
846     vi  = aj + diag[i] - 1 ;
847     nz  = diag[i] - ai[i];
848     while (nz--) {
849       tmp[*vi-- ] -= (*v--)*tmp[i];
850     }
851   }
852 
853   /* copy tmp into x according to permutation */
854   for (i=0; i<n; i++) x[r[i]] += tmp[i];
855 
856   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
857   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
858   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
859   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
860 
861   PetscLogFlops(2*a->nz);
862   PetscFunctionReturn(0);
863 }
864 /* ----------------------------------------------------------------*/
865 EXTERN int MatMissingDiagonal_SeqAIJ(Mat);
866 
867 #undef __FUNCT__
868 #define __FUNCT__ "MatILUFactorSymbolic_SeqAIJ"
869 int MatILUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatFactorInfo *info,Mat *fact)
870 {
871   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b;
872   IS         isicol;
873   int        *r,*ic,ierr,prow,n = A->m,*ai = a->i,*aj = a->j;
874   int        *ainew,*ajnew,jmax,*fill,*xi,nz,*im,*ajfill,*flev;
875   int        *dloc,idx,row,m,fm,nzf,nzi,len, realloc = 0,dcount = 0;
876   int        incrlev,nnz,i,levels,diagonal_fill;
877   PetscTruth col_identity,row_identity;
878   PetscReal  f;
879 
880   PetscFunctionBegin;
881   f             = info->fill;
882   levels        = (int)info->levels;
883   diagonal_fill = (int)info->diagonal_fill;
884   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
885 
886   /* special case that simply copies fill pattern */
887   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
888   ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
889   if (!levels && row_identity && col_identity) {
890     ierr = MatDuplicate_SeqAIJ(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr);
891     (*fact)->factor = FACTOR_LU;
892     b               = (Mat_SeqAIJ*)(*fact)->data;
893     if (!b->diag) {
894       ierr = MatMarkDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
895     }
896     ierr = MatMissingDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
897     b->row              = isrow;
898     b->col              = iscol;
899     b->icol             = isicol;
900     b->lu_damping       = info->damping;
901     b->lu_zeropivot     = info->zeropivot;
902     b->lu_shift         = info->shift;
903     b->lu_shift_fraction= info->shift_fraction;
904     ierr                = PetscMalloc(((*fact)->m+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
905     (*fact)->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
906     ierr                = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
907     ierr                = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
908     PetscFunctionReturn(0);
909   }
910 
911   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
912   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
913 
914   /* get new row pointers */
915   ierr = PetscMalloc((n+1)*sizeof(int),&ainew);CHKERRQ(ierr);
916   ainew[0] = 0;
917   /* don't know how many column pointers are needed so estimate */
918   jmax = (int)(f*(ai[n]+1));
919   ierr = PetscMalloc((jmax)*sizeof(int),&ajnew);CHKERRQ(ierr);
920   /* ajfill is level of fill for each fill entry */
921   ierr = PetscMalloc((jmax)*sizeof(int),&ajfill);CHKERRQ(ierr);
922   /* fill is a linked list of nonzeros in active row */
923   ierr = PetscMalloc((n+1)*sizeof(int),&fill);CHKERRQ(ierr);
924   /* im is level for each filled value */
925   ierr = PetscMalloc((n+1)*sizeof(int),&im);CHKERRQ(ierr);
926   /* dloc is location of diagonal in factor */
927   ierr = PetscMalloc((n+1)*sizeof(int),&dloc);CHKERRQ(ierr);
928   dloc[0]  = 0;
929   for (prow=0; prow<n; prow++) {
930 
931     /* copy current row into linked list */
932     nzf     = nz  = ai[r[prow]+1] - ai[r[prow]];
933     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
934     xi      = aj + ai[r[prow]] ;
935     fill[n]    = n;
936     fill[prow] = -1; /* marker to indicate if diagonal exists */
937     while (nz--) {
938       fm  = n;
939       idx = ic[*xi++ ];
940       do {
941         m  = fm;
942         fm = fill[m];
943       } while (fm < idx);
944       fill[m]   = idx;
945       fill[idx] = fm;
946       im[idx]   = 0;
947     }
948 
949     /* make sure diagonal entry is included */
950     if (diagonal_fill && fill[prow] == -1) {
951       fm = n;
952       while (fill[fm] < prow) fm = fill[fm];
953       fill[prow] = fill[fm]; /* insert diagonal into linked list */
954       fill[fm]   = prow;
955       im[prow]   = 0;
956       nzf++;
957       dcount++;
958     }
959 
960     nzi = 0;
961     row = fill[n];
962     while (row < prow) {
963       incrlev = im[row] + 1;
964       nz      = dloc[row];
965       xi      = ajnew  + ainew[row]  + nz + 1;
966       flev    = ajfill + ainew[row]  + nz + 1;
967       nnz     = ainew[row+1] - ainew[row] - nz - 1;
968       fm      = row;
969       while (nnz-- > 0) {
970         idx = *xi++ ;
971         if (*flev + incrlev > levels) {
972           flev++;
973           continue;
974         }
975         do {
976           m  = fm;
977           fm = fill[m];
978         } while (fm < idx);
979         if (fm != idx) {
980           im[idx]   = *flev + incrlev;
981           fill[m]   = idx;
982           fill[idx] = fm;
983           fm        = idx;
984           nzf++;
985         } else {
986           if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev;
987         }
988         flev++;
989       }
990       row = fill[row];
991       nzi++;
992     }
993     /* copy new filled row into permanent storage */
994     ainew[prow+1] = ainew[prow] + nzf;
995     if (ainew[prow+1] > jmax) {
996 
997       /* estimate how much additional space we will need */
998       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
999       /* just double the memory each time */
1000       /*  maxadd = (int)((f*(ai[n]+!shift)*(n-prow+5))/n); */
1001       int maxadd = jmax;
1002       if (maxadd < nzf) maxadd = (n-prow)*(nzf+1);
1003       jmax += maxadd;
1004 
1005       /* allocate a longer ajnew and ajfill */
1006       ierr   = PetscMalloc(jmax*sizeof(int),&xi);CHKERRQ(ierr);
1007       ierr   = PetscMemcpy(xi,ajnew,(ainew[prow])*sizeof(int));CHKERRQ(ierr);
1008       ierr   = PetscFree(ajnew);CHKERRQ(ierr);
1009       ajnew  = xi;
1010       ierr   = PetscMalloc(jmax*sizeof(int),&xi);CHKERRQ(ierr);
1011       ierr   = PetscMemcpy(xi,ajfill,(ainew[prow])*sizeof(int));CHKERRQ(ierr);
1012       ierr   = PetscFree(ajfill);CHKERRQ(ierr);
1013       ajfill = xi;
1014       realloc++; /* count how many times we realloc */
1015     }
1016     xi          = ajnew + ainew[prow] ;
1017     flev        = ajfill + ainew[prow] ;
1018     dloc[prow]  = nzi;
1019     fm          = fill[n];
1020     while (nzf--) {
1021       *xi++   = fm ;
1022       *flev++ = im[fm];
1023       fm      = fill[fm];
1024     }
1025     /* make sure row has diagonal entry */
1026     if (ajnew[ainew[prow]+dloc[prow]] != prow) {
1027       SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,"Row %d has missing diagonal in factored matrix\n\
1028     try running with -pc_ilu_nonzeros_along_diagonal or -pc_ilu_diagonal_fill",prow);
1029     }
1030   }
1031   ierr = PetscFree(ajfill);CHKERRQ(ierr);
1032   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
1033   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
1034   ierr = PetscFree(fill);CHKERRQ(ierr);
1035   ierr = PetscFree(im);CHKERRQ(ierr);
1036 
1037   {
1038     PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
1039     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af);
1040     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Run with -[sub_]pc_ilu_fill %g or use \n",af);
1041     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:PCILUSetFill([sub]pc,%g);\n",af);
1042     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:for best performance.\n");
1043     if (diagonal_fill) {
1044       PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Detected and replaced %d missing diagonals",dcount);
1045     }
1046   }
1047 
1048   /* put together the new matrix */
1049   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,fact);CHKERRQ(ierr);
1050   PetscLogObjectParent(*fact,isicol);
1051   b = (Mat_SeqAIJ*)(*fact)->data;
1052   ierr = PetscFree(b->imax);CHKERRQ(ierr);
1053   b->singlemalloc = PETSC_FALSE;
1054   len = (ainew[n] )*sizeof(PetscScalar);
1055   /* the next line frees the default space generated by the Create() */
1056   ierr = PetscFree(b->a);CHKERRQ(ierr);
1057   ierr = PetscFree(b->ilen);CHKERRQ(ierr);
1058   ierr = PetscMalloc(len+1,&b->a);CHKERRQ(ierr);
1059   b->j          = ajnew;
1060   b->i          = ainew;
1061   for (i=0; i<n; i++) dloc[i] += ainew[i];
1062   b->diag       = dloc;
1063   b->ilen       = 0;
1064   b->imax       = 0;
1065   b->row        = isrow;
1066   b->col        = iscol;
1067   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
1068   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
1069   b->icol       = isicol;
1070   ierr = PetscMalloc((n+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
1071   /* In b structure:  Free imax, ilen, old a, old j.
1072      Allocate dloc, solve_work, new a, new j */
1073   PetscLogObjectMemory(*fact,(ainew[n]-n) * (sizeof(int)+sizeof(PetscScalar)));
1074   b->maxnz             = b->nz = ainew[n] ;
1075   b->lu_damping        = info->damping;
1076   b->lu_shift          = info->shift;
1077   b->lu_shift_fraction = info->shift_fraction;
1078   b->lu_zeropivot = info->zeropivot;
1079   (*fact)->factor = FACTOR_LU;
1080   ierr = Mat_AIJ_CheckInode(*fact,PETSC_FALSE);CHKERRQ(ierr);
1081   (*fact)->ops->lufactornumeric =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
1082 
1083   (*fact)->info.factor_mallocs    = realloc;
1084   (*fact)->info.fill_ratio_given  = f;
1085   (*fact)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[prow]);
1086   PetscFunctionReturn(0);
1087 }
1088 
1089 #include "src/mat/impls/sbaij/seq/sbaij.h"
1090 #undef __FUNCT__
1091 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqAIJ"
1092 int MatCholeskyFactorNumeric_SeqAIJ(Mat A,Mat *fact)
1093 {
1094   Mat_SeqAIJ          *a = (Mat_SeqAIJ*)A->data;
1095   int                 ierr;
1096 
1097   PetscFunctionBegin;
1098   if (!a->sbaijMat){
1099     ierr = MatConvert(A,MATSEQSBAIJ,&a->sbaijMat);CHKERRQ(ierr);
1100   }
1101 
1102   ierr = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering(a->sbaijMat,fact);CHKERRQ(ierr);
1103   ierr = MatDestroy(a->sbaijMat);CHKERRQ(ierr);
1104   a->sbaijMat = PETSC_NULL;
1105 
1106   PetscFunctionReturn(0);
1107 }
1108 
1109 #undef __FUNCT__
1110 #define __FUNCT__ "MatICCFactorSymbolic_SeqAIJ"
1111 int MatICCFactorSymbolic_SeqAIJ(Mat A,IS perm,MatFactorInfo *info,Mat *fact)
1112 {
1113   Mat_SeqAIJ          *a = (Mat_SeqAIJ*)A->data;
1114   int                 ierr;
1115   PetscTruth          perm_identity;
1116 
1117   PetscFunctionBegin;
1118   ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr);
1119   if (!perm_identity){
1120     SETERRQ(1,"Non-identity permutation is not supported yet");
1121   }
1122   if (!a->sbaijMat){
1123     ierr = MatConvert(A,MATSEQSBAIJ,&a->sbaijMat);CHKERRQ(ierr);
1124   }
1125 
1126   ierr = MatICCFactorSymbolic(a->sbaijMat,perm,info,fact);CHKERRQ(ierr);
1127   (*fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ;
1128 
1129   PetscFunctionReturn(0);
1130 }
1131 
1132 #undef __FUNCT__
1133 #define __FUNCT__ "MatCholeskyFactorSymbolic_SeqAIJ"
1134 int MatCholeskyFactorSymbolic_SeqAIJ(Mat A,IS perm,MatFactorInfo *info,Mat *fact)
1135 {
1136   Mat_SeqAIJ          *a = (Mat_SeqAIJ*)A->data;
1137   int                 ierr;
1138   PetscTruth          perm_identity;
1139 
1140   PetscFunctionBegin;
1141   ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr);
1142   if (!perm_identity){
1143     SETERRQ(1,"Non-identity permutation is not supported yet");
1144   }
1145   if (!a->sbaijMat){
1146     ierr = MatConvert(A,MATSEQSBAIJ,&a->sbaijMat);CHKERRQ(ierr);
1147   }
1148 
1149   ierr = MatCholeskyFactorSymbolic(a->sbaijMat,perm,info,fact);CHKERRQ(ierr);
1150   (*fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ;
1151 
1152   PetscFunctionReturn(0);
1153 }
1154