xref: /petsc/src/mat/impls/aij/seq/aijfact.c (revision ef66eb6987ddfdf4e414d6b820cbc8d8d7d17bc2) 
1 /*$Id: aijfact.c,v 1.165 2001/08/07 03:02:47 balay 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 __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   PetscFunctionReturn(0);
536 }
537 /* ----------------------------------------------------------- */
538 #undef __FUNCT__
539 #define __FUNCT__ "MatSolve_SeqAIJ"
540 int MatSolve_SeqAIJ(Mat A,Vec bb,Vec xx)
541 {
542   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
543   IS           iscol = a->col,isrow = a->row;
544   int          *r,*c,ierr,i, n = A->m,*vi,*ai = a->i,*aj = a->j;
545   int          nz,shift = a->indexshift,*rout,*cout;
546   PetscScalar  *x,*b,*tmp,*tmps,*aa = a->a,sum,*v;
547 
548   PetscFunctionBegin;
549   if (!n) PetscFunctionReturn(0);
550 
551   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
552   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
553   tmp  = a->solve_work;
554 
555   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
556   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
557 
558   /* forward solve the lower triangular */
559   tmp[0] = b[*r++];
560   tmps   = tmp + shift;
561   for (i=1; i<n; i++) {
562     v   = aa + ai[i] + shift;
563     vi  = aj + ai[i] + shift;
564     nz  = a->diag[i] - ai[i];
565     sum = b[*r++];
566     while (nz--) sum -= *v++ * tmps[*vi++];
567     tmp[i] = sum;
568   }
569 
570   /* backward solve the upper triangular */
571   for (i=n-1; i>=0; i--){
572     v   = aa + a->diag[i] + (!shift);
573     vi  = aj + a->diag[i] + (!shift);
574     nz  = ai[i+1] - a->diag[i] - 1;
575     sum = tmp[i];
576     while (nz--) sum -= *v++ * tmps[*vi++];
577     x[*c--] = tmp[i] = sum*aa[a->diag[i]+shift];
578   }
579 
580   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
581   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
582   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
583   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
584   PetscLogFlops(2*a->nz - A->n);
585   PetscFunctionReturn(0);
586 }
587 
588 /* ----------------------------------------------------------- */
589 #undef __FUNCT__
590 #define __FUNCT__ "MatSolve_SeqAIJ_NaturalOrdering"
591 int MatSolve_SeqAIJ_NaturalOrdering(Mat A,Vec bb,Vec xx)
592 {
593   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
594   int          n = A->m,*ai = a->i,*aj = a->j,*adiag = a->diag,ierr;
595   PetscScalar  *x,*b,*aa = a->a,sum;
596 #if !defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
597   int          adiag_i,i,*vi,nz,ai_i;
598   PetscScalar  *v;
599 #endif
600 
601   PetscFunctionBegin;
602   if (!n) PetscFunctionReturn(0);
603   if (a->indexshift) {
604      ierr = MatSolve_SeqAIJ(A,bb,xx);CHKERRQ(ierr);
605      PetscFunctionReturn(0);
606   }
607 
608   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
609   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
610 
611 #if defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
612   fortransolveaij_(&n,x,ai,aj,adiag,aa,b);
613 #else
614   /* forward solve the lower triangular */
615   x[0] = b[0];
616   for (i=1; i<n; i++) {
617     ai_i = ai[i];
618     v    = aa + ai_i;
619     vi   = aj + ai_i;
620     nz   = adiag[i] - ai_i;
621     sum  = b[i];
622     while (nz--) sum -= *v++ * x[*vi++];
623     x[i] = sum;
624   }
625 
626   /* backward solve the upper triangular */
627   for (i=n-1; i>=0; i--){
628     adiag_i = adiag[i];
629     v       = aa + adiag_i + 1;
630     vi      = aj + adiag_i + 1;
631     nz      = ai[i+1] - adiag_i - 1;
632     sum     = x[i];
633     while (nz--) sum -= *v++ * x[*vi++];
634     x[i]    = sum*aa[adiag_i];
635   }
636 #endif
637   PetscLogFlops(2*a->nz - A->n);
638   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
639   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
640   PetscFunctionReturn(0);
641 }
642 
643 #undef __FUNCT__
644 #define __FUNCT__ "MatSolveAdd_SeqAIJ"
645 int MatSolveAdd_SeqAIJ(Mat A,Vec bb,Vec yy,Vec xx)
646 {
647   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
648   IS           iscol = a->col,isrow = a->row;
649   int          *r,*c,ierr,i, n = A->m,*vi,*ai = a->i,*aj = a->j;
650   int          nz,shift = a->indexshift,*rout,*cout;
651   PetscScalar  *x,*b,*tmp,*aa = a->a,sum,*v;
652 
653   PetscFunctionBegin;
654   if (yy != xx) {ierr = VecCopy(yy,xx);CHKERRQ(ierr);}
655 
656   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
657   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
658   tmp  = a->solve_work;
659 
660   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
661   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
662 
663   /* forward solve the lower triangular */
664   tmp[0] = b[*r++];
665   for (i=1; i<n; i++) {
666     v   = aa + ai[i] + shift;
667     vi  = aj + ai[i] + shift;
668     nz  = a->diag[i] - ai[i];
669     sum = b[*r++];
670     while (nz--) sum -= *v++ * tmp[*vi++ + shift];
671     tmp[i] = sum;
672   }
673 
674   /* backward solve the upper triangular */
675   for (i=n-1; i>=0; i--){
676     v   = aa + a->diag[i] + (!shift);
677     vi  = aj + a->diag[i] + (!shift);
678     nz  = ai[i+1] - a->diag[i] - 1;
679     sum = tmp[i];
680     while (nz--) sum -= *v++ * tmp[*vi++ + shift];
681     tmp[i] = sum*aa[a->diag[i]+shift];
682     x[*c--] += tmp[i];
683   }
684 
685   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
686   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
687   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
688   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
689   PetscLogFlops(2*a->nz);
690 
691   PetscFunctionReturn(0);
692 }
693 /* -------------------------------------------------------------------*/
694 #undef __FUNCT__
695 #define __FUNCT__ "MatSolveTranspose_SeqAIJ"
696 int MatSolveTranspose_SeqAIJ(Mat A,Vec bb,Vec xx)
697 {
698   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
699   IS           iscol = a->col,isrow = a->row;
700   int          *r,*c,ierr,i,n = A->m,*vi,*ai = a->i,*aj = a->j;
701   int          nz,shift = a->indexshift,*rout,*cout,*diag = a->diag;
702   PetscScalar  *x,*b,*tmp,*aa = a->a,*v,s1;
703 
704   PetscFunctionBegin;
705   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
706   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
707   tmp  = a->solve_work;
708 
709   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
710   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
711 
712   /* copy the b into temp work space according to permutation */
713   for (i=0; i<n; i++) tmp[i] = b[c[i]];
714 
715   /* forward solve the U^T */
716   for (i=0; i<n; i++) {
717     v   = aa + diag[i] + shift;
718     vi  = aj + diag[i] + (!shift);
719     nz  = ai[i+1] - diag[i] - 1;
720     s1  = tmp[i];
721     s1 *= (*v++);  /* multiply by inverse of diagonal entry */
722     while (nz--) {
723       tmp[*vi++ + shift] -= (*v++)*s1;
724     }
725     tmp[i] = s1;
726   }
727 
728   /* backward solve the L^T */
729   for (i=n-1; i>=0; i--){
730     v   = aa + diag[i] - 1 + shift;
731     vi  = aj + diag[i] - 1 + shift;
732     nz  = diag[i] - ai[i];
733     s1  = tmp[i];
734     while (nz--) {
735       tmp[*vi-- + shift] -= (*v--)*s1;
736     }
737   }
738 
739   /* copy tmp into x according to permutation */
740   for (i=0; i<n; i++) x[r[i]] = tmp[i];
741 
742   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
743   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
744   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
745   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
746 
747   PetscLogFlops(2*a->nz-A->n);
748   PetscFunctionReturn(0);
749 }
750 
751 #undef __FUNCT__
752 #define __FUNCT__ "MatSolveTransposeAdd_SeqAIJ"
753 int MatSolveTransposeAdd_SeqAIJ(Mat A,Vec bb,Vec zz,Vec xx)
754 {
755   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
756   IS           iscol = a->col,isrow = a->row;
757   int          *r,*c,ierr,i,n = A->m,*vi,*ai = a->i,*aj = a->j;
758   int          nz,shift = a->indexshift,*rout,*cout,*diag = a->diag;
759   PetscScalar  *x,*b,*tmp,*aa = a->a,*v;
760 
761   PetscFunctionBegin;
762   if (zz != xx) VecCopy(zz,xx);
763 
764   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
765   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
766   tmp = a->solve_work;
767 
768   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
769   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
770 
771   /* copy the b into temp work space according to permutation */
772   for (i=0; i<n; i++) tmp[i] = b[c[i]];
773 
774   /* forward solve the U^T */
775   for (i=0; i<n; i++) {
776     v   = aa + diag[i] + shift;
777     vi  = aj + diag[i] + (!shift);
778     nz  = ai[i+1] - diag[i] - 1;
779     tmp[i] *= *v++;
780     while (nz--) {
781       tmp[*vi++ + shift] -= (*v++)*tmp[i];
782     }
783   }
784 
785   /* backward solve the L^T */
786   for (i=n-1; i>=0; i--){
787     v   = aa + diag[i] - 1 + shift;
788     vi  = aj + diag[i] - 1 + shift;
789     nz  = diag[i] - ai[i];
790     while (nz--) {
791       tmp[*vi-- + shift] -= (*v--)*tmp[i];
792     }
793   }
794 
795   /* copy tmp into x according to permutation */
796   for (i=0; i<n; i++) x[r[i]] += tmp[i];
797 
798   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
799   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
800   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
801   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
802 
803   PetscLogFlops(2*a->nz);
804   PetscFunctionReturn(0);
805 }
806 /* ----------------------------------------------------------------*/
807 EXTERN int MatMissingDiagonal_SeqAIJ(Mat);
808 
809 #undef __FUNCT__
810 #define __FUNCT__ "MatILUFactorSymbolic_SeqAIJ"
811 int MatILUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatILUInfo *info,Mat *fact)
812 {
813   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b;
814   IS         isicol;
815   int        *r,*ic,ierr,prow,n = A->m,*ai = a->i,*aj = a->j;
816   int        *ainew,*ajnew,jmax,*fill,*xi,nz,*im,*ajfill,*flev;
817   int        *dloc,idx,row,m,fm,nzf,nzi,len, realloc = 0,dcount = 0;
818   int        incrlev,nnz,i,shift = a->indexshift,levels,diagonal_fill;
819   PetscTruth col_identity,row_identity;
820   PetscReal  f;
821 
822   PetscFunctionBegin;
823   if (info) {
824     f             = info->fill;
825     levels        = (int)info->levels;
826     diagonal_fill = (int)info->diagonal_fill;
827   } else {
828     f             = 1.0;
829     levels        = 0;
830     diagonal_fill = 0;
831   }
832   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
833 
834   /* special case that simply copies fill pattern */
835   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
836   ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
837   if (!levels && row_identity && col_identity) {
838     ierr = MatDuplicate_SeqAIJ(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr);
839     (*fact)->factor = FACTOR_LU;
840     b               = (Mat_SeqAIJ*)(*fact)->data;
841     if (!b->diag) {
842       ierr = MatMarkDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
843     }
844     ierr = MatMissingDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
845     b->row              = isrow;
846     b->col              = iscol;
847     b->icol             = isicol;
848     if (info) {
849       b->lu_damp      = (PetscTruth)((int)info->damp);
850       b->lu_damping   = info->damping;
851       b->lu_zeropivot = info->zeropivot;
852     } else {
853       b->lu_damp      = PETSC_FALSE;
854       b->lu_damping   = 0.0;
855       b->lu_zeropivot = 1.e-12;
856     }
857     ierr                = PetscMalloc(((*fact)->m+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
858     (*fact)->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
859     ierr                = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
860     ierr                = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
861     PetscFunctionReturn(0);
862   }
863 
864   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
865   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
866 
867   /* get new row pointers */
868   ierr = PetscMalloc((n+1)*sizeof(int),&ainew);CHKERRQ(ierr);
869   ainew[0] = -shift;
870   /* don't know how many column pointers are needed so estimate */
871   jmax = (int)(f*(ai[n]+!shift));
872   ierr = PetscMalloc((jmax)*sizeof(int),&ajnew);CHKERRQ(ierr);
873   /* ajfill is level of fill for each fill entry */
874   ierr = PetscMalloc((jmax)*sizeof(int),&ajfill);CHKERRQ(ierr);
875   /* fill is a linked list of nonzeros in active row */
876   ierr = PetscMalloc((n+1)*sizeof(int),&fill);CHKERRQ(ierr);
877   /* im is level for each filled value */
878   ierr = PetscMalloc((n+1)*sizeof(int),&im);CHKERRQ(ierr);
879   /* dloc is location of diagonal in factor */
880   ierr = PetscMalloc((n+1)*sizeof(int),&dloc);CHKERRQ(ierr);
881   dloc[0]  = 0;
882   for (prow=0; prow<n; prow++) {
883 
884     /* copy current row into linked list */
885     nzf     = nz  = ai[r[prow]+1] - ai[r[prow]];
886     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
887     xi      = aj + ai[r[prow]] + shift;
888     fill[n]    = n;
889     fill[prow] = -1; /* marker to indicate if diagonal exists */
890     while (nz--) {
891       fm  = n;
892       idx = ic[*xi++ + shift];
893       do {
894         m  = fm;
895         fm = fill[m];
896       } while (fm < idx);
897       fill[m]   = idx;
898       fill[idx] = fm;
899       im[idx]   = 0;
900     }
901 
902     /* make sure diagonal entry is included */
903     if (diagonal_fill && fill[prow] == -1) {
904       fm = n;
905       while (fill[fm] < prow) fm = fill[fm];
906       fill[prow] = fill[fm]; /* insert diagonal into linked list */
907       fill[fm]   = prow;
908       im[prow]   = 0;
909       nzf++;
910       dcount++;
911     }
912 
913     nzi = 0;
914     row = fill[n];
915     while (row < prow) {
916       incrlev = im[row] + 1;
917       nz      = dloc[row];
918       xi      = ajnew  + ainew[row] + shift + nz + 1;
919       flev    = ajfill + ainew[row] + shift + nz + 1;
920       nnz     = ainew[row+1] - ainew[row] - nz - 1;
921       fm      = row;
922       while (nnz-- > 0) {
923         idx = *xi++ + shift;
924         if (*flev + incrlev > levels) {
925           flev++;
926           continue;
927         }
928         do {
929           m  = fm;
930           fm = fill[m];
931         } while (fm < idx);
932         if (fm != idx) {
933           im[idx]   = *flev + incrlev;
934           fill[m]   = idx;
935           fill[idx] = fm;
936           fm        = idx;
937           nzf++;
938         } else {
939           if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev;
940         }
941         flev++;
942       }
943       row = fill[row];
944       nzi++;
945     }
946     /* copy new filled row into permanent storage */
947     ainew[prow+1] = ainew[prow] + nzf;
948     if (ainew[prow+1] > jmax-shift) {
949 
950       /* estimate how much additional space we will need */
951       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
952       /* just double the memory each time */
953       /*  maxadd = (int)((f*(ai[n]+!shift)*(n-prow+5))/n); */
954       int maxadd = jmax;
955       if (maxadd < nzf) maxadd = (n-prow)*(nzf+1);
956       jmax += maxadd;
957 
958       /* allocate a longer ajnew and ajfill */
959       ierr   = PetscMalloc(jmax*sizeof(int),&xi);CHKERRQ(ierr);
960       ierr   = PetscMemcpy(xi,ajnew,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr);
961       ierr   = PetscFree(ajnew);CHKERRQ(ierr);
962       ajnew  = xi;
963       ierr   = PetscMalloc(jmax*sizeof(int),&xi);CHKERRQ(ierr);
964       ierr   = PetscMemcpy(xi,ajfill,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr);
965       ierr   = PetscFree(ajfill);CHKERRQ(ierr);
966       ajfill = xi;
967       realloc++; /* count how many times we realloc */
968     }
969     xi          = ajnew + ainew[prow] + shift;
970     flev        = ajfill + ainew[prow] + shift;
971     dloc[prow]  = nzi;
972     fm          = fill[n];
973     while (nzf--) {
974       *xi++   = fm - shift;
975       *flev++ = im[fm];
976       fm      = fill[fm];
977     }
978     /* make sure row has diagonal entry */
979     if (ajnew[ainew[prow]+shift+dloc[prow]]+shift != prow) {
980       SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,"Row %d has missing diagonal in factored matrix\n\
981     try running with -pc_ilu_nonzeros_along_diagonal or -pc_ilu_diagonal_fill",prow);
982     }
983   }
984   ierr = PetscFree(ajfill);CHKERRQ(ierr);
985   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
986   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
987   ierr = PetscFree(fill);CHKERRQ(ierr);
988   ierr = PetscFree(im);CHKERRQ(ierr);
989 
990   {
991     PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
992     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af);
993     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Run with -[sub_]pc_ilu_fill %g or use \n",af);
994     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:PCILUSetFill([sub]pc,%g);\n",af);
995     PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:for best performance.\n");
996     if (diagonal_fill) {
997       PetscLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Detected and replace %d missing diagonals",dcount);
998     }
999   }
1000 
1001   /* put together the new matrix */
1002   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,fact);CHKERRQ(ierr);
1003   PetscLogObjectParent(*fact,isicol);
1004   b = (Mat_SeqAIJ*)(*fact)->data;
1005   ierr = PetscFree(b->imax);CHKERRQ(ierr);
1006   b->singlemalloc = PETSC_FALSE;
1007   len = (ainew[n] + shift)*sizeof(PetscScalar);
1008   /* the next line frees the default space generated by the Create() */
1009   ierr = PetscFree(b->a);CHKERRQ(ierr);
1010   ierr = PetscFree(b->ilen);CHKERRQ(ierr);
1011   ierr = PetscMalloc(len+1,&b->a);CHKERRQ(ierr);
1012   b->j          = ajnew;
1013   b->i          = ainew;
1014   for (i=0; i<n; i++) dloc[i] += ainew[i];
1015   b->diag       = dloc;
1016   b->ilen       = 0;
1017   b->imax       = 0;
1018   b->row        = isrow;
1019   b->col        = iscol;
1020   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
1021   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
1022   b->icol       = isicol;
1023   ierr = PetscMalloc((n+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
1024   /* In b structure:  Free imax, ilen, old a, old j.
1025      Allocate dloc, solve_work, new a, new j */
1026   PetscLogObjectMemory(*fact,(ainew[n]+shift-n) * (sizeof(int)+sizeof(PetscScalar)));
1027   b->maxnz          = b->nz = ainew[n] + shift;
1028   if (info) {
1029     b->lu_damp      = (PetscTruth)((int)info->damp);
1030     b->lu_damping   = info->damping;
1031     b->lu_zeropivot = info->zeropivot;
1032   } else {
1033     b->lu_damp      = PETSC_FALSE;
1034     b->lu_damping   = 0.0;
1035     b->lu_zeropivot = 1.e-12;
1036   }
1037   (*fact)->factor   = FACTOR_LU;
1038   ierr = Mat_AIJ_CheckInode(*fact,PETSC_FALSE);CHKERRQ(ierr);
1039   (*fact)->ops->lufactornumeric   =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
1040 
1041   (*fact)->info.factor_mallocs    = realloc;
1042   (*fact)->info.fill_ratio_given  = f;
1043   (*fact)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[prow]);
1044   (*fact)->factor                 =  FACTOR_LU;
1045   PetscFunctionReturn(0);
1046 }
1047 
1048 
1049 
1050 
1051