xref: /petsc/src/mat/impls/aij/seq/aijfact.c (revision 5ed2d596fd9914de76abbaed56c65f4792ae57ae) 
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 
7 #undef __FUNCT__
8 #define __FUNCT__ "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,PetscTruth);
22 
23 EXTERN int SPARSEKIT2dperm(int*,PetscScalar*,int*,int*,PetscScalar*,int*,int*,int*,int*,int*);
24 EXTERN int SPARSEKIT2ilutp(int*,PetscScalar*,int*,int*,int*,PetscReal,PetscReal*,int*,PetscScalar*,int*,int*,int*,PetscScalar*,int*,int*,int*);
25 EXTERN int SPARSEKIT2msrcsr(int*,PetscScalar*,int*,PetscScalar*,int*,int*,PetscScalar*,int*);
26 
27 #undef __FUNCT__
28 #define __FUNCT__ "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   PetscScalar  *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(PetscScalar),&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(PetscScalar),&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(PetscScalar),&w);CHKERRQ(ierr);
145 
146   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);
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(PetscScalar),&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(PetscScalar),&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,PETSC_FALSE);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 __FUNCT__
267 #define __FUNCT__ "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   if (A->M != A->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
279   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
280   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
281   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
282 
283   /* get new row pointers */
284   ierr = PetscMalloc((n+1)*sizeof(int),&ainew);CHKERRQ(ierr);
285   ainew[0] = -shift;
286   /* don't know how many column pointers are needed so estimate */
287   if (info) f = info->fill; else f = 1.0;
288   jmax  = (int)(f*ai[n]+(!shift));
289   ierr = PetscMalloc((jmax)*sizeof(int),&ajnew);CHKERRQ(ierr);
290   /* fill is a linked list of nonzeros in active row */
291   ierr = PetscMalloc((2*n+1)*sizeof(int),&fill);CHKERRQ(ierr);
292   im   = fill + n + 1;
293   /* idnew is location of diagonal in factor */
294   ierr = PetscMalloc((n+1)*sizeof(int),&idnew);CHKERRQ(ierr);
295   idnew[0] = -shift;
296 
297   for (i=0; i<n; i++) {
298     /* first copy previous fill into linked list */
299     nnz     = nz    = ai[r[i]+1] - ai[r[i]];
300     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
301     ajtmp   = aj + ai[r[i]] + shift;
302     fill[n] = n;
303     while (nz--) {
304       fm  = n;
305       idx = ic[*ajtmp++ + shift];
306       do {
307         m  = fm;
308         fm = fill[m];
309       } while (fm < idx);
310       fill[m]   = idx;
311       fill[idx] = fm;
312     }
313     row = fill[n];
314     while (row < i) {
315       ajtmp = ajnew + idnew[row] + (!shift);
316       nzbd  = 1 + idnew[row] - ainew[row];
317       nz    = im[row] - nzbd;
318       fm    = row;
319       while (nz-- > 0) {
320         idx = *ajtmp++ + shift;
321         nzbd++;
322         if (idx == i) im[row] = nzbd;
323         do {
324           m  = fm;
325           fm = fill[m];
326         } while (fm < idx);
327         if (fm != idx) {
328           fill[m]   = idx;
329           fill[idx] = fm;
330           fm        = idx;
331           nnz++;
332         }
333       }
334       row = fill[row];
335     }
336     /* copy new filled row into permanent storage */
337     ainew[i+1] = ainew[i] + nnz;
338     if (ainew[i+1] > jmax) {
339 
340       /* estimate how much additional space we will need */
341       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
342       /* just double the memory each time */
343       int maxadd = jmax;
344       /* maxadd = (int)((f*(ai[n]+(!shift))*(n-i+5))/n); */
345       if (maxadd < nnz) maxadd = (n-i)*(nnz+1);
346       jmax += maxadd;
347 
348       /* allocate a longer ajnew */
349       ierr = PetscMalloc(jmax*sizeof(int),&ajtmp);CHKERRQ(ierr);
350       ierr  = PetscMemcpy(ajtmp,ajnew,(ainew[i]+shift)*sizeof(int));CHKERRQ(ierr);
351       ierr  = PetscFree(ajnew);CHKERRQ(ierr);
352       ajnew = ajtmp;
353       realloc++; /* count how many times we realloc */
354     }
355     ajtmp = ajnew + ainew[i] + shift;
356     fm    = fill[n];
357     nzi   = 0;
358     im[i] = nnz;
359     while (nnz--) {
360       if (fm < i) nzi++;
361       *ajtmp++ = fm - shift;
362       fm       = fill[fm];
363     }
364     idnew[i] = ainew[i] + nzi;
365   }
366   if (ai[n] != 0) {
367     PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
368     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af);
369     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Run with -pc_lu_fill %g or use \n",af);
370     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:PCLUSetFill(pc,%g);\n",af);
371     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:for best performance.\n");
372   } else {
373     PetscLogInfo(A,"MatLUFactorSymbolic_SeqAIJ: Empty matrix\n");
374   }
375 
376   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
377   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
378 
379   ierr = PetscFree(fill);CHKERRQ(ierr);
380 
381   /* put together the new matrix */
382   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,B);CHKERRQ(ierr);
383   PetscLogObjectParent(*B,isicol);
384   b = (Mat_SeqAIJ*)(*B)->data;
385   ierr = PetscFree(b->imax);CHKERRQ(ierr);
386   b->singlemalloc = PETSC_FALSE;
387   /* the next line frees the default space generated by the Create() */
388   ierr = PetscFree(b->a);CHKERRQ(ierr);
389   ierr = PetscFree(b->ilen);CHKERRQ(ierr);
390   ierr          = PetscMalloc((ainew[n]+shift+1)*sizeof(PetscScalar),&b->a);CHKERRQ(ierr);
391   b->j          = ajnew;
392   b->i          = ainew;
393   b->diag       = idnew;
394   b->ilen       = 0;
395   b->imax       = 0;
396   b->row        = isrow;
397   b->col        = iscol;
398   if (info) {
399     b->lu_damp      = (PetscTruth) ((int)info->damp);
400     b->lu_damping   = info->damping;
401     b->lu_zeropivot = info->zeropivot;
402   } else {
403     b->lu_damp      = PETSC_FALSE;
404     b->lu_damping   = 0.0;
405     b->lu_zeropivot = 1.e-12;
406   }
407   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
408   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
409   b->icol       = isicol;
410   ierr          = PetscMalloc((n+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
411   /* In b structure:  Free imax, ilen, old a, old j.
412      Allocate idnew, solve_work, new a, new j */
413   PetscLogObjectMemory(*B,(ainew[n]+shift-n)*(sizeof(int)+sizeof(PetscScalar)));
414   b->maxnz = b->nz = ainew[n] + shift;
415 
416   (*B)->factor                 =  FACTOR_LU;
417   (*B)->info.factor_mallocs    = realloc;
418   (*B)->info.fill_ratio_given  = f;
419   ierr = Mat_AIJ_CheckInode(*B,PETSC_FALSE);CHKERRQ(ierr);
420   (*B)->ops->lufactornumeric   =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
421 
422   if (ai[n] != 0) {
423     (*B)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
424   } else {
425     (*B)->info.fill_ratio_needed = 0.0;
426   }
427   PetscFunctionReturn(0);
428 }
429 /* ----------------------------------------------------------- */
430 EXTERN int Mat_AIJ_CheckInode(Mat,PetscTruth);
431 
432 #undef __FUNCT__
433 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJ"
434 int MatLUFactorNumeric_SeqAIJ(Mat A,Mat *B)
435 {
436   Mat          C = *B;
437   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ *)C->data;
438   IS           isrow = b->row,isicol = b->icol;
439   int          *r,*ic,ierr,i,j,n = A->m,*ai = b->i,*aj = b->j;
440   int          *ajtmpold,*ajtmp,nz,row,*ics,shift = a->indexshift;
441   int          *diag_offset = b->diag,diag;
442   int          *pj,ndamp = 0;
443   PetscScalar  *rtmp,*v,*pc,multiplier,*pv,*rtmps;
444   PetscReal    damping = b->lu_damping, zeropivot = b->lu_zeropivot;
445   PetscTruth   damp;
446 
447   PetscFunctionBegin;
448   ierr  = ISGetIndices(isrow,&r);CHKERRQ(ierr);
449   ierr  = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
450   ierr  = PetscMalloc((n+1)*sizeof(PetscScalar),&rtmp);CHKERRQ(ierr);
451   ierr  = PetscMemzero(rtmp,(n+1)*sizeof(PetscScalar));CHKERRQ(ierr);
452   rtmps = rtmp + shift; ics = ic + shift;
453 
454   do {
455     damp = PETSC_FALSE;
456     for (i=0; i<n; i++) {
457       nz    = ai[i+1] - ai[i];
458       ajtmp = aj + ai[i] + shift;
459       for  (j=0; j<nz; j++) rtmps[ajtmp[j]] = 0.0;
460 
461       /* load in initial (unfactored row) */
462       nz       = a->i[r[i]+1] - a->i[r[i]];
463       ajtmpold = a->j + a->i[r[i]] + shift;
464       v        = a->a + a->i[r[i]] + shift;
465       for (j=0; j<nz; j++) {
466         rtmp[ics[ajtmpold[j]]] = v[j];
467         if (ajtmpold[j] == r[i]) { /* damp the diagonal of the matrix */
468           rtmp[ics[ajtmpold[j]]] += damping;
469         }
470       }
471 
472       row = *ajtmp++ + shift;
473       while  (row < i) {
474         pc = rtmp + row;
475         if (*pc != 0.0) {
476           pv         = b->a + diag_offset[row] + shift;
477           pj         = b->j + diag_offset[row] + (!shift);
478           multiplier = *pc / *pv++;
479           *pc        = multiplier;
480           nz         = ai[row+1] - diag_offset[row] - 1;
481           for (j=0; j<nz; j++) rtmps[pj[j]] -= multiplier * pv[j];
482           PetscLogFlops(2*nz);
483         }
484         row = *ajtmp++ + shift;
485       }
486       /* finished row so stick it into b->a */
487       pv = b->a + ai[i] + shift;
488       pj = b->j + ai[i] + shift;
489       nz = ai[i+1] - ai[i];
490       for (j=0; j<nz; j++) {pv[j] = rtmps[pj[j]];}
491       diag = diag_offset[i] - ai[i];
492       if (PetscAbsScalar(pv[diag]) < zeropivot) {
493         if (b->lu_damp) {
494           damp = PETSC_TRUE;
495           if (damping) damping *= 2.0;
496           else damping = 1.e-12;
497           ndamp++;
498           break;
499         } else {
500           SETERRQ3(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %d value %g tolerance %g",i,PetscAbsScalar(pv[diag]),zeropivot);
501         }
502       }
503     }
504   } while (damp);
505 
506   /* invert diagonal entries for simplier triangular solves */
507   for (i=0; i<n; i++) {
508     b->a[diag_offset[i]+shift] = 1.0/b->a[diag_offset[i]+shift];
509   }
510 
511   ierr = PetscFree(rtmp);CHKERRQ(ierr);
512   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
513   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
514   C->factor = FACTOR_LU;
515   (*B)->ops->lufactornumeric   =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
516   C->assembled = PETSC_TRUE;
517   PetscLogFlops(C->n);
518   if (ndamp || b->lu_damping) {
519     PetscLogInfo(0,"MatLUFactorNumerical_SeqAIJ: number of damping tries %d damping value %g\n",ndamp,damping);
520   }
521   PetscFunctionReturn(0);
522 }
523 /* ----------------------------------------------------------- */
524 #undef __FUNCT__
525 #define __FUNCT__ "MatLUFactor_SeqAIJ"
526 int MatLUFactor_SeqAIJ(Mat A,IS row,IS col,MatLUInfo *info)
527 {
528   int ierr;
529   Mat C;
530 
531   PetscFunctionBegin;
532   ierr = MatLUFactorSymbolic(A,row,col,info,&C);CHKERRQ(ierr);
533   ierr = MatLUFactorNumeric(A,&C);CHKERRQ(ierr);
534   ierr = MatHeaderCopy(A,C);CHKERRQ(ierr);
535   PetscLogObjectParent(A,((Mat_SeqAIJ*)(A->data))->icol);
536   PetscFunctionReturn(0);
537 }
538 /* ----------------------------------------------------------- */
539 #undef __FUNCT__
540 #define __FUNCT__ "MatSolve_SeqAIJ"
541 int MatSolve_SeqAIJ(Mat A,Vec bb,Vec xx)
542 {
543   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
544   IS           iscol = a->col,isrow = a->row;
545   int          *r,*c,ierr,i, n = A->m,*vi,*ai = a->i,*aj = a->j;
546   int          nz,shift = a->indexshift,*rout,*cout;
547   PetscScalar  *x,*b,*tmp,*tmps,*aa = a->a,sum,*v;
548 
549   PetscFunctionBegin;
550   if (!n) PetscFunctionReturn(0);
551 
552   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
553   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
554   tmp  = a->solve_work;
555 
556   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
557   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
558 
559   /* forward solve the lower triangular */
560   tmp[0] = b[*r++];
561   tmps   = tmp + shift;
562   for (i=1; i<n; i++) {
563     v   = aa + ai[i] + shift;
564     vi  = aj + ai[i] + shift;
565     nz  = a->diag[i] - ai[i];
566     sum = b[*r++];
567     while (nz--) sum -= *v++ * tmps[*vi++];
568     tmp[i] = sum;
569   }
570 
571   /* backward solve the upper triangular */
572   for (i=n-1; i>=0; i--){
573     v   = aa + a->diag[i] + (!shift);
574     vi  = aj + a->diag[i] + (!shift);
575     nz  = ai[i+1] - a->diag[i] - 1;
576     sum = tmp[i];
577     while (nz--) sum -= *v++ * tmps[*vi++];
578     x[*c--] = tmp[i] = sum*aa[a->diag[i]+shift];
579   }
580 
581   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
582   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
583   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
584   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
585   PetscLogFlops(2*a->nz - A->n);
586   PetscFunctionReturn(0);
587 }
588 
589 /* ----------------------------------------------------------- */
590 #undef __FUNCT__
591 #define __FUNCT__ "MatSolve_SeqAIJ_NaturalOrdering"
592 int MatSolve_SeqAIJ_NaturalOrdering(Mat A,Vec bb,Vec xx)
593 {
594   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
595   int          n = A->m,*ai = a->i,*aj = a->j,*adiag = a->diag,ierr;
596   PetscScalar  *x,*b,*aa = a->a,sum;
597 #if !defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
598   int          adiag_i,i,*vi,nz,ai_i;
599   PetscScalar  *v;
600 #endif
601 
602   PetscFunctionBegin;
603   if (!n) PetscFunctionReturn(0);
604   if (a->indexshift) {
605      ierr = MatSolve_SeqAIJ(A,bb,xx);CHKERRQ(ierr);
606      PetscFunctionReturn(0);
607   }
608 
609   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
610   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
611 
612 #if defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
613   fortransolveaij_(&n,x,ai,aj,adiag,aa,b);
614 #else
615   /* forward solve the lower triangular */
616   x[0] = b[0];
617   for (i=1; i<n; i++) {
618     ai_i = ai[i];
619     v    = aa + ai_i;
620     vi   = aj + ai_i;
621     nz   = adiag[i] - ai_i;
622     sum  = b[i];
623     while (nz--) sum -= *v++ * x[*vi++];
624     x[i] = sum;
625   }
626 
627   /* backward solve the upper triangular */
628   for (i=n-1; i>=0; i--){
629     adiag_i = adiag[i];
630     v       = aa + adiag_i + 1;
631     vi      = aj + adiag_i + 1;
632     nz      = ai[i+1] - adiag_i - 1;
633     sum     = x[i];
634     while (nz--) sum -= *v++ * x[*vi++];
635     x[i]    = sum*aa[adiag_i];
636   }
637 #endif
638   PetscLogFlops(2*a->nz - A->n);
639   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
640   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
641   PetscFunctionReturn(0);
642 }
643 
644 #undef __FUNCT__
645 #define __FUNCT__ "MatSolveAdd_SeqAIJ"
646 int MatSolveAdd_SeqAIJ(Mat A,Vec bb,Vec yy,Vec xx)
647 {
648   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
649   IS           iscol = a->col,isrow = a->row;
650   int          *r,*c,ierr,i, n = A->m,*vi,*ai = a->i,*aj = a->j;
651   int          nz,shift = a->indexshift,*rout,*cout;
652   PetscScalar  *x,*b,*tmp,*aa = a->a,sum,*v;
653 
654   PetscFunctionBegin;
655   if (yy != xx) {ierr = VecCopy(yy,xx);CHKERRQ(ierr);}
656 
657   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
658   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
659   tmp  = a->solve_work;
660 
661   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
662   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
663 
664   /* forward solve the lower triangular */
665   tmp[0] = b[*r++];
666   for (i=1; i<n; i++) {
667     v   = aa + ai[i] + shift;
668     vi  = aj + ai[i] + shift;
669     nz  = a->diag[i] - ai[i];
670     sum = b[*r++];
671     while (nz--) sum -= *v++ * tmp[*vi++ + shift];
672     tmp[i] = sum;
673   }
674 
675   /* backward solve the upper triangular */
676   for (i=n-1; i>=0; i--){
677     v   = aa + a->diag[i] + (!shift);
678     vi  = aj + a->diag[i] + (!shift);
679     nz  = ai[i+1] - a->diag[i] - 1;
680     sum = tmp[i];
681     while (nz--) sum -= *v++ * tmp[*vi++ + shift];
682     tmp[i] = sum*aa[a->diag[i]+shift];
683     x[*c--] += tmp[i];
684   }
685 
686   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
687   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
688   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
689   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
690   PetscLogFlops(2*a->nz);
691 
692   PetscFunctionReturn(0);
693 }
694 /* -------------------------------------------------------------------*/
695 #undef __FUNCT__
696 #define __FUNCT__ "MatSolveTranspose_SeqAIJ"
697 int MatSolveTranspose_SeqAIJ(Mat A,Vec bb,Vec xx)
698 {
699   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
700   IS           iscol = a->col,isrow = a->row;
701   int          *r,*c,ierr,i,n = A->m,*vi,*ai = a->i,*aj = a->j;
702   int          nz,shift = a->indexshift,*rout,*cout,*diag = a->diag;
703   PetscScalar  *x,*b,*tmp,*aa = a->a,*v,s1;
704 
705   PetscFunctionBegin;
706   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
707   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
708   tmp  = a->solve_work;
709 
710   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
711   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
712 
713   /* copy the b into temp work space according to permutation */
714   for (i=0; i<n; i++) tmp[i] = b[c[i]];
715 
716   /* forward solve the U^T */
717   for (i=0; i<n; i++) {
718     v   = aa + diag[i] + shift;
719     vi  = aj + diag[i] + (!shift);
720     nz  = ai[i+1] - diag[i] - 1;
721     s1  = tmp[i];
722     s1 *= (*v++);  /* multiply by inverse of diagonal entry */
723     while (nz--) {
724       tmp[*vi++ + shift] -= (*v++)*s1;
725     }
726     tmp[i] = s1;
727   }
728 
729   /* backward solve the L^T */
730   for (i=n-1; i>=0; i--){
731     v   = aa + diag[i] - 1 + shift;
732     vi  = aj + diag[i] - 1 + shift;
733     nz  = diag[i] - ai[i];
734     s1  = tmp[i];
735     while (nz--) {
736       tmp[*vi-- + shift] -= (*v--)*s1;
737     }
738   }
739 
740   /* copy tmp into x according to permutation */
741   for (i=0; i<n; i++) x[r[i]] = tmp[i];
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 
748   PetscLogFlops(2*a->nz-A->n);
749   PetscFunctionReturn(0);
750 }
751 
752 #undef __FUNCT__
753 #define __FUNCT__ "MatSolveTransposeAdd_SeqAIJ"
754 int MatSolveTransposeAdd_SeqAIJ(Mat A,Vec bb,Vec zz,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,shift = a->indexshift,*rout,*cout,*diag = a->diag;
760   PetscScalar  *x,*b,*tmp,*aa = a->a,*v;
761 
762   PetscFunctionBegin;
763   if (zz != xx) VecCopy(zz,xx);
764 
765   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
766   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
767   tmp = a->solve_work;
768 
769   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
770   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
771 
772   /* copy the b into temp work space according to permutation */
773   for (i=0; i<n; i++) tmp[i] = b[c[i]];
774 
775   /* forward solve the U^T */
776   for (i=0; i<n; i++) {
777     v   = aa + diag[i] + shift;
778     vi  = aj + diag[i] + (!shift);
779     nz  = ai[i+1] - diag[i] - 1;
780     tmp[i] *= *v++;
781     while (nz--) {
782       tmp[*vi++ + shift] -= (*v++)*tmp[i];
783     }
784   }
785 
786   /* backward solve the L^T */
787   for (i=n-1; i>=0; i--){
788     v   = aa + diag[i] - 1 + shift;
789     vi  = aj + diag[i] - 1 + shift;
790     nz  = diag[i] - ai[i];
791     while (nz--) {
792       tmp[*vi-- + shift] -= (*v--)*tmp[i];
793     }
794   }
795 
796   /* copy tmp into x according to permutation */
797   for (i=0; i<n; i++) x[r[i]] += tmp[i];
798 
799   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
800   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
801   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
802   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
803 
804   PetscLogFlops(2*a->nz);
805   PetscFunctionReturn(0);
806 }
807 /* ----------------------------------------------------------------*/
808 EXTERN int MatMissingDiagonal_SeqAIJ(Mat);
809 
810 #undef __FUNCT__
811 #define __FUNCT__ "MatILUFactorSymbolic_SeqAIJ"
812 int MatILUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatILUInfo *info,Mat *fact)
813 {
814   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b;
815   IS         isicol;
816   int        *r,*ic,ierr,prow,n = A->m,*ai = a->i,*aj = a->j;
817   int        *ainew,*ajnew,jmax,*fill,*xi,nz,*im,*ajfill,*flev;
818   int        *dloc,idx,row,m,fm,nzf,nzi,len, realloc = 0,dcount = 0;
819   int        incrlev,nnz,i,shift = a->indexshift,levels,diagonal_fill;
820   PetscTruth col_identity,row_identity;
821   PetscReal  f;
822 
823   PetscFunctionBegin;
824   if (info) {
825     f             = info->fill;
826     levels        = (int)info->levels;
827     diagonal_fill = (int)info->diagonal_fill;
828   } else {
829     f             = 1.0;
830     levels        = 0;
831     diagonal_fill = 0;
832   }
833   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
834 
835   /* special case that simply copies fill pattern */
836   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
837   ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
838   if (!levels && row_identity && col_identity) {
839     ierr = MatDuplicate_SeqAIJ(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr);
840     (*fact)->factor = FACTOR_LU;
841     b               = (Mat_SeqAIJ*)(*fact)->data;
842     if (!b->diag) {
843       ierr = MatMarkDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
844     }
845     ierr = MatMissingDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
846     b->row              = isrow;
847     b->col              = iscol;
848     b->icol             = isicol;
849     if (info) {
850       b->lu_damp      = (PetscTruth)((int)info->damp);
851       b->lu_damping   = info->damping;
852       b->lu_zeropivot = info->zeropivot;
853     } else {
854       b->lu_damp      = PETSC_FALSE;
855       b->lu_damping   = 0.0;
856       b->lu_zeropivot = 1.e-12;
857     }
858     ierr                = PetscMalloc(((*fact)->m+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
859     (*fact)->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
860     ierr                = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
861     ierr                = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
862     PetscFunctionReturn(0);
863   }
864 
865   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
866   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
867 
868   /* get new row pointers */
869   ierr = PetscMalloc((n+1)*sizeof(int),&ainew);CHKERRQ(ierr);
870   ainew[0] = -shift;
871   /* don't know how many column pointers are needed so estimate */
872   jmax = (int)(f*(ai[n]+!shift));
873   ierr = PetscMalloc((jmax)*sizeof(int),&ajnew);CHKERRQ(ierr);
874   /* ajfill is level of fill for each fill entry */
875   ierr = PetscMalloc((jmax)*sizeof(int),&ajfill);CHKERRQ(ierr);
876   /* fill is a linked list of nonzeros in active row */
877   ierr = PetscMalloc((n+1)*sizeof(int),&fill);CHKERRQ(ierr);
878   /* im is level for each filled value */
879   ierr = PetscMalloc((n+1)*sizeof(int),&im);CHKERRQ(ierr);
880   /* dloc is location of diagonal in factor */
881   ierr = PetscMalloc((n+1)*sizeof(int),&dloc);CHKERRQ(ierr);
882   dloc[0]  = 0;
883   for (prow=0; prow<n; prow++) {
884 
885     /* copy current row into linked list */
886     nzf     = nz  = ai[r[prow]+1] - ai[r[prow]];
887     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
888     xi      = aj + ai[r[prow]] + shift;
889     fill[n]    = n;
890     fill[prow] = -1; /* marker to indicate if diagonal exists */
891     while (nz--) {
892       fm  = n;
893       idx = ic[*xi++ + shift];
894       do {
895         m  = fm;
896         fm = fill[m];
897       } while (fm < idx);
898       fill[m]   = idx;
899       fill[idx] = fm;
900       im[idx]   = 0;
901     }
902 
903     /* make sure diagonal entry is included */
904     if (diagonal_fill && fill[prow] == -1) {
905       fm = n;
906       while (fill[fm] < prow) fm = fill[fm];
907       fill[prow] = fill[fm]; /* insert diagonal into linked list */
908       fill[fm]   = prow;
909       im[prow]   = 0;
910       nzf++;
911       dcount++;
912     }
913 
914     nzi = 0;
915     row = fill[n];
916     while (row < prow) {
917       incrlev = im[row] + 1;
918       nz      = dloc[row];
919       xi      = ajnew  + ainew[row] + shift + nz + 1;
920       flev    = ajfill + ainew[row] + shift + nz + 1;
921       nnz     = ainew[row+1] - ainew[row] - nz - 1;
922       fm      = row;
923       while (nnz-- > 0) {
924         idx = *xi++ + shift;
925         if (*flev + incrlev > levels) {
926           flev++;
927           continue;
928         }
929         do {
930           m  = fm;
931           fm = fill[m];
932         } while (fm < idx);
933         if (fm != idx) {
934           im[idx]   = *flev + incrlev;
935           fill[m]   = idx;
936           fill[idx] = fm;
937           fm        = idx;
938           nzf++;
939         } else {
940           if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev;
941         }
942         flev++;
943       }
944       row = fill[row];
945       nzi++;
946     }
947     /* copy new filled row into permanent storage */
948     ainew[prow+1] = ainew[prow] + nzf;
949     if (ainew[prow+1] > jmax-shift) {
950 
951       /* estimate how much additional space we will need */
952       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
953       /* just double the memory each time */
954       /*  maxadd = (int)((f*(ai[n]+!shift)*(n-prow+5))/n); */
955       int maxadd = jmax;
956       if (maxadd < nzf) maxadd = (n-prow)*(nzf+1);
957       jmax += maxadd;
958 
959       /* allocate a longer ajnew and ajfill */
960       ierr   = PetscMalloc(jmax*sizeof(int),&xi);CHKERRQ(ierr);
961       ierr   = PetscMemcpy(xi,ajnew,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr);
962       ierr   = PetscFree(ajnew);CHKERRQ(ierr);
963       ajnew  = xi;
964       ierr   = PetscMalloc(jmax*sizeof(int),&xi);CHKERRQ(ierr);
965       ierr   = PetscMemcpy(xi,ajfill,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr);
966       ierr   = PetscFree(ajfill);CHKERRQ(ierr);
967       ajfill = xi;
968       realloc++; /* count how many times we realloc */
969     }
970     xi          = ajnew + ainew[prow] + shift;
971     flev        = ajfill + ainew[prow] + shift;
972     dloc[prow]  = nzi;
973     fm          = fill[n];
974     while (nzf--) {
975       *xi++   = fm - shift;
976       *flev++ = im[fm];
977       fm      = fill[fm];
978     }
979     /* make sure row has diagonal entry */
980     if (ajnew[ainew[prow]+shift+dloc[prow]]+shift != prow) {
981       SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,"Row %d has missing diagonal in factored matrix\n\
982     try running with -pc_ilu_nonzeros_along_diagonal or -pc_ilu_diagonal_fill",prow);
983     }
984   }
985   ierr = PetscFree(ajfill);CHKERRQ(ierr);
986   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
987   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
988   ierr = PetscFree(fill);CHKERRQ(ierr);
989   ierr = PetscFree(im);CHKERRQ(ierr);
990 
991   {
992     PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
993     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af);
994     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Run with -[sub_]pc_ilu_fill %g or use \n",af);
995     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:PCILUSetFill([sub]pc,%g);\n",af);
996     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:for best performance.\n");
997     if (diagonal_fill) {
998       PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Detected and replaced %d missing diagonals",dcount);
999     }
1000   }
1001 
1002   /* put together the new matrix */
1003   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,fact);CHKERRQ(ierr);
1004   PetscLogObjectParent(*fact,isicol);
1005   b = (Mat_SeqAIJ*)(*fact)->data;
1006   ierr = PetscFree(b->imax);CHKERRQ(ierr);
1007   b->singlemalloc = PETSC_FALSE;
1008   len = (ainew[n] + shift)*sizeof(PetscScalar);
1009   /* the next line frees the default space generated by the Create() */
1010   ierr = PetscFree(b->a);CHKERRQ(ierr);
1011   ierr = PetscFree(b->ilen);CHKERRQ(ierr);
1012   ierr = PetscMalloc(len+1,&b->a);CHKERRQ(ierr);
1013   b->j          = ajnew;
1014   b->i          = ainew;
1015   for (i=0; i<n; i++) dloc[i] += ainew[i];
1016   b->diag       = dloc;
1017   b->ilen       = 0;
1018   b->imax       = 0;
1019   b->row        = isrow;
1020   b->col        = iscol;
1021   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
1022   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
1023   b->icol       = isicol;
1024   ierr = PetscMalloc((n+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
1025   /* In b structure:  Free imax, ilen, old a, old j.
1026      Allocate dloc, solve_work, new a, new j */
1027   PetscLogObjectMemory(*fact,(ainew[n]+shift-n) * (sizeof(int)+sizeof(PetscScalar)));
1028   b->maxnz          = b->nz = ainew[n] + shift;
1029   if (info) {
1030     b->lu_damp      = (PetscTruth)((int)info->damp);
1031     b->lu_damping   = info->damping;
1032     b->lu_zeropivot = info->zeropivot;
1033   } else {
1034     b->lu_damp      = PETSC_FALSE;
1035     b->lu_damping   = 0.0;
1036     b->lu_zeropivot = 1.e-12;
1037   }
1038   (*fact)->factor   = FACTOR_LU;
1039   ierr = Mat_AIJ_CheckInode(*fact,PETSC_FALSE);CHKERRQ(ierr);
1040   (*fact)->ops->lufactornumeric   =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
1041 
1042   (*fact)->info.factor_mallocs    = realloc;
1043   (*fact)->info.fill_ratio_given  = f;
1044   (*fact)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[prow]);
1045   (*fact)->factor                 =  FACTOR_LU;
1046   PetscFunctionReturn(0);
1047 }
1048 
1049 
1050 
1051 
1052