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