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