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