xref: /petsc/src/mat/impls/aij/seq/aijfact.c (revision 35d8aa7f675fa4ee5c80aa031e65d3ef307a4891) 
1 /*$Id: aijfact.c,v 1.156 2000/09/28 21:11:00 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 __FUNC__
8 #define __FUNC__ /*<a name=""></a>*/"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 __FUNC__
28 #define __FUNC__ /*<a name=""></a>*/"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   new_i = (int*)   PetscMalloc((n+1)*sizeof(int));CHKPTRQ(new_i);
94   new_j = (int*)   PetscMalloc(jmax*sizeof(int));CHKPTRQ(new_j);
95   new_a = (Scalar*)PetscMalloc(jmax*sizeof(Scalar));CHKPTRQ(new_a);
96 
97   ordcol = (int*)PetscMalloc(n*sizeof(int));CHKPTRQ(ordcol);
98 
99   /* ------------------------------------------------------------
100      Make sure that everything is Fortran formatted (1-Based)
101      ------------------------------------------------------------
102   */
103   if (ishift) {
104     for (i=old_i[0];i<old_i[n];i++) {
105       old_j[i]++;
106     }
107     for(i=0;i<n+1;i++) {
108       old_i[i]++;
109     };
110   };
111 
112   if (reorder) {
113     ierr = ISGetIndices(iscol,&c);           CHKERRQ(ierr);
114     ierr = ISGetIndices(isrow,&r);           CHKERRQ(ierr);
115     for(i=0;i<n;i++) {
116       r[i]  = r[i]+1;
117       c[i]  = c[i]+1;
118     }
119     old_i2 = (int*)PetscMalloc((n+1)*sizeof(int));CHKPTRQ(old_i2);
120     old_j2 = (int*)PetscMalloc((old_i[n]-old_i[0]+1)*sizeof(int));CHKPTRQ(old_j2);
121     old_a2 = (Scalar*)PetscMalloc((old_i[n]-old_i[0]+1)*sizeof(Scalar));CHKPTRQ(old_a2);
122     job = 3; SPARSEKIT2dperm(&n,old_a,old_j,old_i,old_a2,old_j2,old_i2,r,c,&job);
123     for (i=0;i<n;i++) {
124       r[i]  = r[i]-1;
125       c[i]  = c[i]-1;
126     }
127     ierr = ISRestoreIndices(iscol,&c);CHKERRQ(ierr);
128     ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
129     o_a = old_a2;
130     o_j = old_j2;
131     o_i = old_i2;
132   } else {
133     o_a = old_a;
134     o_j = old_j;
135     o_i = old_i;
136   }
137 
138   /* ------------------------------------------------------------
139      Call Saad's ilutp() routine to generate the factorization
140      ------------------------------------------------------------
141   */
142 
143   iperm   = (int*)   PetscMalloc(2*n*sizeof(int));CHKPTRQ(iperm);
144   jw      = (int*)   PetscMalloc(2*n*sizeof(int));CHKPTRQ(jw);
145   w       = (Scalar*)PetscMalloc(n*sizeof(Scalar));CHKPTRQ(w);
146 
147   SPARSEKIT2ilutp(&n,o_a,o_j,o_i,&lfill,&info->dt,&permtol,&n,new_a,new_j,new_i,&jmax,w,jw,iperm,&ierr);
148   if (ierr) {
149     switch (ierr) {
150       case -3: SETERRQ2(1,"ilutp(), matrix U overflows, need larger info->fill current fill %g space allocated %d",info->fill,jmax);
151       case -2: SETERRQ2(1,"ilutp(), matrix L overflows, need larger info->fill current fill %g space allocated %d",info->fill,jmax);
152       case -5: SETERRQ(1,"ilutp(), zero row encountered");
153       case -1: SETERRQ(1,"ilutp(), input matrix may be wrong");
154       case -4: SETERRQ1(1,"ilutp(), illegal info->fill value %d",jmax);
155       default: SETERRQ1(1,"ilutp(), zero pivot detected on row %d",ierr);
156     }
157   }
158 
159   ierr = PetscFree(w);CHKERRQ(ierr);
160   ierr = PetscFree(jw);CHKERRQ(ierr);
161 
162   /* ------------------------------------------------------------
163      Saad's routine gives the result in Modified Sparse Row (msr)
164      Convert to Compressed Sparse Row format (csr)
165      ------------------------------------------------------------
166   */
167 
168   wk  = (Scalar*)PetscMalloc(n*sizeof(Scalar));CHKPTRQ(wk);
169   iwk = (int*)   PetscMalloc((n+1)*sizeof(int));CHKPTRQ(iwk);
170 
171   SPARSEKIT2msrcsr(&n,new_a,new_j,new_a,new_j,new_i,wk,iwk);
172 
173   ierr = PetscFree(iwk);CHKERRQ(ierr);
174   ierr = PetscFree(wk);CHKERRQ(ierr);
175 
176   if (reorder) {
177     ierr = PetscFree(old_a2);CHKERRQ(ierr);
178     ierr = PetscFree(old_j2);CHKERRQ(ierr);
179     ierr = PetscFree(old_i2);CHKERRQ(ierr);
180   } else {
181     /* fix permutation of old_j that the factorization introduced */
182     for (i=old_i[0]; i<old_i[n]; i++) {
183       old_j[i-1] = iperm[old_j[i-1]-1];
184     }
185   }
186 
187   /* get rid of the shift to indices starting at 1 */
188   if (ishift) {
189     for (i=0; i<n+1; i++) {
190       old_i[i]--;
191     }
192     for (i=old_i[0];i<old_i[n];i++) {
193       old_j[i]--;
194     }
195   }
196 
197   /* Make the factored matrix 0-based */
198   if (ishift) {
199     for (i=0; i<n+1; i++) {
200       new_i[i]--;
201     }
202     for (i=new_i[0];i<new_i[n];i++) {
203       new_j[i]--;
204     }
205   }
206 
207   /*-- due to the pivoting, we need to reorder iscol to correctly --*/
208   /*-- permute the right-hand-side and solution vectors           --*/
209   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
210   ierr = ISInvertPermutation(isrow,PETSC_DECIDE,&isirow);CHKERRQ(ierr);
211   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
212   for(i=0; i<n; i++) {
213     ordcol[i] = ic[iperm[i]-1];
214   };
215   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
216   ierr = ISDestroy(isicol);CHKERRQ(ierr);
217 
218   ierr = PetscFree(iperm);CHKERRQ(ierr);
219 
220   ierr = ISCreateGeneral(PETSC_COMM_SELF,n,ordcol,&iscolf);CHKERRQ(ierr);
221   ierr = PetscFree(ordcol);CHKERRQ(ierr);
222 
223   /*----- put together the new matrix -----*/
224 
225   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,fact);CHKERRQ(ierr);
226   (*fact)->factor    = FACTOR_LU;
227   (*fact)->assembled = PETSC_TRUE;
228 
229   b = (Mat_SeqAIJ*)(*fact)->data;
230   ierr = PetscFree(b->imax);CHKERRQ(ierr);
231   b->sorted        = PETSC_FALSE;
232   b->singlemalloc  = PETSC_FALSE;
233   /* the next line frees the default space generated by the MatCreate() */
234   ierr             = PetscFree(b->a);CHKERRQ(ierr);
235   ierr             = PetscFree(b->ilen);CHKERRQ(ierr);
236   b->a             = new_a;
237   b->j             = new_j;
238   b->i             = new_i;
239   b->ilen          = 0;
240   b->imax          = 0;
241   /*  I am not sure why these are the inverses of the row and column permutations; but the other way is NO GOOD */
242   b->row           = isirow;
243   b->col           = iscolf;
244   b->solve_work    =  (Scalar*)PetscMalloc((n+1)*sizeof(Scalar));CHKPTRQ(b->solve_work);
245   b->maxnz = b->nz = new_i[n];
246   b->indexshift    = a->indexshift;
247   ierr = MatMarkDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
248   (*fact)->info.factor_mallocs = 0;
249 
250   ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);
251 
252   /* check out for identical nodes. If found, use inode functions */
253   ierr = Mat_AIJ_CheckInode(*fact);CHKERRQ(ierr);
254 
255   af = ((double)b->nz)/((double)a->nz) + .001;
256   PLogInfo(A,"MatILUDTFactor_SeqAIJ:Fill ratio:given %g needed %g\n",info->fill,af);
257   PLogInfo(A,"MatILUDTFactor_SeqAIJ:Run with -pc_ilu_fill %g or use \n",af);
258   PLogInfo(A,"MatILUDTFactor_SeqAIJ:PCILUSetFill(pc,%g);\n",af);
259   PLogInfo(A,"MatILUDTFactor_SeqAIJ:for best performance.\n");
260 
261   PetscFunctionReturn(0);
262 }
263 
264 /*
265     Factorization code for AIJ format.
266 */
267 #undef __FUNC__
268 #define __FUNC__ /*<a name="MatLUFactorSymbolic_SeqAIJ""></a>*/"MatLUFactorSymbolic_SeqAIJ"
269 int MatLUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatLUInfo *info,Mat *B)
270 {
271   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b;
272   IS         isicol;
273   int        *r,*ic,ierr,i,n = A->m,*ai = a->i,*aj = a->j;
274   int        *ainew,*ajnew,jmax,*fill,*ajtmp,nz,shift = a->indexshift;
275   int        *idnew,idx,row,m,fm,nnz,nzi,realloc = 0,nzbd,*im;
276   PetscReal  f;
277 
278   PetscFunctionBegin;
279   PetscValidHeaderSpecific(isrow,IS_COOKIE);
280   PetscValidHeaderSpecific(iscol,IS_COOKIE);
281   if (A->M != A->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
282 
283   if (!isrow) {
284     ierr = ISCreateStride(PETSC_COMM_SELF,A->M,0,1,&isrow);CHKERRQ(ierr);
285   }
286   if (!iscol) {
287     ierr = ISCreateStride(PETSC_COMM_SELF,A->M,0,1,&iscol);CHKERRQ(ierr);
288   }
289 
290   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
291   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
292   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
293 
294   /* get new row pointers */
295   ainew    = (int*)PetscMalloc((n+1)*sizeof(int));CHKPTRQ(ainew);
296   ainew[0] = -shift;
297   /* don't know how many column pointers are needed so estimate */
298   if (info) f = info->fill; else f = 1.0;
299   jmax  = (int)(f*ai[n]+(!shift));
300   ajnew = (int*)PetscMalloc((jmax)*sizeof(int));CHKPTRQ(ajnew);
301   /* fill is a linked list of nonzeros in active row */
302   fill = (int*)PetscMalloc((2*n+1)*sizeof(int));CHKPTRQ(fill);
303   im   = fill + n + 1;
304   /* idnew is location of diagonal in factor */
305   idnew    = (int*)PetscMalloc((n+1)*sizeof(int));CHKPTRQ(idnew);
306   idnew[0] = -shift;
307 
308   for (i=0; i<n; i++) {
309     /* first copy previous fill into linked list */
310     nnz     = nz    = ai[r[i]+1] - ai[r[i]];
311     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
312     ajtmp   = aj + ai[r[i]] + shift;
313     fill[n] = n;
314     while (nz--) {
315       fm  = n;
316       idx = ic[*ajtmp++ + shift];
317       do {
318         m  = fm;
319         fm = fill[m];
320       } while (fm < idx);
321       fill[m]   = idx;
322       fill[idx] = fm;
323     }
324     row = fill[n];
325     while (row < i) {
326       ajtmp = ajnew + idnew[row] + (!shift);
327       nzbd  = 1 + idnew[row] - ainew[row];
328       nz    = im[row] - nzbd;
329       fm    = row;
330       while (nz-- > 0) {
331         idx = *ajtmp++ + shift;
332         nzbd++;
333         if (idx == i) im[row] = nzbd;
334         do {
335           m  = fm;
336           fm = fill[m];
337         } while (fm < idx);
338         if (fm != idx) {
339           fill[m]   = idx;
340           fill[idx] = fm;
341           fm        = idx;
342           nnz++;
343         }
344       }
345       row = fill[row];
346     }
347     /* copy new filled row into permanent storage */
348     ainew[i+1] = ainew[i] + nnz;
349     if (ainew[i+1] > jmax) {
350 
351       /* estimate how much additional space we will need */
352       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
353       /* just double the memory each time */
354       int maxadd = jmax;
355       /* maxadd = (int)((f*(ai[n]+(!shift))*(n-i+5))/n); */
356       if (maxadd < nnz) maxadd = (n-i)*(nnz+1);
357       jmax += maxadd;
358 
359       /* allocate a longer ajnew */
360       ajtmp = (int*)PetscMalloc(jmax*sizeof(int));CHKPTRQ(ajtmp);
361       ierr  = PetscMemcpy(ajtmp,ajnew,(ainew[i]+shift)*sizeof(int));CHKERRQ(ierr);
362       ierr  = PetscFree(ajnew);CHKERRQ(ierr);
363       ajnew = ajtmp;
364       realloc++; /* count how many times we realloc */
365     }
366     ajtmp = ajnew + ainew[i] + shift;
367     fm    = fill[n];
368     nzi   = 0;
369     im[i] = nnz;
370     while (nnz--) {
371       if (fm < i) nzi++;
372       *ajtmp++ = fm - shift;
373       fm       = fill[fm];
374     }
375     idnew[i] = ainew[i] + nzi;
376   }
377   if (ai[n] != 0) {
378     PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
379     PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af);
380     PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Run with -pc_lu_fill %g or use \n",af);
381     PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:PCLUSetFill(pc,%g);\n",af);
382     PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:for best performance.\n");
383   } else {
384     PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ: Empty matrix\n");
385   }
386 
387   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
388   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
389 
390   ierr = PetscFree(fill);CHKERRQ(ierr);
391 
392   /* put together the new matrix */
393   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,B);CHKERRQ(ierr);
394   PLogObjectParent(*B,isicol);
395   b = (Mat_SeqAIJ*)(*B)->data;
396   ierr = PetscFree(b->imax);CHKERRQ(ierr);
397   b->singlemalloc = PETSC_FALSE;
398   /* the next line frees the default space generated by the Create() */
399   ierr = PetscFree(b->a);CHKERRQ(ierr);
400   ierr = PetscFree(b->ilen);CHKERRQ(ierr);
401   b->a          = (Scalar*)PetscMalloc((ainew[n]+shift+1)*sizeof(Scalar));CHKPTRQ(b->a);
402   b->j          = ajnew;
403   b->i          = ainew;
404   b->diag       = idnew;
405   b->ilen       = 0;
406   b->imax       = 0;
407   b->row        = isrow;
408   b->col        = iscol;
409   if (info) {
410     b->lu_damp    = (PetscTruth) ((int)info->damp);
411     b->lu_damping = info->damping;
412   } else {
413     b->lu_damp    = PETSC_FALSE;
414     b->lu_damping = 0.0;
415   }
416   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
417   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
418   b->icol       = isicol;
419   b->solve_work = (Scalar*)PetscMalloc((n+1)*sizeof(Scalar));CHKPTRQ(b->solve_work);
420   /* In b structure:  Free imax, ilen, old a, old j.
421      Allocate idnew, solve_work, new a, new j */
422   PLogObjectMemory(*B,(ainew[n]+shift-n)*(sizeof(int)+sizeof(Scalar)));
423   b->maxnz = b->nz = ainew[n] + shift;
424 
425   (*B)->factor                 =  FACTOR_LU;
426   (*B)->info.factor_mallocs    = realloc;
427   (*B)->info.fill_ratio_given  = f;
428   (*B)->ops->lufactornumeric   =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
429 
430   if (ai[n] != 0) {
431     (*B)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
432   } else {
433     (*B)->info.fill_ratio_needed = 0.0;
434   }
435   PetscFunctionReturn(0);
436 }
437 /* ----------------------------------------------------------- */
438 EXTERN int Mat_AIJ_CheckInode(Mat);
439 
440 #undef __FUNC__
441 #define __FUNC__ /*<a name="MatLUFactorNumeric_SeqAIJ"></a>*/"MatLUFactorNumeric_SeqAIJ"
442 int MatLUFactorNumeric_SeqAIJ(Mat A,Mat *B)
443 {
444   Mat        C = *B;
445   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ *)C->data;
446   IS         isrow = b->row,isicol = b->icol;
447   int        *r,*ic,ierr,i,j,n = A->m,*ai = b->i,*aj = b->j;
448   int        *ajtmpold,*ajtmp,nz,row,*ics,shift = a->indexshift;
449   int        *diag_offset = b->diag,diag;
450   int        *pj,ndamp = 0;
451   Scalar     *rtmp,*v,*pc,multiplier,*pv,*rtmps;
452   PetscReal  damping = b->lu_damping;
453   PetscTruth damp;
454 
455   PetscFunctionBegin;
456 
457   ierr  = ISGetIndices(isrow,&r);CHKERRQ(ierr);
458   ierr  = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
459   rtmp  = (Scalar*)PetscMalloc((n+1)*sizeof(Scalar));CHKPTRQ(rtmp);
460   ierr  = PetscMemzero(rtmp,(n+1)*sizeof(Scalar));CHKERRQ(ierr);
461   rtmps = rtmp + shift; ics = ic + shift;
462 
463   do {
464     damp = PETSC_FALSE;
465     for (i=0; i<n; i++) {
466       nz    = ai[i+1] - ai[i];
467       ajtmp = aj + ai[i] + shift;
468       for  (j=0; j<nz; j++) rtmps[ajtmp[j]] = 0.0;
469 
470       /* load in initial (unfactored row) */
471       nz       = a->i[r[i]+1] - a->i[r[i]];
472       ajtmpold = a->j + a->i[r[i]] + shift;
473       v        = a->a + a->i[r[i]] + shift;
474       for (j=0; j<nz; j++) {
475         rtmp[ics[ajtmpold[j]]] = v[j];
476         if (ajtmpold[j] == r[i]) { /* damp the diagonal of the matrix */
477           rtmp[ics[ajtmpold[j]]] += damping;
478         }
479       }
480 
481       row = *ajtmp++ + shift;
482       while  (row < i) {
483         pc = rtmp + row;
484         if (*pc != 0.0) {
485           pv         = b->a + diag_offset[row] + shift;
486           pj         = b->j + diag_offset[row] + (!shift);
487           multiplier = *pc / *pv++;
488           *pc        = multiplier;
489           nz         = ai[row+1] - diag_offset[row] - 1;
490           for (j=0; j<nz; j++) rtmps[pj[j]] -= multiplier * pv[j];
491           PLogFlops(2*nz);
492         }
493         row = *ajtmp++ + shift;
494       }
495       /* finished row so stick it into b->a */
496       pv = b->a + ai[i] + shift;
497       pj = b->j + ai[i] + shift;
498       nz = ai[i+1] - ai[i];
499       for (j=0; j<nz; j++) {pv[j] = rtmps[pj[j]];}
500       diag = diag_offset[i] - ai[i];
501       if (PetscAbsScalar(pv[diag]) < 1.e-12) {
502         if (b->lu_damp) {
503           damp = PETSC_TRUE;
504           if (damping) damping *= 2.0;
505           else damping = 1.e-12;
506           ndamp++;
507           break;
508         } else {
509           SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %d",i);
510         }
511       }
512     }
513   } while (damp);
514 
515   /* invert diagonal entries for simplier triangular solves */
516   for (i=0; i<n; i++) {
517     b->a[diag_offset[i]+shift] = 1.0/b->a[diag_offset[i]+shift];
518   }
519 
520   ierr = PetscFree(rtmp);CHKERRQ(ierr);
521   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
522   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
523   C->factor = FACTOR_LU;
524   ierr = Mat_AIJ_CheckInode(C);CHKERRQ(ierr);
525   (*B)->ops->lufactornumeric   =  A->ops->lufactornumeric; /* Use Inode variant ONLY if A has inodes */
526   C->assembled = PETSC_TRUE;
527   PLogFlops(C->n);
528   if (ndamp || b->lu_damping) {
529     PLogInfo(0,"MatLUFactorNumerical_SeqAIJ: number of damping tries %d damping value %g\n",ndamp,damping);
530   }
531   PetscFunctionReturn(0);
532 }
533 /* ----------------------------------------------------------- */
534 #undef __FUNC__
535 #define __FUNC__ /*<a name="MatLUFactor_SeqAIJ"></a>*/"MatLUFactor_SeqAIJ"
536 int MatLUFactor_SeqAIJ(Mat A,IS row,IS col,MatLUInfo *info)
537 {
538   int ierr;
539   Mat C;
540 
541   PetscFunctionBegin;
542   ierr = MatLUFactorSymbolic(A,row,col,info,&C);CHKERRQ(ierr);
543   ierr = MatLUFactorNumeric(A,&C);CHKERRQ(ierr);
544   ierr = MatHeaderCopy(A,C);CHKERRQ(ierr);
545   PetscFunctionReturn(0);
546 }
547 /* ----------------------------------------------------------- */
548 #undef __FUNC__
549 #define __FUNC__ /*<a name=""></a>*/"MatSolve_SeqAIJ"
550 int MatSolve_SeqAIJ(Mat A,Vec bb,Vec xx)
551 {
552   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
553   IS         iscol = a->col,isrow = a->row;
554   int        *r,*c,ierr,i, n = A->m,*vi,*ai = a->i,*aj = a->j;
555   int        nz,shift = a->indexshift,*rout,*cout;
556   Scalar     *x,*b,*tmp,*tmps,*aa = a->a,sum,*v;
557 
558   PetscFunctionBegin;
559   if (!n) PetscFunctionReturn(0);
560 
561   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
562   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
563   tmp  = a->solve_work;
564 
565   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
566   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
567 
568   /* forward solve the lower triangular */
569   tmp[0] = b[*r++];
570   tmps   = tmp + shift;
571   for (i=1; i<n; i++) {
572     v   = aa + ai[i] + shift;
573     vi  = aj + ai[i] + shift;
574     nz  = a->diag[i] - ai[i];
575     sum = b[*r++];
576     while (nz--) sum -= *v++ * tmps[*vi++];
577     tmp[i] = sum;
578   }
579 
580   /* backward solve the upper triangular */
581   for (i=n-1; i>=0; i--){
582     v   = aa + a->diag[i] + (!shift);
583     vi  = aj + a->diag[i] + (!shift);
584     nz  = ai[i+1] - a->diag[i] - 1;
585     sum = tmp[i];
586     while (nz--) sum -= *v++ * tmps[*vi++];
587     x[*c--] = tmp[i] = sum*aa[a->diag[i]+shift];
588   }
589 
590   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
591   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
592   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
593   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
594   PLogFlops(2*a->nz - A->n);
595   PetscFunctionReturn(0);
596 }
597 
598 /* ----------------------------------------------------------- */
599 #undef __FUNC__
600 #define __FUNC__ /*<a name=""></a>*/"MatSolve_SeqAIJ_NaturalOrdering"
601 int MatSolve_SeqAIJ_NaturalOrdering(Mat A,Vec bb,Vec xx)
602 {
603   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
604   int        n = A->m,*ai = a->i,*aj = a->j,*adiag = a->diag,ierr;
605   Scalar     *x,*b,*aa = a->a,sum;
606 #if !defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
607   int        adiag_i,i,*vi,nz,ai_i;
608   Scalar     *v;
609 #endif
610 
611   PetscFunctionBegin;
612   if (!n) PetscFunctionReturn(0);
613   if (a->indexshift) {
614      ierr = MatSolve_SeqAIJ(A,bb,xx);CHKERRQ(ierr);
615      PetscFunctionReturn(0);
616   }
617 
618   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
619   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
620 
621 #if defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
622   fortransolveaij_(&n,x,ai,aj,adiag,aa,b);
623 #else
624   /* forward solve the lower triangular */
625   x[0] = b[0];
626   for (i=1; i<n; i++) {
627     ai_i = ai[i];
628     v    = aa + ai_i;
629     vi   = aj + ai_i;
630     nz   = adiag[i] - ai_i;
631     sum  = b[i];
632     while (nz--) sum -= *v++ * x[*vi++];
633     x[i] = sum;
634   }
635 
636   /* backward solve the upper triangular */
637   for (i=n-1; i>=0; i--){
638     adiag_i = adiag[i];
639     v       = aa + adiag_i + 1;
640     vi      = aj + adiag_i + 1;
641     nz      = ai[i+1] - adiag_i - 1;
642     sum     = x[i];
643     while (nz--) sum -= *v++ * x[*vi++];
644     x[i]    = sum*aa[adiag_i];
645   }
646 #endif
647   PLogFlops(2*a->nz - A->n);
648   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
649   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
650   PetscFunctionReturn(0);
651 }
652 
653 #undef __FUNC__
654 #define __FUNC__ /*<a name=""></a>*/"MatSolveAdd_SeqAIJ"
655 int MatSolveAdd_SeqAIJ(Mat A,Vec bb,Vec yy,Vec xx)
656 {
657   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
658   IS         iscol = a->col,isrow = a->row;
659   int        *r,*c,ierr,i, n = A->m,*vi,*ai = a->i,*aj = a->j;
660   int        nz,shift = a->indexshift,*rout,*cout;
661   Scalar     *x,*b,*tmp,*aa = a->a,sum,*v;
662 
663   PetscFunctionBegin;
664   if (yy != xx) {ierr = VecCopy(yy,xx);CHKERRQ(ierr);}
665 
666   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
667   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
668   tmp  = a->solve_work;
669 
670   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
671   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
672 
673   /* forward solve the lower triangular */
674   tmp[0] = b[*r++];
675   for (i=1; i<n; i++) {
676     v   = aa + ai[i] + shift;
677     vi  = aj + ai[i] + shift;
678     nz  = a->diag[i] - ai[i];
679     sum = b[*r++];
680     while (nz--) sum -= *v++ * tmp[*vi++ + shift];
681     tmp[i] = sum;
682   }
683 
684   /* backward solve the upper triangular */
685   for (i=n-1; i>=0; i--){
686     v   = aa + a->diag[i] + (!shift);
687     vi  = aj + a->diag[i] + (!shift);
688     nz  = ai[i+1] - a->diag[i] - 1;
689     sum = tmp[i];
690     while (nz--) sum -= *v++ * tmp[*vi++ + shift];
691     tmp[i] = sum*aa[a->diag[i]+shift];
692     x[*c--] += tmp[i];
693   }
694 
695   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
696   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
697   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
698   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
699   PLogFlops(2*a->nz);
700 
701   PetscFunctionReturn(0);
702 }
703 /* -------------------------------------------------------------------*/
704 #undef __FUNC__
705 #define __FUNC__ /*<a name=""></a>*/"MatSolveTranspose_SeqAIJ"
706 int MatSolveTranspose_SeqAIJ(Mat A,Vec bb,Vec xx)
707 {
708   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
709   IS         iscol = a->col,isrow = a->row;
710   int        *r,*c,ierr,i,n = A->m,*vi,*ai = a->i,*aj = a->j;
711   int        nz,shift = a->indexshift,*rout,*cout,*diag = a->diag;
712   Scalar     *x,*b,*tmp,*aa = a->a,*v,s1;
713 
714   PetscFunctionBegin;
715   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
716   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
717   tmp  = a->solve_work;
718 
719   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
720   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
721 
722   /* copy the b into temp work space according to permutation */
723   for (i=0; i<n; i++) tmp[i] = b[c[i]];
724 
725   /* forward solve the U^T */
726   for (i=0; i<n; i++) {
727     v   = aa + diag[i] + shift;
728     vi  = aj + diag[i] + (!shift);
729     nz  = ai[i+1] - diag[i] - 1;
730     s1  = tmp[i];
731     s1 *= *(v++);  /* multiply by inverse of diagonal entry */
732     while (nz--) {
733       tmp[*vi++ + shift] -= (*v++)*s1;
734     }
735     tmp[i] = s1;
736   }
737 
738   /* backward solve the L^T */
739   for (i=n-1; i>=0; i--){
740     v   = aa + diag[i] - 1 + shift;
741     vi  = aj + diag[i] - 1 + shift;
742     nz  = diag[i] - ai[i];
743     s1  = tmp[i];
744     while (nz--) {
745       tmp[*vi-- + shift] -= (*v--)*s1;
746     }
747   }
748 
749   /* copy tmp into x according to permutation */
750   for (i=0; i<n; i++) x[r[i]] = tmp[i];
751 
752   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
753   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
754   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
755   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
756 
757   PLogFlops(2*a->nz-A->n);
758   PetscFunctionReturn(0);
759 }
760 
761 #undef __FUNC__
762 #define __FUNC__ /*<a name=""></a>*/"MatSolveTransposeAdd_SeqAIJ"
763 int MatSolveTransposeAdd_SeqAIJ(Mat A,Vec bb,Vec zz,Vec xx)
764 {
765   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
766   IS         iscol = a->col,isrow = a->row;
767   int        *r,*c,ierr,i,n = A->m,*vi,*ai = a->i,*aj = a->j;
768   int        nz,shift = a->indexshift,*rout,*cout,*diag = a->diag;
769   Scalar     *x,*b,*tmp,*aa = a->a,*v;
770 
771   PetscFunctionBegin;
772   if (zz != xx) VecCopy(zz,xx);
773 
774   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
775   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
776   tmp = a->solve_work;
777 
778   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
779   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
780 
781   /* copy the b into temp work space according to permutation */
782   for (i=0; i<n; i++) tmp[i] = b[c[i]];
783 
784   /* forward solve the U^T */
785   for (i=0; i<n; i++) {
786     v   = aa + diag[i] + shift;
787     vi  = aj + diag[i] + (!shift);
788     nz  = ai[i+1] - diag[i] - 1;
789     tmp[i] *= *v++;
790     while (nz--) {
791       tmp[*vi++ + shift] -= (*v++)*tmp[i];
792     }
793   }
794 
795   /* backward solve the L^T */
796   for (i=n-1; i>=0; i--){
797     v   = aa + diag[i] - 1 + shift;
798     vi  = aj + diag[i] - 1 + shift;
799     nz  = diag[i] - ai[i];
800     while (nz--) {
801       tmp[*vi-- + shift] -= (*v--)*tmp[i];
802     }
803   }
804 
805   /* copy tmp into x according to permutation */
806   for (i=0; i<n; i++) x[r[i]] += tmp[i];
807 
808   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
809   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
810   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
811   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
812 
813   PLogFlops(2*a->nz);
814   PetscFunctionReturn(0);
815 }
816 /* ----------------------------------------------------------------*/
817 EXTERN int MatMissingDiagonal_SeqAIJ(Mat);
818 
819 #undef __FUNC__
820 #define __FUNC__ /*<a name=""></a>*/"MatILUFactorSymbolic_SeqAIJ"
821 int MatILUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatILUInfo *info,Mat *fact)
822 {
823   Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b;
824   IS         isicol;
825   int        *r,*ic,ierr,prow,n = A->m,*ai = a->i,*aj = a->j;
826   int        *ainew,*ajnew,jmax,*fill,*xi,nz,*im,*ajfill,*flev;
827   int        *dloc,idx,row,m,fm,nzf,nzi,len, realloc = 0,dcount = 0;
828   int        incrlev,nnz,i,shift = a->indexshift,levels,diagonal_fill;
829   PetscTruth col_identity,row_identity;
830   PetscReal  f;
831 
832   PetscFunctionBegin;
833   if (info) {
834     f             = info->fill;
835     levels        = (int)info->levels;
836     diagonal_fill = (int)info->diagonal_fill;
837   } else {
838     f             = 1.0;
839     levels        = 0;
840     diagonal_fill = 0;
841   }
842   if (!isrow) {
843     ierr = ISCreateStride(PETSC_COMM_SELF,A->M,0,1,&isrow);CHKERRQ(ierr);
844   }
845   if (!iscol) {
846     ierr = ISCreateStride(PETSC_COMM_SELF,A->M,0,1,&iscol);CHKERRQ(ierr);
847   }
848 
849   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
850 
851   /* special case that simply copies fill pattern */
852   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
853   ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
854   if (!levels && row_identity && col_identity) {
855     ierr = MatDuplicate_SeqAIJ(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr);
856     (*fact)->factor = FACTOR_LU;
857     b               = (Mat_SeqAIJ*)(*fact)->data;
858     if (!b->diag) {
859       ierr = MatMarkDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
860     }
861     ierr = MatMissingDiagonal_SeqAIJ(*fact);CHKERRQ(ierr);
862     b->row              = isrow;
863     b->col              = iscol;
864     b->icol             = isicol;
865     if (info) {
866       b->lu_damp    = (PetscTruth)((int)info->damp);
867       b->lu_damping = info->damping;
868     } else {
869       b->lu_damp    = PETSC_FALSE;
870       b->lu_damping = 0.0;
871     }
872     b->solve_work       = (Scalar*)PetscMalloc(((*fact)->m+1)*sizeof(Scalar));CHKPTRQ(b->solve_work);
873     (*fact)->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
874     ierr                = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
875     ierr                = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
876     PetscFunctionReturn(0);
877   }
878 
879   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
880   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
881 
882   /* get new row pointers */
883   ainew = (int*)PetscMalloc((n+1)*sizeof(int));CHKPTRQ(ainew);
884   ainew[0] = -shift;
885   /* don't know how many column pointers are needed so estimate */
886   jmax = (int)(f*(ai[n]+!shift));
887   ajnew = (int*)PetscMalloc((jmax)*sizeof(int));CHKPTRQ(ajnew);
888   /* ajfill is level of fill for each fill entry */
889   ajfill = (int*)PetscMalloc((jmax)*sizeof(int));CHKPTRQ(ajfill);
890   /* fill is a linked list of nonzeros in active row */
891   fill = (int*)PetscMalloc((n+1)*sizeof(int));CHKPTRQ(fill);
892   /* im is level for each filled value */
893   im = (int*)PetscMalloc((n+1)*sizeof(int));CHKPTRQ(im);
894   /* dloc is location of diagonal in factor */
895   dloc = (int*)PetscMalloc((n+1)*sizeof(int));CHKPTRQ(dloc);
896   dloc[0]  = 0;
897   for (prow=0; prow<n; prow++) {
898 
899     /* copy current row into linked list */
900     nzf     = nz  = ai[r[prow]+1] - ai[r[prow]];
901     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
902     xi      = aj + ai[r[prow]] + shift;
903     fill[n]    = n;
904     fill[prow] = -1; /* marker to indicate if diagonal exists */
905     while (nz--) {
906       fm  = n;
907       idx = ic[*xi++ + shift];
908       do {
909         m  = fm;
910         fm = fill[m];
911       } while (fm < idx);
912       fill[m]   = idx;
913       fill[idx] = fm;
914       im[idx]   = 0;
915     }
916 
917     /* make sure diagonal entry is included */
918     if (diagonal_fill && fill[prow] == -1) {
919       fm = n;
920       while (fill[fm] < prow) fm = fill[fm];
921       fill[prow] = fill[fm]; /* insert diagonal into linked list */
922       fill[fm]   = prow;
923       im[prow]   = 0;
924       nzf++;
925       dcount++;
926     }
927 
928     nzi = 0;
929     row = fill[n];
930     while (row < prow) {
931       incrlev = im[row] + 1;
932       nz      = dloc[row];
933       xi      = ajnew  + ainew[row] + shift + nz + 1;
934       flev    = ajfill + ainew[row] + shift + nz + 1;
935       nnz     = ainew[row+1] - ainew[row] - nz - 1;
936       fm      = row;
937       while (nnz-- > 0) {
938         idx = *xi++ + shift;
939         if (*flev + incrlev > levels) {
940           flev++;
941           continue;
942         }
943         do {
944           m  = fm;
945           fm = fill[m];
946         } while (fm < idx);
947         if (fm != idx) {
948           im[idx]   = *flev + incrlev;
949           fill[m]   = idx;
950           fill[idx] = fm;
951           fm        = idx;
952           nzf++;
953         } else {
954           if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev;
955         }
956         flev++;
957       }
958       row = fill[row];
959       nzi++;
960     }
961     /* copy new filled row into permanent storage */
962     ainew[prow+1] = ainew[prow] + nzf;
963     if (ainew[prow+1] > jmax-shift) {
964 
965       /* estimate how much additional space we will need */
966       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
967       /* just double the memory each time */
968       /*  maxadd = (int)((f*(ai[n]+!shift)*(n-prow+5))/n); */
969       int maxadd = jmax;
970       if (maxadd < nzf) maxadd = (n-prow)*(nzf+1);
971       jmax += maxadd;
972 
973       /* allocate a longer ajnew and ajfill */
974       xi     = (int*)PetscMalloc(jmax*sizeof(int));CHKPTRQ(xi);
975       ierr   = PetscMemcpy(xi,ajnew,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr);
976       ierr = PetscFree(ajnew);CHKERRQ(ierr);
977       ajnew  = xi;
978       xi     = (int*)PetscMalloc(jmax*sizeof(int));CHKPTRQ(xi);
979       ierr   = PetscMemcpy(xi,ajfill,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr);
980       ierr = PetscFree(ajfill);CHKERRQ(ierr);
981       ajfill = xi;
982       realloc++; /* count how many times we realloc */
983     }
984     xi          = ajnew + ainew[prow] + shift;
985     flev        = ajfill + ainew[prow] + shift;
986     dloc[prow]  = nzi;
987     fm          = fill[n];
988     while (nzf--) {
989       *xi++   = fm - shift;
990       *flev++ = im[fm];
991       fm      = fill[fm];
992     }
993     /* make sure row has diagonal entry */
994     if (ajnew[ainew[prow]+shift+dloc[prow]]+shift != prow) {
995       SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,"Row %d has missing diagonal in factored matrix\n\
996     try running with -pc_ilu_nonzeros_along_diagonal or -pc_ilu_diagonal_fill",prow);
997     }
998   }
999   ierr = PetscFree(ajfill);CHKERRQ(ierr);
1000   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
1001   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
1002   ierr = PetscFree(fill);CHKERRQ(ierr);
1003   ierr = PetscFree(im);CHKERRQ(ierr);
1004 
1005   {
1006     PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
1007     PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af);
1008     PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Run with -pc_ilu_fill %g or use \n",af);
1009     PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:PCILUSetFill(pc,%g);\n",af);
1010     PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:for best performance.\n");
1011     if (diagonal_fill) {
1012       PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Detected and replace %d missing diagonals",dcount);
1013     }
1014   }
1015 
1016   /* put together the new matrix */
1017   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,fact);CHKERRQ(ierr);
1018   PLogObjectParent(*fact,isicol);
1019   b = (Mat_SeqAIJ*)(*fact)->data;
1020   ierr = PetscFree(b->imax);CHKERRQ(ierr);
1021   b->singlemalloc = PETSC_FALSE;
1022   len = (ainew[n] + shift)*sizeof(Scalar);
1023   /* the next line frees the default space generated by the Create() */
1024   ierr = PetscFree(b->a);CHKERRQ(ierr);
1025   ierr = PetscFree(b->ilen);CHKERRQ(ierr);
1026   b->a          = (Scalar*)PetscMalloc(len+1);CHKPTRQ(b->a);
1027   b->j          = ajnew;
1028   b->i          = ainew;
1029   for (i=0; i<n; i++) dloc[i] += ainew[i];
1030   b->diag       = dloc;
1031   b->ilen       = 0;
1032   b->imax       = 0;
1033   b->row        = isrow;
1034   b->col        = iscol;
1035   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
1036   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
1037   b->icol       = isicol;
1038   b->solve_work = (Scalar*)PetscMalloc((n+1)*sizeof(Scalar));CHKPTRQ(b->solve_work);
1039   /* In b structure:  Free imax, ilen, old a, old j.
1040      Allocate dloc, solve_work, new a, new j */
1041   PLogObjectMemory(*fact,(ainew[n]+shift-n) * (sizeof(int)+sizeof(Scalar)));
1042   b->maxnz          = b->nz = ainew[n] + shift;
1043   if (info) {
1044     b->lu_damp    = (PetscTruth)((int)info->damp);
1045     b->lu_damping = info->damping;
1046   } else {
1047     b->lu_damp    = PETSC_FALSE;
1048     b->lu_damping = 0.0;
1049   }
1050   (*fact)->factor   = FACTOR_LU;
1051 
1052   (*fact)->info.factor_mallocs    = realloc;
1053   (*fact)->info.fill_ratio_given  = f;
1054   (*fact)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[prow]);
1055   (*fact)->factor                 =  FACTOR_LU;
1056 
1057   PetscFunctionReturn(0);
1058 }
1059 
1060 
1061 
1062 
1063