xref: /petsc/src/mat/impls/aij/seq/aijfact.c (revision b0a32e0c6855ee6a6cd3495fa7da12ea9885bc5d) 
1 /*$Id: aijfact.c,v 1.157 2000/10/24 20:25:32 bsmith Exp bsmith $*/
2 
3 #include "src/mat/impls/aij/seq/aij.h"
4 #include "src/vec/vecimpl.h"
5 #include "src/inline/dot.h"
6 
7 #undef __FUNC__
8 #define __FUNC__ "MatOrdering_Flow_SeqAIJ"
9 int MatOrdering_Flow_SeqAIJ(Mat mat,MatOrderingType type,IS *irow,IS *icol)
10 {
11   PetscFunctionBegin;
12 
13   SETERRQ(PETSC_ERR_SUP,"Code not written");
14 #if !defined(PETSC_USE_DEBUG)
15   PetscFunctionReturn(0);
16 #endif
17 }
18 
19 
20 EXTERN int MatMarkDiagonal_SeqAIJ(Mat);
21 EXTERN int Mat_AIJ_CheckInode(Mat);
22 
23 EXTERN int SPARSEKIT2dperm(int*,Scalar*,int*,int*,Scalar*,int*,int*,int*,int*,int*);
24 EXTERN int SPARSEKIT2ilutp(int*,Scalar*,int*,int*,int*,PetscReal*,PetscReal*,int*,Scalar*,int*,int*,int*,Scalar*,int*,int*,int*);
25 EXTERN int SPARSEKIT2msrcsr(int*,Scalar*,int*,Scalar*,int*,int*,Scalar*,int*);
26 
27 #undef __FUNC__
28 #define __FUNC__ "MatILUDTFactor_SeqAIJ"
29   /* ------------------------------------------------------------
30 
31           This interface was contribed by Tony Caola
32 
33      This routine is an interface to the pivoting drop-tolerance
34      ILU routine written by Yousef Saad (saad@cs.umn.edu) as part of
35      SPARSEKIT2.
36 
37      The SPARSEKIT2 routines used here are covered by the GNU
38      copyright; see the file gnu in this directory.
39 
40      Thanks to Prof. Saad, Dr. Hysom, and Dr. Smith for their
41      help in getting this routine ironed out.
42 
43      The major drawback to this routine is that if info->fill is
44      not large enough it fails rather than allocating more space;
45      this can be fixed by hacking/improving the f2c version of
46      Yousef Saad's code.
47 
48      ishift = 0, for indices start at 1
49      ishift = 1, for indices starting at 0
50      ------------------------------------------------------------
51   */
52 
53 int MatILUDTFactor_SeqAIJ(Mat A,MatILUInfo *info,IS isrow,IS iscol,Mat *fact)
54 {
55   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b;
56   IS         iscolf,isicol,isirow;
57   PetscTruth reorder;
58   int        *c,*r,*ic,ierr,i,n = A->m;
59   int        *old_i = a->i,*old_j = a->j,*new_i,*old_i2 = 0,*old_j2 = 0,*new_j;
60   int        *ordcol,*iwk,*iperm,*jw;
61   int        ishift = !a->indexshift;
62   int        jmax,lfill,job,*o_i,*o_j;
63   Scalar     *old_a = a->a,*w,*new_a,*old_a2 = 0,*wk,*o_a;
64   PetscReal  permtol,af;
65 
66   PetscFunctionBegin;
67 
68   if (info->dt == PETSC_DEFAULT)      info->dt      = .005;
69   if (info->dtcount == PETSC_DEFAULT) info->dtcount = (int)(1.5*a->rmax);
70   if (info->dtcol == PETSC_DEFAULT)   info->dtcol   = .01;
71   if (info->fill == PETSC_DEFAULT)    info->fill    = ((double)(n*(info->dtcount+1)))/a->nz;
72   lfill   = (int)(info->dtcount/2.0);
73   jmax    = (int)(info->fill*a->nz);
74   permtol = info->dtcol;
75 
76 
77   /* ------------------------------------------------------------
78      If reorder=.TRUE., then the original matrix has to be
79      reordered to reflect the user selected ordering scheme, and
80      then de-reordered so it is in it's original format.
81      Because Saad's dperm() is NOT in place, we have to copy
82      the original matrix and allocate more storage. . .
83      ------------------------------------------------------------
84   */
85 
86   /* set reorder to true if either isrow or iscol is not identity */
87   ierr = ISIdentity(isrow,&reorder);CHKERRQ(ierr);
88   if (reorder) {ierr = ISIdentity(iscol,&reorder);CHKERRQ(ierr);}
89   reorder = PetscNot(reorder);
90 
91 
92   /* storage for ilu factor */
93   ierr = PetscMalloc((n+1)*sizeof(int),&new_i);CHKERRQ(ierr);
94   ierr = PetscMalloc(jmax*sizeof(int),&new_j);CHKERRQ(ierr);
95   ierr = PetscMalloc(jmax*sizeof(Scalar),&new_a);CHKERRQ(ierr);
96   ierr = PetscMalloc(n*sizeof(int),&ordcol);CHKERRQ(ierr);
97 
98   /* ------------------------------------------------------------
99      Make sure that everything is Fortran formatted (1-Based)
100      ------------------------------------------------------------
101   */
102   if (ishift) {
103     for (i=old_i[0];i<old_i[n];i++) {
104       old_j[i]++;
105     }
106     for(i=0;i<n+1;i++) {
107       old_i[i]++;
108     };
109   };
110 
111   if (reorder) {
112     ierr = ISGetIndices(iscol,&c);           CHKERRQ(ierr);
113     ierr = ISGetIndices(isrow,&r);           CHKERRQ(ierr);
114     for(i=0;i<n;i++) {
115       r[i]  = r[i]+1;
116       c[i]  = c[i]+1;
117     }
118     ierr = PetscMalloc((n+1)*sizeof(int),&old_i2);CHKERRQ(ierr);
119     ierr = PetscMalloc((old_i[n]-old_i[0]+1)*sizeof(int),&old_j2);CHKERRQ(ierr);
120     ierr = PetscMalloc((old_i[n]-old_i[0]+1)*sizeof(Scalar),old_a2);CHKERRQ(ierr);
121     job  = 3; SPARSEKIT2dperm(&n,old_a,old_j,old_i,old_a2,old_j2,old_i2,r,c,&job);
122     for (i=0;i<n;i++) {
123       r[i]  = r[i]-1;
124       c[i]  = c[i]-1;
125     }
126     ierr = ISRestoreIndices(iscol,&c);CHKERRQ(ierr);
127     ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
128     o_a = old_a2;
129     o_j = old_j2;
130     o_i = old_i2;
131   } else {
132     o_a = old_a;
133     o_j = old_j;
134     o_i = old_i;
135   }
136 
137   /* ------------------------------------------------------------
138      Call Saad's ilutp() routine to generate the factorization
139      ------------------------------------------------------------
140   */
141 
142   ierr = PetscMalloc(2*n*sizeof(int),&iperm);CHKERRQ(ierr);
143   ierr = PetscMalloc(2*n*sizeof(int),&jw);CHKERRQ(ierr);
144   ierr = PetscMalloc(n*sizeof(Scalar),&w);CHKERRQ(ierr);
145 
146   SPARSEKIT2ilutp(&n,o_a,o_j,o_i,&lfill,&info->dt,&permtol,&n,new_a,new_j,new_i,&jmax,w,jw,iperm,&ierr);
147   if (ierr) {
148     switch (ierr) {
149       case -3: SETERRQ2(1,"ilutp(), matrix U overflows, need larger info->fill current fill %g space allocated %d",info->fill,jmax);
150       case -2: SETERRQ2(1,"ilutp(), matrix L overflows, need larger info->fill current fill %g space allocated %d",info->fill,jmax);
151       case -5: SETERRQ(1,"ilutp(), zero row encountered");
152       case -1: SETERRQ(1,"ilutp(), input matrix may be wrong");
153       case -4: SETERRQ1(1,"ilutp(), illegal info->fill value %d",jmax);
154       default: SETERRQ1(1,"ilutp(), zero pivot detected on row %d",ierr);
155     }
156   }
157 
158   ierr = PetscFree(w);CHKERRQ(ierr);
159   ierr = PetscFree(jw);CHKERRQ(ierr);
160 
161   /* ------------------------------------------------------------
162      Saad's routine gives the result in Modified Sparse Row (msr)
163      Convert to Compressed Sparse Row format (csr)
164      ------------------------------------------------------------
165   */
166 
167   ierr = PetscMalloc(n*sizeof(Scalar),&wk);CHKERRQ(ierr);
168   ierr = PetscMalloc((n+1)*sizeof(int),&iwk);CHKERRQ(ierr);
169 
170   SPARSEKIT2msrcsr(&n,new_a,new_j,new_a,new_j,new_i,wk,iwk);
171 
172   ierr = PetscFree(iwk);CHKERRQ(ierr);
173   ierr = PetscFree(wk);CHKERRQ(ierr);
174 
175   if (reorder) {
176     ierr = PetscFree(old_a2);CHKERRQ(ierr);
177     ierr = PetscFree(old_j2);CHKERRQ(ierr);
178     ierr = PetscFree(old_i2);CHKERRQ(ierr);
179   } else {
180     /* fix permutation of old_j that the factorization introduced */
181     for (i=old_i[0]; i<old_i[n]; i++) {
182       old_j[i-1] = iperm[old_j[i-1]-1];
183     }
184   }
185 
186   /* get rid of the shift to indices starting at 1 */
187   if (ishift) {
188     for (i=0; i<n+1; i++) {
189       old_i[i]--;
190     }
191     for (i=old_i[0];i<old_i[n];i++) {
192       old_j[i]--;
193     }
194   }
195 
196   /* Make the factored matrix 0-based */
197   if (ishift) {
198     for (i=0; i<n+1; i++) {
199       new_i[i]--;
200     }
201     for (i=new_i[0];i<new_i[n];i++) {
202       new_j[i]--;
203     }
204   }
205 
206   /*-- due to the pivoting, we need to reorder iscol to correctly --*/
207   /*-- permute the right-hand-side and solution vectors           --*/
208   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
209   ierr = ISInvertPermutation(isrow,PETSC_DECIDE,&isirow);CHKERRQ(ierr);
210   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
211   for(i=0; i<n; i++) {
212     ordcol[i] = ic[iperm[i]-1];
213   };
214   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
215   ierr = ISDestroy(isicol);CHKERRQ(ierr);
216 
217   ierr = PetscFree(iperm);CHKERRQ(ierr);
218 
219   ierr = ISCreateGeneral(PETSC_COMM_SELF,n,ordcol,&iscolf);CHKERRQ(ierr);
220   ierr = PetscFree(ordcol);CHKERRQ(ierr);
221 
222   /*----- put together the new matrix -----*/
223 
224   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,fact);CHKERRQ(ierr);
225   (*fact)->factor    = FACTOR_LU;
226   (*fact)->assembled = PETSC_TRUE;
227 
228   b = (Mat_SeqAIJ*)(*fact)->data;
229   ierr = PetscFree(b->imax);CHKERRQ(ierr);
230   b->sorted        = PETSC_FALSE;
231   b->singlemalloc  = PETSC_FALSE;
232   /* the next line frees the default space generated by the MatCreate() */
233   ierr             = PetscFree(b->a);CHKERRQ(ierr);
234   ierr             = PetscFree(b->ilen);CHKERRQ(ierr);
235   b->a             = new_a;
236   b->j             = new_j;
237   b->i             = new_i;
238   b->ilen          = 0;
239   b->imax          = 0;
240   /*  I am not sure why these are the inverses of the row and column permutations; but the other way is NO GOOD */
241   b->row           = isirow;
242   b->col           = iscolf;
243   ierr = PetscMalloc((n+1)*sizeof(Scalar),&b->solve_work);CHKERRQ(ierr);
244   b->maxnz = b->nz = new_i[n];
245   b->indexshift    = a->indexshift;
246   ierr = MatMarkDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
247   (*fact)->info.factor_mallocs = 0;
248 
249   ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);
250 
251   /* check out for identical nodes. If found, use inode functions */
252   ierr = Mat_AIJ_CheckInode(*fact);CHKERRQ(ierr);
253 
254   af = ((double)b->nz)/((double)a->nz) + .001;
255   PetscLogInfo(A,"MatILUDTFactor_SeqAIJ:Fill ratio:given %g needed %g\n",info->fill,af);
256   PetscLogInfo(A,"MatILUDTFactor_SeqAIJ:Run with -pc_ilu_fill %g or use \n",af);
257   PetscLogInfo(A,"MatILUDTFactor_SeqAIJ:PCILUSetFill(pc,%g);\n",af);
258   PetscLogInfo(A,"MatILUDTFactor_SeqAIJ:for best performance.\n");
259 
260   PetscFunctionReturn(0);
261 }
262 
263 /*
264     Factorization code for AIJ format.
265 */
266 #undef __FUNC__
267 #define __FUNC__ /*<a name="MatLUFactorSymbolic_SeqAIJ""></a>*/"MatLUFactorSymbolic_SeqAIJ"
268 int MatLUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatLUInfo *info,Mat *B)
269 {
270   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b;
271   IS         isicol;
272   int        *r,*ic,ierr,i,n = A->m,*ai = a->i,*aj = a->j;
273   int        *ainew,*ajnew,jmax,*fill,*ajtmp,nz,shift = a->indexshift;
274   int        *idnew,idx,row,m,fm,nnz,nzi,realloc = 0,nzbd,*im;
275   PetscReal  f;
276 
277   PetscFunctionBegin;
278   PetscValidHeaderSpecific(isrow,IS_COOKIE);
279   PetscValidHeaderSpecific(iscol,IS_COOKIE);
280   if (A->M != A->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
281 
282   if (!isrow) {
283     ierr = ISCreateStride(PETSC_COMM_SELF,A->M,0,1,&isrow);CHKERRQ(ierr);
284   }
285   if (!iscol) {
286     ierr = ISCreateStride(PETSC_COMM_SELF,A->M,0,1,&iscol);CHKERRQ(ierr);
287   }
288 
289   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
290   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
291   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
292 
293   /* get new row pointers */
294   ierr = PetscMalloc((n+1)*sizeof(int),&ainew);CHKERRQ(ierr);
295   ainew[0] = -shift;
296   /* don't know how many column pointers are needed so estimate */
297   if (info) f = info->fill; else f = 1.0;
298   jmax  = (int)(f*ai[n]+(!shift));
299   ierr = PetscMalloc((jmax)*sizeof(int),&ajnew);CHKERRQ(ierr);
300   /* fill is a linked list of nonzeros in active row */
301   ierr = PetscMalloc((2*n+1)*sizeof(int),&fill);CHKERRQ(ierr);
302   im   = fill + n + 1;
303   /* idnew is location of diagonal in factor */
304   ierr = PetscMalloc((n+1)*sizeof(int),&idnew);CHKERRQ(ierr);
305   idnew[0] = -shift;
306 
307   for (i=0; i<n; i++) {
308     /* first copy previous fill into linked list */
309     nnz     = nz    = ai[r[i]+1] - ai[r[i]];
310     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
311     ajtmp   = aj + ai[r[i]] + shift;
312     fill[n] = n;
313     while (nz--) {
314       fm  = n;
315       idx = ic[*ajtmp++ + shift];
316       do {
317         m  = fm;
318         fm = fill[m];
319       } while (fm < idx);
320       fill[m]   = idx;
321       fill[idx] = fm;
322     }
323     row = fill[n];
324     while (row < i) {
325       ajtmp = ajnew + idnew[row] + (!shift);
326       nzbd  = 1 + idnew[row] - ainew[row];
327       nz    = im[row] - nzbd;
328       fm    = row;
329       while (nz-- > 0) {
330         idx = *ajtmp++ + shift;
331         nzbd++;
332         if (idx == i) im[row] = nzbd;
333         do {
334           m  = fm;
335           fm = fill[m];
336         } while (fm < idx);
337         if (fm != idx) {
338           fill[m]   = idx;
339           fill[idx] = fm;
340           fm        = idx;
341           nnz++;
342         }
343       }
344       row = fill[row];
345     }
346     /* copy new filled row into permanent storage */
347     ainew[i+1] = ainew[i] + nnz;
348     if (ainew[i+1] > jmax) {
349 
350       /* estimate how much additional space we will need */
351       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
352       /* just double the memory each time */
353       int maxadd = jmax;
354       /* maxadd = (int)((f*(ai[n]+(!shift))*(n-i+5))/n); */
355       if (maxadd < nnz) maxadd = (n-i)*(nnz+1);
356       jmax += maxadd;
357 
358       /* allocate a longer ajnew */
359       ierr = PetscMalloc(jmax*sizeof(int),&ajtmp);CHKERRQ(ierr);
360       ierr  = PetscMemcpy(ajtmp,ajnew,(ainew[i]+shift)*sizeof(int));CHKERRQ(ierr);
361       ierr  = PetscFree(ajnew);CHKERRQ(ierr);
362       ajnew = ajtmp;
363       realloc++; /* count how many times we realloc */
364     }
365     ajtmp = ajnew + ainew[i] + shift;
366     fm    = fill[n];
367     nzi   = 0;
368     im[i] = nnz;
369     while (nnz--) {
370       if (fm < i) nzi++;
371       *ajtmp++ = fm - shift;
372       fm       = fill[fm];
373     }
374     idnew[i] = ainew[i] + nzi;
375   }
376   if (ai[n] != 0) {
377     PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
378     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af);
379     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Run with -pc_lu_fill %g or use \n",af);
380     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:PCLUSetFill(pc,%g);\n",af);
381     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:for best performance.\n");
382   } else {
383     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ: Empty matrix\n");
384   }
385 
386   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
387   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
388 
389   ierr = PetscFree(fill);CHKERRQ(ierr);
390 
391   /* put together the new matrix */
392   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,B);CHKERRQ(ierr);
393   PetscLogObjectParent(*B,isicol);
394   b = (Mat_SeqAIJ*)(*B)->data;
395   ierr = PetscFree(b->imax);CHKERRQ(ierr);
396   b->singlemalloc = PETSC_FALSE;
397   /* the next line frees the default space generated by the Create() */
398   ierr = PetscFree(b->a);CHKERRQ(ierr);
399   ierr = PetscFree(b->ilen);CHKERRQ(ierr);
400   ierr          = PetscMalloc((ainew[n]+shift+1)*sizeof(Scalar),&b->a);CHKERRQ(ierr);
401   b->j          = ajnew;
402   b->i          = ainew;
403   b->diag       = idnew;
404   b->ilen       = 0;
405   b->imax       = 0;
406   b->row        = isrow;
407   b->col        = iscol;
408   if (info) {
409     b->lu_damp    = (PetscTruth) ((int)info->damp);
410     b->lu_damping = info->damping;
411   } else {
412     b->lu_damp    = PETSC_FALSE;
413     b->lu_damping = 0.0;
414   }
415   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
416   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
417   b->icol       = isicol;
418   ierr          = PetscMalloc((n+1)*sizeof(Scalar),&b->solve_work);CHKERRQ(ierr);
419   /* In b structure:  Free imax, ilen, old a, old j.
420      Allocate idnew, solve_work, new a, new j */
421   PetscLogObjectMemory(*B,(ainew[n]+shift-n)*(sizeof(int)+sizeof(Scalar)));
422   b->maxnz = b->nz = ainew[n] + shift;
423 
424   (*B)->factor                 =  FACTOR_LU;
425   (*B)->info.factor_mallocs    = realloc;
426   (*B)->info.fill_ratio_given  = f;
427   (*B)->ops->lufactornumeric   =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
428 
429   if (ai[n] != 0) {
430     (*B)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
431   } else {
432     (*B)->info.fill_ratio_needed = 0.0;
433   }
434   PetscFunctionReturn(0);
435 }
436 /* ----------------------------------------------------------- */
437 EXTERN int Mat_AIJ_CheckInode(Mat);
438 
439 #undef __FUNC__
440 #define __FUNC__ "MatLUFactorNumeric_SeqAIJ"
441 int MatLUFactorNumeric_SeqAIJ(Mat A,Mat *B)
442 {
443   Mat        C = *B;
444   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ *)C->data;
445   IS         isrow = b->row,isicol = b->icol;
446   int        *r,*ic,ierr,i,j,n = A->m,*ai = b->i,*aj = b->j;
447   int        *ajtmpold,*ajtmp,nz,row,*ics,shift = a->indexshift;
448   int        *diag_offset = b->diag,diag;
449   int        *pj,ndamp = 0;
450   Scalar     *rtmp,*v,*pc,multiplier,*pv,*rtmps;
451   PetscReal  damping = b->lu_damping;
452   PetscTruth damp;
453 
454   PetscFunctionBegin;
455 
456   ierr  = ISGetIndices(isrow,&r);CHKERRQ(ierr);
457   ierr  = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
458   ierr = PetscMalloc((n+1)*sizeof(Scalar),&rtmp);CHKERRQ(ierr);
459   ierr  = PetscMemzero(rtmp,(n+1)*sizeof(Scalar));CHKERRQ(ierr);
460   rtmps = rtmp + shift; ics = ic + shift;
461 
462   do {
463     damp = PETSC_FALSE;
464     for (i=0; i<n; i++) {
465       nz    = ai[i+1] - ai[i];
466       ajtmp = aj + ai[i] + shift;
467       for  (j=0; j<nz; j++) rtmps[ajtmp[j]] = 0.0;
468 
469       /* load in initial (unfactored row) */
470       nz       = a->i[r[i]+1] - a->i[r[i]];
471       ajtmpold = a->j + a->i[r[i]] + shift;
472       v        = a->a + a->i[r[i]] + shift;
473       for (j=0; j<nz; j++) {
474         rtmp[ics[ajtmpold[j]]] = v[j];
475         if (ajtmpold[j] == r[i]) { /* damp the diagonal of the matrix */
476           rtmp[ics[ajtmpold[j]]] += damping;
477         }
478       }
479 
480       row = *ajtmp++ + shift;
481       while  (row < i) {
482         pc = rtmp + row;
483         if (*pc != 0.0) {
484           pv         = b->a + diag_offset[row] + shift;
485           pj         = b->j + diag_offset[row] + (!shift);
486           multiplier = *pc / *pv++;
487           *pc        = multiplier;
488           nz         = ai[row+1] - diag_offset[row] - 1;
489           for (j=0; j<nz; j++) rtmps[pj[j]] -= multiplier * pv[j];
490           PetscLogFlops(2*nz);
491         }
492         row = *ajtmp++ + shift;
493       }
494       /* finished row so stick it into b->a */
495       pv = b->a + ai[i] + shift;
496       pj = b->j + ai[i] + shift;
497       nz = ai[i+1] - ai[i];
498       for (j=0; j<nz; j++) {pv[j] = rtmps[pj[j]];}
499       diag = diag_offset[i] - ai[i];
500       if (PetscAbsScalar(pv[diag]) < 1.e-12) {
501         if (b->lu_damp) {
502           damp = PETSC_TRUE;
503           if (damping) damping *= 2.0;
504           else damping = 1.e-12;
505           ndamp++;
506           break;
507         } else {
508           SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %d",i);
509         }
510       }
511     }
512   } while (damp);
513 
514   /* invert diagonal entries for simplier triangular solves */
515   for (i=0; i<n; i++) {
516     b->a[diag_offset[i]+shift] = 1.0/b->a[diag_offset[i]+shift];
517   }
518 
519   ierr = PetscFree(rtmp);CHKERRQ(ierr);
520   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
521   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
522   C->factor = FACTOR_LU;
523   ierr = Mat_AIJ_CheckInode(C);CHKERRQ(ierr);
524   (*B)->ops->lufactornumeric   =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
525   C->assembled = PETSC_TRUE;
526   PetscLogFlops(C->n);
527   if (ndamp || b->lu_damping) {
528     PetscLogInfo(0,"MatLUFactorNumerical_SeqAIJ: number of damping tries %d damping value %g\n",ndamp,damping);
529   }
530   PetscFunctionReturn(0);
531 }
532 /* ----------------------------------------------------------- */
533 #undef __FUNC__
534 #define __FUNC__ "MatLUFactor_SeqAIJ"
535 int MatLUFactor_SeqAIJ(Mat A,IS row,IS col,MatLUInfo *info)
536 {
537   int ierr;
538   Mat C;
539 
540   PetscFunctionBegin;
541   ierr = MatLUFactorSymbolic(A,row,col,info,&C);CHKERRQ(ierr);
542   ierr = MatLUFactorNumeric(A,&C);CHKERRQ(ierr);
543   ierr = MatHeaderCopy(A,C);CHKERRQ(ierr);
544   PetscFunctionReturn(0);
545 }
546 /* ----------------------------------------------------------- */
547 #undef __FUNC__
548 #define __FUNC__ "MatSolve_SeqAIJ"
549 int MatSolve_SeqAIJ(Mat A,Vec bb,Vec xx)
550 {
551   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
552   IS         iscol = a->col,isrow = a->row;
553   int        *r,*c,ierr,i, n = A->m,*vi,*ai = a->i,*aj = a->j;
554   int        nz,shift = a->indexshift,*rout,*cout;
555   Scalar     *x,*b,*tmp,*tmps,*aa = a->a,sum,*v;
556 
557   PetscFunctionBegin;
558   if (!n) PetscFunctionReturn(0);
559 
560   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
561   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
562   tmp  = a->solve_work;
563 
564   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
565   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
566 
567   /* forward solve the lower triangular */
568   tmp[0] = b[*r++];
569   tmps   = tmp + shift;
570   for (i=1; i<n; i++) {
571     v   = aa + ai[i] + shift;
572     vi  = aj + ai[i] + shift;
573     nz  = a->diag[i] - ai[i];
574     sum = b[*r++];
575     while (nz--) sum -= *v++ * tmps[*vi++];
576     tmp[i] = sum;
577   }
578 
579   /* backward solve the upper triangular */
580   for (i=n-1; i>=0; i--){
581     v   = aa + a->diag[i] + (!shift);
582     vi  = aj + a->diag[i] + (!shift);
583     nz  = ai[i+1] - a->diag[i] - 1;
584     sum = tmp[i];
585     while (nz--) sum -= *v++ * tmps[*vi++];
586     x[*c--] = tmp[i] = sum*aa[a->diag[i]+shift];
587   }
588 
589   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
590   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
591   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
592   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
593   PetscLogFlops(2*a->nz - A->n);
594   PetscFunctionReturn(0);
595 }
596 
597 /* ----------------------------------------------------------- */
598 #undef __FUNC__
599 #define __FUNC__ "MatSolve_SeqAIJ_NaturalOrdering"
600 int MatSolve_SeqAIJ_NaturalOrdering(Mat A,Vec bb,Vec xx)
601 {
602   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
603   int        n = A->m,*ai = a->i,*aj = a->j,*adiag = a->diag,ierr;
604   Scalar     *x,*b,*aa = a->a,sum;
605 #if !defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
606   int        adiag_i,i,*vi,nz,ai_i;
607   Scalar     *v;
608 #endif
609 
610   PetscFunctionBegin;
611   if (!n) PetscFunctionReturn(0);
612   if (a->indexshift) {
613      ierr = MatSolve_SeqAIJ(A,bb,xx);CHKERRQ(ierr);
614      PetscFunctionReturn(0);
615   }
616 
617   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
618   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
619 
620 #if defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
621   fortransolveaij_(&n,x,ai,aj,adiag,aa,b);
622 #else
623   /* forward solve the lower triangular */
624   x[0] = b[0];
625   for (i=1; i<n; i++) {
626     ai_i = ai[i];
627     v    = aa + ai_i;
628     vi   = aj + ai_i;
629     nz   = adiag[i] - ai_i;
630     sum  = b[i];
631     while (nz--) sum -= *v++ * x[*vi++];
632     x[i] = sum;
633   }
634 
635   /* backward solve the upper triangular */
636   for (i=n-1; i>=0; i--){
637     adiag_i = adiag[i];
638     v       = aa + adiag_i + 1;
639     vi      = aj + adiag_i + 1;
640     nz      = ai[i+1] - adiag_i - 1;
641     sum     = x[i];
642     while (nz--) sum -= *v++ * x[*vi++];
643     x[i]    = sum*aa[adiag_i];
644   }
645 #endif
646   PetscLogFlops(2*a->nz - A->n);
647   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
648   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
649   PetscFunctionReturn(0);
650 }
651 
652 #undef __FUNC__
653 #define __FUNC__ "MatSolveAdd_SeqAIJ"
654 int MatSolveAdd_SeqAIJ(Mat A,Vec bb,Vec yy,Vec xx)
655 {
656   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
657   IS         iscol = a->col,isrow = a->row;
658   int        *r,*c,ierr,i, n = A->m,*vi,*ai = a->i,*aj = a->j;
659   int        nz,shift = a->indexshift,*rout,*cout;
660   Scalar     *x,*b,*tmp,*aa = a->a,sum,*v;
661 
662   PetscFunctionBegin;
663   if (yy != xx) {ierr = VecCopy(yy,xx);CHKERRQ(ierr);}
664 
665   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
666   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
667   tmp  = a->solve_work;
668 
669   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
670   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
671 
672   /* forward solve the lower triangular */
673   tmp[0] = b[*r++];
674   for (i=1; i<n; i++) {
675     v   = aa + ai[i] + shift;
676     vi  = aj + ai[i] + shift;
677     nz  = a->diag[i] - ai[i];
678     sum = b[*r++];
679     while (nz--) sum -= *v++ * tmp[*vi++ + shift];
680     tmp[i] = sum;
681   }
682 
683   /* backward solve the upper triangular */
684   for (i=n-1; i>=0; i--){
685     v   = aa + a->diag[i] + (!shift);
686     vi  = aj + a->diag[i] + (!shift);
687     nz  = ai[i+1] - a->diag[i] - 1;
688     sum = tmp[i];
689     while (nz--) sum -= *v++ * tmp[*vi++ + shift];
690     tmp[i] = sum*aa[a->diag[i]+shift];
691     x[*c--] += tmp[i];
692   }
693 
694   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
695   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
696   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
697   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
698   PetscLogFlops(2*a->nz);
699 
700   PetscFunctionReturn(0);
701 }
702 /* -------------------------------------------------------------------*/
703 #undef __FUNC__
704 #define __FUNC__ "MatSolveTranspose_SeqAIJ"
705 int MatSolveTranspose_SeqAIJ(Mat A,Vec bb,Vec xx)
706 {
707   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
708   IS         iscol = a->col,isrow = a->row;
709   int        *r,*c,ierr,i,n = A->m,*vi,*ai = a->i,*aj = a->j;
710   int        nz,shift = a->indexshift,*rout,*cout,*diag = a->diag;
711   Scalar     *x,*b,*tmp,*aa = a->a,*v,s1;
712 
713   PetscFunctionBegin;
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;
720 
721   /* copy the b into temp work space according to permutation */
722   for (i=0; i<n; i++) tmp[i] = b[c[i]];
723 
724   /* forward solve the U^T */
725   for (i=0; i<n; i++) {
726     v   = aa + diag[i] + shift;
727     vi  = aj + diag[i] + (!shift);
728     nz  = ai[i+1] - diag[i] - 1;
729     s1  = tmp[i];
730     s1 *= *(v++);  /* multiply by inverse of diagonal entry */
731     while (nz--) {
732       tmp[*vi++ + shift] -= (*v++)*s1;
733     }
734     tmp[i] = s1;
735   }
736 
737   /* backward solve the L^T */
738   for (i=n-1; i>=0; i--){
739     v   = aa + diag[i] - 1 + shift;
740     vi  = aj + diag[i] - 1 + shift;
741     nz  = diag[i] - ai[i];
742     s1  = tmp[i];
743     while (nz--) {
744       tmp[*vi-- + shift] -= (*v--)*s1;
745     }
746   }
747 
748   /* copy tmp into x according to permutation */
749   for (i=0; i<n; i++) x[r[i]] = tmp[i];
750 
751   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
752   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
753   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
754   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
755 
756   PetscLogFlops(2*a->nz-A->n);
757   PetscFunctionReturn(0);
758 }
759 
760 #undef __FUNC__
761 #define __FUNC__ "MatSolveTransposeAdd_SeqAIJ"
762 int MatSolveTransposeAdd_SeqAIJ(Mat A,Vec bb,Vec zz,Vec xx)
763 {
764   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
765   IS         iscol = a->col,isrow = a->row;
766   int        *r,*c,ierr,i,n = A->m,*vi,*ai = a->i,*aj = a->j;
767   int        nz,shift = a->indexshift,*rout,*cout,*diag = a->diag;
768   Scalar     *x,*b,*tmp,*aa = a->a,*v;
769 
770   PetscFunctionBegin;
771   if (zz != xx) VecCopy(zz,xx);
772 
773   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
774   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
775   tmp = a->solve_work;
776 
777   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
778   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
779 
780   /* copy the b into temp work space according to permutation */
781   for (i=0; i<n; i++) tmp[i] = b[c[i]];
782 
783   /* forward solve the U^T */
784   for (i=0; i<n; i++) {
785     v   = aa + diag[i] + shift;
786     vi  = aj + diag[i] + (!shift);
787     nz  = ai[i+1] - diag[i] - 1;
788     tmp[i] *= *v++;
789     while (nz--) {
790       tmp[*vi++ + shift] -= (*v++)*tmp[i];
791     }
792   }
793 
794   /* backward solve the L^T */
795   for (i=n-1; i>=0; i--){
796     v   = aa + diag[i] - 1 + shift;
797     vi  = aj + diag[i] - 1 + shift;
798     nz  = diag[i] - ai[i];
799     while (nz--) {
800       tmp[*vi-- + shift] -= (*v--)*tmp[i];
801     }
802   }
803 
804   /* copy tmp into x according to permutation */
805   for (i=0; i<n; i++) x[r[i]] += tmp[i];
806 
807   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
808   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
809   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
810   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
811 
812   PetscLogFlops(2*a->nz);
813   PetscFunctionReturn(0);
814 }
815 /* ----------------------------------------------------------------*/
816 EXTERN int MatMissingDiagonal_SeqAIJ(Mat);
817 
818 #undef __FUNC__
819 #define __FUNC__ "MatILUFactorSymbolic_SeqAIJ"
820 int MatILUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatILUInfo *info,Mat *fact)
821 {
822   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b;
823   IS         isicol;
824   int        *r,*ic,ierr,prow,n = A->m,*ai = a->i,*aj = a->j;
825   int        *ainew,*ajnew,jmax,*fill,*xi,nz,*im,*ajfill,*flev;
826   int        *dloc,idx,row,m,fm,nzf,nzi,len, realloc = 0,dcount = 0;
827   int        incrlev,nnz,i,shift = a->indexshift,levels,diagonal_fill;
828   PetscTruth col_identity,row_identity;
829   PetscReal  f;
830 
831   PetscFunctionBegin;
832   if (info) {
833     f             = info->fill;
834     levels        = (int)info->levels;
835     diagonal_fill = (int)info->diagonal_fill;
836   } else {
837     f             = 1.0;
838     levels        = 0;
839     diagonal_fill = 0;
840   }
841   if (!isrow) {
842     ierr = ISCreateStride(PETSC_COMM_SELF,A->M,0,1,&isrow);CHKERRQ(ierr);
843   }
844   if (!iscol) {
845     ierr = ISCreateStride(PETSC_COMM_SELF,A->M,0,1,&iscol);CHKERRQ(ierr);
846   }
847 
848   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
849 
850   /* special case that simply copies fill pattern */
851   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
852   ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
853   if (!levels && row_identity && col_identity) {
854     ierr = MatDuplicate_SeqAIJ(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr);
855     (*fact)->factor = FACTOR_LU;
856     b               = (Mat_SeqAIJ*)(*fact)->data;
857     if (!b->diag) {
858       ierr = MatMarkDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
859     }
860     ierr = MatMissingDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
861     b->row              = isrow;
862     b->col              = iscol;
863     b->icol             = isicol;
864     if (info) {
865       b->lu_damp    = (PetscTruth)((int)info->damp);
866       b->lu_damping = info->damping;
867     } else {
868       b->lu_damp    = PETSC_FALSE;
869       b->lu_damping = 0.0;
870     }
871     ierr                = PetscMalloc(((*fact)->m+1)*sizeof(Scalar),&b->solve_work);CHKERRQ(ierr);
872     (*fact)->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
873     ierr                = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
874     ierr                = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
875     PetscFunctionReturn(0);
876   }
877 
878   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
879   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
880 
881   /* get new row pointers */
882   ierr = PetscMalloc((n+1)*sizeof(int),&ainew);CHKERRQ(ierr);
883   ainew[0] = -shift;
884   /* don't know how many column pointers are needed so estimate */
885   jmax = (int)(f*(ai[n]+!shift));
886   ierr = PetscMalloc((jmax)*sizeof(int),&ajnew);CHKERRQ(ierr);
887   /* ajfill is level of fill for each fill entry */
888   ierr = PetscMalloc((jmax)*sizeof(int),&ajfill);CHKERRQ(ierr);
889   /* fill is a linked list of nonzeros in active row */
890   ierr = PetscMalloc((n+1)*sizeof(int),&fill);CHKERRQ(ierr);
891   /* im is level for each filled value */
892   ierr = PetscMalloc((n+1)*sizeof(int),&im);CHKERRQ(ierr);
893   /* dloc is location of diagonal in factor */
894   ierr = PetscMalloc((n+1)*sizeof(int),&dloc);CHKERRQ(ierr);
895   dloc[0]  = 0;
896   for (prow=0; prow<n; prow++) {
897 
898     /* copy current row into linked list */
899     nzf     = nz  = ai[r[prow]+1] - ai[r[prow]];
900     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
901     xi      = aj + ai[r[prow]] + shift;
902     fill[n]    = n;
903     fill[prow] = -1; /* marker to indicate if diagonal exists */
904     while (nz--) {
905       fm  = n;
906       idx = ic[*xi++ + shift];
907       do {
908         m  = fm;
909         fm = fill[m];
910       } while (fm < idx);
911       fill[m]   = idx;
912       fill[idx] = fm;
913       im[idx]   = 0;
914     }
915 
916     /* make sure diagonal entry is included */
917     if (diagonal_fill && fill[prow] == -1) {
918       fm = n;
919       while (fill[fm] < prow) fm = fill[fm];
920       fill[prow] = fill[fm]; /* insert diagonal into linked list */
921       fill[fm]   = prow;
922       im[prow]   = 0;
923       nzf++;
924       dcount++;
925     }
926 
927     nzi = 0;
928     row = fill[n];
929     while (row < prow) {
930       incrlev = im[row] + 1;
931       nz      = dloc[row];
932       xi      = ajnew  + ainew[row] + shift + nz + 1;
933       flev    = ajfill + ainew[row] + shift + nz + 1;
934       nnz     = ainew[row+1] - ainew[row] - nz - 1;
935       fm      = row;
936       while (nnz-- > 0) {
937         idx = *xi++ + shift;
938         if (*flev + incrlev > levels) {
939           flev++;
940           continue;
941         }
942         do {
943           m  = fm;
944           fm = fill[m];
945         } while (fm < idx);
946         if (fm != idx) {
947           im[idx]   = *flev + incrlev;
948           fill[m]   = idx;
949           fill[idx] = fm;
950           fm        = idx;
951           nzf++;
952         } else {
953           if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev;
954         }
955         flev++;
956       }
957       row = fill[row];
958       nzi++;
959     }
960     /* copy new filled row into permanent storage */
961     ainew[prow+1] = ainew[prow] + nzf;
962     if (ainew[prow+1] > jmax-shift) {
963 
964       /* estimate how much additional space we will need */
965       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
966       /* just double the memory each time */
967       /*  maxadd = (int)((f*(ai[n]+!shift)*(n-prow+5))/n); */
968       int maxadd = jmax;
969       if (maxadd < nzf) maxadd = (n-prow)*(nzf+1);
970       jmax += maxadd;
971 
972       /* allocate a longer ajnew and ajfill */
973       ierr   = PetscMalloc(jmax*sizeof(int),&xi);CHKERRQ(ierr);
974       ierr   = PetscMemcpy(xi,ajnew,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr);
975       ierr   = PetscFree(ajnew);CHKERRQ(ierr);
976       ajnew  = xi;
977       ierr   = PetscMalloc(jmax*sizeof(int),&xi);CHKERRQ(ierr);
978       ierr   = PetscMemcpy(xi,ajfill,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr);
979       ierr   = PetscFree(ajfill);CHKERRQ(ierr);
980       ajfill = xi;
981       realloc++; /* count how many times we realloc */
982     }
983     xi          = ajnew + ainew[prow] + shift;
984     flev        = ajfill + ainew[prow] + shift;
985     dloc[prow]  = nzi;
986     fm          = fill[n];
987     while (nzf--) {
988       *xi++   = fm - shift;
989       *flev++ = im[fm];
990       fm      = fill[fm];
991     }
992     /* make sure row has diagonal entry */
993     if (ajnew[ainew[prow]+shift+dloc[prow]]+shift != prow) {
994       SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,"Row %d has missing diagonal in factored matrix\n\
995     try running with -pc_ilu_nonzeros_along_diagonal or -pc_ilu_diagonal_fill",prow);
996     }
997   }
998   ierr = PetscFree(ajfill);CHKERRQ(ierr);
999   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
1000   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
1001   ierr = PetscFree(fill);CHKERRQ(ierr);
1002   ierr = PetscFree(im);CHKERRQ(ierr);
1003 
1004   {
1005     PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
1006     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af);
1007     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Run with -pc_ilu_fill %g or use \n",af);
1008     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:PCILUSetFill(pc,%g);\n",af);
1009     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:for best performance.\n");
1010     if (diagonal_fill) {
1011       PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Detected and replace %d missing diagonals",dcount);
1012     }
1013   }
1014 
1015   /* put together the new matrix */
1016   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,fact);CHKERRQ(ierr);
1017   PetscLogObjectParent(*fact,isicol);
1018   b = (Mat_SeqAIJ*)(*fact)->data;
1019   ierr = PetscFree(b->imax);CHKERRQ(ierr);
1020   b->singlemalloc = PETSC_FALSE;
1021   len = (ainew[n] + shift)*sizeof(Scalar);
1022   /* the next line frees the default space generated by the Create() */
1023   ierr = PetscFree(b->a);CHKERRQ(ierr);
1024   ierr = PetscFree(b->ilen);CHKERRQ(ierr);
1025   ierr = PetscMalloc(len+1,&b->a);CHKERRQ(ierr);
1026   b->j          = ajnew;
1027   b->i          = ainew;
1028   for (i=0; i<n; i++) dloc[i] += ainew[i];
1029   b->diag       = dloc;
1030   b->ilen       = 0;
1031   b->imax       = 0;
1032   b->row        = isrow;
1033   b->col        = iscol;
1034   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
1035   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
1036   b->icol       = isicol;
1037   ierr = PetscMalloc((n+1)*sizeof(Scalar),&b->solve_work);CHKERRQ(ierr);
1038   /* In b structure:  Free imax, ilen, old a, old j.
1039      Allocate dloc, solve_work, new a, new j */
1040   PetscLogObjectMemory(*fact,(ainew[n]+shift-n) * (sizeof(int)+sizeof(Scalar)));
1041   b->maxnz          = b->nz = ainew[n] + shift;
1042   if (info) {
1043     b->lu_damp    = (PetscTruth)((int)info->damp);
1044     b->lu_damping = info->damping;
1045   } else {
1046     b->lu_damp    = PETSC_FALSE;
1047     b->lu_damping = 0.0;
1048   }
1049   (*fact)->factor   = FACTOR_LU;
1050 
1051   (*fact)->info.factor_mallocs    = realloc;
1052   (*fact)->info.fill_ratio_given  = f;
1053   (*fact)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[prow]);
1054   (*fact)->factor                 =  FACTOR_LU;
1055 
1056   PetscFunctionReturn(0);
1057 }
1058 
1059 
1060 
1061 
1062