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