xref: /petsc/src/mat/impls/aij/seq/aijfact.c (revision c38d4ed214221df9ea04de46f7761bef149d00ff) 
1 /*$Id: aijfact.c,v 1.125 1999/11/05 14:45:18 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__ "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,0,"Code not written");
14 #if !defined(PETSC_USE_DEBUG)
15   PetscFunctionReturn(0);
16 #endif
17 }
18 
19 /*
20     Factorization code for AIJ format.
21 */
22 #undef __FUNC__
23 #define __FUNC__ "MatLUFactorSymbolic_SeqAIJ"
24 int MatLUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,double f,Mat *B)
25 {
26   Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data, *b;
27   IS         isicol;
28   int        *r,*ic, ierr, i, n = a->m, *ai = a->i, *aj = a->j;
29   int        *ainew,*ajnew, jmax,*fill, *ajtmp, nz,shift = a->indexshift;
30   int        *idnew, idx, row,m,fm, nnz, nzi, realloc = 0,nzbd,*im;
31 
32   PetscFunctionBegin;
33   PetscValidHeaderSpecific(isrow,IS_COOKIE);
34   PetscValidHeaderSpecific(iscol,IS_COOKIE);
35   if (A->M != A->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square");
36 
37   ierr = ISInvertPermutation(iscol,&isicol);CHKERRQ(ierr);
38   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
39   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
40 
41   /* get new row pointers */
42   ainew    = (int *) PetscMalloc( (n+1)*sizeof(int) );CHKPTRQ(ainew);
43   ainew[0] = -shift;
44   /* don't know how many column pointers are needed so estimate */
45   jmax  = (int) (f*ai[n]+(!shift));
46   ajnew = (int *) PetscMalloc( (jmax)*sizeof(int) );CHKPTRQ(ajnew);
47   /* fill is a linked list of nonzeros in active row */
48   fill = (int *) PetscMalloc( (2*n+1)*sizeof(int));CHKPTRQ(fill);
49   im   = fill + n + 1;
50   /* idnew is location of diagonal in factor */
51   idnew    = (int *) PetscMalloc( (n+1)*sizeof(int));CHKPTRQ(idnew);
52   idnew[0] = -shift;
53 
54   for ( i=0; i<n; i++ ) {
55     /* first copy previous fill into linked list */
56     nnz     = nz    = ai[r[i]+1] - ai[r[i]];
57     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,1,"Empty row in matrix");
58     ajtmp   = aj + ai[r[i]] + shift;
59     fill[n] = n;
60     while (nz--) {
61       fm  = n;
62       idx = ic[*ajtmp++ + shift];
63       do {
64         m  = fm;
65         fm = fill[m];
66       } while (fm < idx);
67       fill[m]   = idx;
68       fill[idx] = fm;
69     }
70     row = fill[n];
71     while ( row < i ) {
72       ajtmp = ajnew + idnew[row] + (!shift);
73       nzbd  = 1 + idnew[row] - ainew[row];
74       nz    = im[row] - nzbd;
75       fm    = row;
76       while (nz-- > 0) {
77         idx = *ajtmp++ + shift;
78         nzbd++;
79         if (idx == i) im[row] = nzbd;
80         do {
81           m  = fm;
82           fm = fill[m];
83         } while (fm < idx);
84         if (fm != idx) {
85           fill[m]   = idx;
86           fill[idx] = fm;
87           fm        = idx;
88           nnz++;
89         }
90       }
91       row = fill[row];
92     }
93     /* copy new filled row into permanent storage */
94     ainew[i+1] = ainew[i] + nnz;
95     if (ainew[i+1] > jmax) {
96 
97       /* estimate how much additional space we will need */
98       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
99       /* just double the memory each time */
100       int maxadd = jmax;
101       /* maxadd = (int) ((f*(ai[n]+(!shift))*(n-i+5))/n); */
102       if (maxadd < nnz) maxadd = (n-i)*(nnz+1);
103       jmax += maxadd;
104 
105       /* allocate a longer ajnew */
106       ajtmp = (int *) PetscMalloc( jmax*sizeof(int) );CHKPTRQ(ajtmp);
107       ierr  = PetscMemcpy(ajtmp,ajnew,(ainew[i]+shift)*sizeof(int));CHKERRQ(ierr);
108       ierr = PetscFree(ajnew);CHKERRQ(ierr);
109       ajnew = ajtmp;
110       realloc++; /* count how many times we realloc */
111     }
112     ajtmp = ajnew + ainew[i] + shift;
113     fm    = fill[n];
114     nzi   = 0;
115     im[i] = nnz;
116     while (nnz--) {
117       if (fm < i) nzi++;
118       *ajtmp++ = fm - shift;
119       fm       = fill[fm];
120     }
121     idnew[i] = ainew[i] + nzi;
122   }
123   if (ai[n] != 0) {
124     double af = ((double)ainew[n])/((double)ai[n]);
125     PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af);
126     PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Run with -pc_lu_fill %g or use \n",af);
127     PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:PCLUSetFill(pc,%g);\n",af);
128     PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:for best performance.\n");
129   } else {
130     PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ: Empty matrix\n");
131   }
132 
133   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
134   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
135 
136   ierr = PetscFree(fill);CHKERRQ(ierr);
137 
138   /* put together the new matrix */
139   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,B);CHKERRQ(ierr);
140   PLogObjectParent(*B,isicol);
141   b = (Mat_SeqAIJ *) (*B)->data;
142   ierr = PetscFree(b->imax);CHKERRQ(ierr);
143   b->singlemalloc = 0;
144   /* the next line frees the default space generated by the Create() */
145   ierr = PetscFree(b->a);CHKERRQ(ierr);
146   ierr = PetscFree(b->ilen);CHKERRQ(ierr);
147   b->a          = (Scalar *) PetscMalloc((ainew[n]+shift+1)*sizeof(Scalar));CHKPTRQ(b->a);
148   b->j          = ajnew;
149   b->i          = ainew;
150   b->diag       = idnew;
151   b->ilen       = 0;
152   b->imax       = 0;
153   b->row        = isrow;
154   b->col        = iscol;
155   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
156   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
157   b->icol       = isicol;
158   b->solve_work = (Scalar *) PetscMalloc( (n+1)*sizeof(Scalar));CHKPTRQ(b->solve_work);
159   /* In b structure:  Free imax, ilen, old a, old j.
160      Allocate idnew, solve_work, new a, new j */
161   PLogObjectMemory(*B,(ainew[n]+shift-n)*(sizeof(int)+sizeof(Scalar)));
162   b->maxnz = b->nz = ainew[n] + shift;
163 
164   (*B)->factor                 =  FACTOR_LU;;
165   (*B)->info.factor_mallocs    = realloc;
166   (*B)->info.fill_ratio_given  = f;
167   (*B)->ops->lufactornumeric   =  A->ops->lufactornumeric; /* Use Inode variant if A has inodes */
168 
169   if (ai[n] != 0) {
170     (*B)->info.fill_ratio_needed = ((double)ainew[n])/((double)ai[n]);
171   } else {
172     (*B)->info.fill_ratio_needed = 0.0;
173   }
174   PetscFunctionReturn(0);
175 }
176 /* ----------------------------------------------------------- */
177 extern int Mat_AIJ_CheckInode(Mat);
178 
179 #undef __FUNC__
180 #define __FUNC__ "MatLUFactorNumeric_SeqAIJ"
181 int MatLUFactorNumeric_SeqAIJ(Mat A,Mat *B)
182 {
183   Mat        C = *B;
184   Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data, *b = (Mat_SeqAIJ *)C->data;
185   IS         isrow = b->row, isicol = b->icol;
186   int        *r,*ic, ierr, i, j, n = a->m, *ai = b->i, *aj = b->j;
187   int        *ajtmpold, *ajtmp, nz, row, *ics, shift = a->indexshift;
188   int        *diag_offset = b->diag,diag,k;
189   int        preserve_row_sums = (int) a->ilu_preserve_row_sums;
190   register   int    *pj;
191   Scalar     *rtmp,*v, *pc, multiplier,sum,inner_sum,*rowsums = 0;
192   double     ssum;
193   register   Scalar *pv, *rtmps,*u_values;
194 
195   PetscFunctionBegin;
196 
197   ierr  = ISGetIndices(isrow,&r);CHKERRQ(ierr);
198   ierr  = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
199   rtmp  = (Scalar *) PetscMalloc( (n+1)*sizeof(Scalar) );CHKPTRQ(rtmp);
200   ierr  = PetscMemzero(rtmp,(n+1)*sizeof(Scalar));CHKERRQ(ierr);
201   rtmps = rtmp + shift; ics = ic + shift;
202 
203   /* precalculate row sums */
204   if (preserve_row_sums) {
205     rowsums = (Scalar *) PetscMalloc( n*sizeof(Scalar) );CHKPTRQ(rowsums);
206     for ( i=0; i<n; i++ ) {
207       nz  = a->i[r[i]+1] - a->i[r[i]];
208       v   = a->a + a->i[r[i]] + shift;
209       sum = 0.0;
210       for ( j=0; j<nz; j++ ) sum += v[j];
211       rowsums[i] = sum;
212     }
213   }
214 
215   for ( i=0; i<n; i++ ) {
216     nz    = ai[i+1] - ai[i];
217     ajtmp = aj + ai[i] + shift;
218     for  ( j=0; j<nz; j++ ) rtmps[ajtmp[j]] = 0.0;
219 
220     /* load in initial (unfactored row) */
221     nz       = a->i[r[i]+1] - a->i[r[i]];
222     ajtmpold = a->j + a->i[r[i]] + shift;
223     v        = a->a + a->i[r[i]] + shift;
224     for ( j=0; j<nz; j++ ) rtmp[ics[ajtmpold[j]]] =  v[j];
225 
226     row = *ajtmp++ + shift;
227     while  (row < i ) {
228       pc = rtmp + row;
229       if (*pc != 0.0) {
230         pv         = b->a + diag_offset[row] + shift;
231         pj         = b->j + diag_offset[row] + (!shift);
232         multiplier = *pc / *pv++;
233         *pc        = multiplier;
234         nz         = ai[row+1] - diag_offset[row] - 1;
235         for (j=0; j<nz; j++) rtmps[pj[j]] -= multiplier * pv[j];
236         PLogFlops(2*nz);
237       }
238       row = *ajtmp++ + shift;
239     }
240     /* finished row so stick it into b->a */
241     pv = b->a + ai[i] + shift;
242     pj = b->j + ai[i] + shift;
243     nz = ai[i+1] - ai[i];
244     for ( j=0; j<nz; j++ ) {pv[j] = rtmps[pj[j]];}
245     diag = diag_offset[i] - ai[i];
246     /*
247           Possibly adjust diagonal entry on current row to force
248         LU matrix to have same row sum as initial matrix.
249     */
250     if (pv[diag] == 0.0) {
251       SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,0,"Zero pivot row %d",i);
252     }
253     if (preserve_row_sums) {
254       pj  = b->j + ai[i] + shift;
255       sum = rowsums[i];
256       for ( j=0; j<diag; j++ ) {
257         u_values  = b->a + diag_offset[pj[j]] + shift;
258         nz        = ai[pj[j]+1] - diag_offset[pj[j]];
259         inner_sum = 0.0;
260         for ( k=0; k<nz; k++ ) {
261           inner_sum += u_values[k];
262         }
263         sum -= pv[j]*inner_sum;
264 
265       }
266       nz       = ai[i+1] - diag_offset[i] - 1;
267       u_values = b->a + diag_offset[i] + 1 + shift;
268       for ( k=0; k<nz; k++ ) {
269         sum -= u_values[k];
270       }
271       ssum = PetscAbsScalar(sum/pv[diag]);
272       if (ssum < 1000. && ssum > .001) pv[diag] = sum;
273     }
274     /* check pivot entry for current row */
275   }
276 
277   /* invert diagonal entries for simplier triangular solves */
278   for ( i=0; i<n; i++ ) {
279     b->a[diag_offset[i]+shift] = 1.0/b->a[diag_offset[i]+shift];
280   }
281 
282   if (preserve_row_sums) {ierr = PetscFree(rowsums);CHKERRQ(ierr);}
283   ierr = PetscFree(rtmp);CHKERRQ(ierr);
284   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
285   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
286   C->factor = FACTOR_LU;
287   ierr = Mat_AIJ_CheckInode(C);CHKERRQ(ierr);
288   C->assembled = PETSC_TRUE;
289   PLogFlops(b->n);
290   PetscFunctionReturn(0);
291 }
292 /* ----------------------------------------------------------- */
293 #undef __FUNC__
294 #define __FUNC__ "MatLUFactor_SeqAIJ"
295 int MatLUFactor_SeqAIJ(Mat A,IS row,IS col,double f)
296 {
297   Mat_SeqAIJ     *mat = (Mat_SeqAIJ *) A->data;
298   int            ierr;
299   Mat            C;
300   PetscOps       *Abops;
301   MatOps         Aops;
302 
303   PetscFunctionBegin;
304   ierr = MatLUFactorSymbolic(A,row,col,f,&C);CHKERRQ(ierr);
305   ierr = MatLUFactorNumeric(A,&C);CHKERRQ(ierr);
306 
307   /* free all the data structures from mat */
308   ierr = PetscFree(mat->a);CHKERRQ(ierr);
309   if (!mat->singlemalloc) {
310     ierr = PetscFree(mat->i);CHKERRQ(ierr);
311     ierr = PetscFree(mat->j);CHKERRQ(ierr);
312   }
313   if (mat->diag) {ierr = PetscFree(mat->diag);CHKERRQ(ierr);}
314   if (mat->ilen) {ierr = PetscFree(mat->ilen);CHKERRQ(ierr);}
315   if (mat->imax) {ierr = PetscFree(mat->imax);CHKERRQ(ierr);}
316   if (mat->solve_work) {ierr = PetscFree(mat->solve_work);CHKERRQ(ierr);}
317   if (mat->inode.size) {ierr = PetscFree(mat->inode.size);CHKERRQ(ierr);}
318   if (mat->icol) {ierr = ISDestroy(mat->icol);CHKERRQ(ierr);}
319   ierr = PetscFree(mat);CHKERRQ(ierr);
320 
321   ierr = MapDestroy(A->rmap);CHKERRQ(ierr);
322   ierr = MapDestroy(A->cmap);CHKERRQ(ierr);
323 
324   /*
325        This is horrible, horrible code. We need to keep the
326     A pointers for the bops and ops but copy everything
327     else from C.
328   */
329   Abops = A->bops;
330   Aops  = A->ops;
331   ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr);
332   mat  = (Mat_SeqAIJ *) A->data;
333   PLogObjectParent(A,mat->icol);
334 
335   A->bops  = Abops;
336   A->ops   = Aops;
337   A->qlist = 0;
338   /* copy over the type_name and name */
339   ierr     = PetscStrallocpy(C->type_name,&A->type_name);CHKERRQ(ierr);
340   ierr     = PetscStrallocpy(C->name,&A->name);CHKERRQ(ierr);
341 
342   PetscHeaderDestroy(C);
343 
344   PetscFunctionReturn(0);
345 }
346 /* ----------------------------------------------------------- */
347 #undef __FUNC__
348 #define __FUNC__ "MatSolve_SeqAIJ"
349 int MatSolve_SeqAIJ(Mat A,Vec bb, Vec xx)
350 {
351   Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data;
352   IS         iscol = a->col, isrow = a->row;
353   int        *r,*c, ierr, i,  n = a->m, *vi, *ai = a->i, *aj = a->j;
354   int        nz,shift = a->indexshift,*rout,*cout;
355   Scalar     *x,*b,*tmp, *tmps, *aa = a->a, sum, *v;
356 
357   PetscFunctionBegin;
358   if (!n) PetscFunctionReturn(0);
359 
360   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
361   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
362   tmp  = a->solve_work;
363 
364   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
365   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
366 
367   /* forward solve the lower triangular */
368   tmp[0] = b[*r++];
369   tmps   = tmp + shift;
370   for ( i=1; i<n; i++ ) {
371     v   = aa + ai[i] + shift;
372     vi  = aj + ai[i] + shift;
373     nz  = a->diag[i] - ai[i];
374     sum = b[*r++];
375     while (nz--) sum -= *v++ * tmps[*vi++];
376     tmp[i] = sum;
377   }
378 
379   /* backward solve the upper triangular */
380   for ( i=n-1; i>=0; i-- ){
381     v   = aa + a->diag[i] + (!shift);
382     vi  = aj + a->diag[i] + (!shift);
383     nz  = ai[i+1] - a->diag[i] - 1;
384     sum = tmp[i];
385     while (nz--) sum -= *v++ * tmps[*vi++];
386     x[*c--] = tmp[i] = sum*aa[a->diag[i]+shift];
387   }
388 
389   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
390   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
391   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
392   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
393   PLogFlops(2*a->nz - a->n);
394   PetscFunctionReturn(0);
395 }
396 
397 /* ----------------------------------------------------------- */
398 #undef __FUNC__
399 #define __FUNC__ "MatSolve_SeqAIJ_NaturalOrdering"
400 int MatSolve_SeqAIJ_NaturalOrdering(Mat A,Vec bb, Vec xx)
401 {
402   Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data;
403   int        n = a->m, *ai = a->i, *aj = a->j, *adiag = a->diag,ierr;
404   Scalar     *x,*b, *aa = a->a, sum;
405 #if !defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
406   int        adiag_i,i,*vi,nz,ai_i;
407   Scalar     *v;
408 #endif
409 
410   PetscFunctionBegin;
411   if (!n) PetscFunctionReturn(0);
412   if (a->indexshift) {
413      ierr = MatSolve_SeqAIJ(A,bb,xx);CHKERRQ(ierr);
414      PetscFunctionReturn(0);
415   }
416 
417   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
418   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
419 
420 #if defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ)
421   fortransolveaij_(&n,x,ai,aj,adiag,aa,b);
422 #else
423   /* forward solve the lower triangular */
424   x[0] = b[0];
425   for ( i=1; i<n; i++ ) {
426     ai_i = ai[i];
427     v    = aa + ai_i;
428     vi   = aj + ai_i;
429     nz   = adiag[i] - ai_i;
430     sum  = b[i];
431     while (nz--) sum -= *v++ * x[*vi++];
432     x[i] = sum;
433   }
434 
435   /* backward solve the upper triangular */
436   for ( i=n-1; i>=0; i-- ){
437     adiag_i = adiag[i];
438     v       = aa + adiag_i + 1;
439     vi      = aj + adiag_i + 1;
440     nz      = ai[i+1] - adiag_i - 1;
441     sum     = x[i];
442     while (nz--) sum -= *v++ * x[*vi++];
443     x[i]    = sum*aa[adiag_i];
444   }
445 #endif
446   PLogFlops(2*a->nz - a->n);
447   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
448   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
449   PetscFunctionReturn(0);
450 }
451 
452 #undef __FUNC__
453 #define __FUNC__ "MatSolveAdd_SeqAIJ"
454 int MatSolveAdd_SeqAIJ(Mat A,Vec bb, Vec yy, Vec xx)
455 {
456   Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data;
457   IS         iscol = a->col, isrow = a->row;
458   int        *r,*c, ierr, i,  n = a->m, *vi, *ai = a->i, *aj = a->j;
459   int        nz, shift = a->indexshift,*rout,*cout;
460   Scalar     *x,*b,*tmp, *aa = a->a, sum, *v;
461 
462   PetscFunctionBegin;
463   if (yy != xx) {ierr = VecCopy(yy,xx);CHKERRQ(ierr);}
464 
465   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
466   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
467   tmp  = a->solve_work;
468 
469   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
470   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1);
471 
472   /* forward solve the lower triangular */
473   tmp[0] = b[*r++];
474   for ( i=1; i<n; i++ ) {
475     v   = aa + ai[i] + shift;
476     vi  = aj + ai[i] + shift;
477     nz  = a->diag[i] - ai[i];
478     sum = b[*r++];
479     while (nz--) sum -= *v++ * tmp[*vi++ + shift];
480     tmp[i] = sum;
481   }
482 
483   /* backward solve the upper triangular */
484   for ( i=n-1; i>=0; i-- ){
485     v   = aa + a->diag[i] + (!shift);
486     vi  = aj + a->diag[i] + (!shift);
487     nz  = ai[i+1] - a->diag[i] - 1;
488     sum = tmp[i];
489     while (nz--) sum -= *v++ * tmp[*vi++ + shift];
490     tmp[i] = sum*aa[a->diag[i]+shift];
491     x[*c--] += tmp[i];
492   }
493 
494   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
495   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
496   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
497   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
498   PLogFlops(2*a->nz);
499 
500   PetscFunctionReturn(0);
501 }
502 /* -------------------------------------------------------------------*/
503 #undef __FUNC__
504 #define __FUNC__ "MatSolveTrans_SeqAIJ"
505 int MatSolveTrans_SeqAIJ(Mat A,Vec bb, Vec xx)
506 {
507   Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data;
508   IS         iscol = a->col, isrow = a->row;
509   int        *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j;
510   int        nz,shift = a->indexshift,*rout,*cout, *diag = a->diag;
511   Scalar     *x,*b,*tmp, *aa = a->a, *v, s1;
512 
513   PetscFunctionBegin;
514   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
515   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
516   tmp  = a->solve_work;
517 
518   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
519   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
520 
521   /* copy the b into temp work space according to permutation */
522   for ( i=0; i<n; i++ ) tmp[i] = b[c[i]];
523 
524   /* forward solve the U^T */
525   for ( i=0; i<n; i++ ) {
526     v   = aa + diag[i] + shift;
527     vi  = aj + diag[i] + (!shift);
528     nz  = ai[i+1] - diag[i] - 1;
529     s1  = tmp[i];
530     s1 *= *(v++);  /* multiply by inverse of diagonal entry */
531     while (nz--) {
532       tmp[*vi++ + shift] -= (*v++)*s1;
533     }
534     tmp[i] = s1;
535   }
536 
537   /* backward solve the L^T */
538   for ( i=n-1; i>=0; i-- ){
539     v   = aa + diag[i] - 1 + shift;
540     vi  = aj + diag[i] - 1 + shift;
541     nz  = diag[i] - ai[i];
542     s1  = tmp[i];
543     while (nz--) {
544       tmp[*vi-- + shift] -= (*v--)*s1;
545     }
546   }
547 
548   /* copy tmp into x according to permutation */
549   for ( i=0; i<n; i++ ) x[r[i]] = tmp[i];
550 
551   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
552   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
553   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
554   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
555 
556   PLogFlops(2*a->nz-a->n);
557   PetscFunctionReturn(0);
558 }
559 
560 #undef __FUNC__
561 #define __FUNC__ "MatSolveTransAdd_SeqAIJ"
562 int MatSolveTransAdd_SeqAIJ(Mat A,Vec bb, Vec zz,Vec xx)
563 {
564   Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data;
565   IS         iscol = a->col, isrow = a->row;
566   int        *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j;
567   int        nz,shift = a->indexshift, *rout, *cout, *diag = a->diag;
568   Scalar     *x,*b,*tmp, *aa = a->a, *v;
569 
570   PetscFunctionBegin;
571   if (zz != xx) VecCopy(zz,xx);
572 
573   ierr = VecGetArray(bb,&b);CHKERRQ(ierr);
574   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
575   tmp = a->solve_work;
576 
577   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
578   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
579 
580   /* copy the b into temp work space according to permutation */
581   for ( i=0; i<n; i++ ) tmp[i] = b[c[i]];
582 
583   /* forward solve the U^T */
584   for ( i=0; i<n; i++ ) {
585     v   = aa + diag[i] + shift;
586     vi  = aj + diag[i] + (!shift);
587     nz  = ai[i+1] - diag[i] - 1;
588     tmp[i] *= *v++;
589     while (nz--) {
590       tmp[*vi++ + shift] -= (*v++)*tmp[i];
591     }
592   }
593 
594   /* backward solve the L^T */
595   for ( i=n-1; i>=0; i-- ){
596     v   = aa + diag[i] - 1 + shift;
597     vi  = aj + diag[i] - 1 + shift;
598     nz  = diag[i] - ai[i];
599     while (nz--) {
600       tmp[*vi-- + shift] -= (*v--)*tmp[i];
601     }
602   }
603 
604   /* copy tmp into x according to permutation */
605   for ( i=0; i<n; i++ ) x[r[i]] += tmp[i];
606 
607   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
608   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
609   ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);
610   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
611 
612   PLogFlops(2*a->nz);
613   PetscFunctionReturn(0);
614 }
615 /* ----------------------------------------------------------------*/
616 extern int MatMissingDiag_SeqAIJ(Mat);
617 
618 #undef __FUNC__
619 #define __FUNC__ "MatILUFactorSymbolic_SeqAIJ"
620 int MatILUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatILUInfo *info,Mat *fact)
621 {
622   Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data, *b;
623   IS         isicol;
624   int        *r,*ic, ierr, prow, n = a->m, *ai = a->i, *aj = a->j;
625   int        *ainew,*ajnew, jmax,*fill, *xi, nz, *im,*ajfill,*flev;
626   int        *dloc, idx, row,m,fm, nzf, nzi,len,  realloc = 0, dcount = 0;
627   int        incrlev,nnz,i,shift = a->indexshift,levels,diagonal_fill;
628   PetscTruth col_identity, row_identity;
629   double     f;
630 
631   PetscFunctionBegin;
632   if (info) {
633     f             = info->fill;
634     levels        = (int) info->levels;
635     diagonal_fill = (int) info->diagonal_fill;
636   } else {
637     f             = 1.0;
638     levels        = 0;
639     diagonal_fill = 0;
640   }
641   ierr = ISInvertPermutation(iscol,&isicol);CHKERRQ(ierr);
642 
643   /* special case that simply copies fill pattern */
644   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
645   ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
646   if (levels == 0 && row_identity && col_identity) {
647     ierr = MatDuplicate_SeqAIJ(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr);
648     (*fact)->factor = FACTOR_LU;
649     b               = (Mat_SeqAIJ *) (*fact)->data;
650     if (!b->diag) {
651       ierr = MatMarkDiag_SeqAIJ(*fact);CHKERRQ(ierr);
652     }
653     ierr = MatMissingDiag_SeqAIJ(*fact);CHKERRQ(ierr);
654     b->row              = isrow;
655     b->col              = iscol;
656     b->icol             = isicol;
657     b->solve_work       = (Scalar *) PetscMalloc((b->m+1)*sizeof(Scalar));CHKPTRQ(b->solve_work);
658     (*fact)->ops->solve = MatSolve_SeqAIJ_NaturalOrdering;
659     ierr                = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
660     ierr                = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
661     PetscFunctionReturn(0);
662   }
663 
664   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
665   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
666 
667   /* get new row pointers */
668   ainew = (int *) PetscMalloc( (n+1)*sizeof(int) );CHKPTRQ(ainew);
669   ainew[0] = -shift;
670   /* don't know how many column pointers are needed so estimate */
671   jmax = (int) (f*(ai[n]+!shift));
672   ajnew = (int *) PetscMalloc( (jmax)*sizeof(int) );CHKPTRQ(ajnew);
673   /* ajfill is level of fill for each fill entry */
674   ajfill = (int *) PetscMalloc( (jmax)*sizeof(int) );CHKPTRQ(ajfill);
675   /* fill is a linked list of nonzeros in active row */
676   fill = (int *) PetscMalloc( (n+1)*sizeof(int));CHKPTRQ(fill);
677   /* im is level for each filled value */
678   im = (int *) PetscMalloc( (n+1)*sizeof(int));CHKPTRQ(im);
679   /* dloc is location of diagonal in factor */
680   dloc = (int *) PetscMalloc( (n+1)*sizeof(int));CHKPTRQ(dloc);
681   dloc[0]  = 0;
682   for ( prow=0; prow<n; prow++ ) {
683 
684     /* copy current row into linked list */
685     nzf     = nz  = ai[r[prow]+1] - ai[r[prow]];
686     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,1,"Empty row in matrix");
687     xi      = aj + ai[r[prow]] + shift;
688     fill[n]    = n;
689     fill[prow] = -1; /* marker to indicate if diagonal exists */
690     while (nz--) {
691       fm  = n;
692       idx = ic[*xi++ + shift];
693       do {
694         m  = fm;
695         fm = fill[m];
696       } while (fm < idx);
697       fill[m]   = idx;
698       fill[idx] = fm;
699       im[idx]   = 0;
700     }
701 
702     /* make sure diagonal entry is included */
703     if (diagonal_fill && fill[prow] == -1) {
704       fm = n;
705       while (fill[fm] < prow) fm = fill[fm];
706       fill[prow] = fill[fm]; /* insert diagonal into linked list */
707       fill[fm]   = prow;
708       im[prow]   = 0;
709       nzf++;
710       dcount++;
711     }
712 
713     nzi = 0;
714     row = fill[n];
715     while ( row < prow ) {
716       incrlev = im[row] + 1;
717       nz      = dloc[row];
718       xi      = ajnew  + ainew[row] + shift + nz + 1;
719       flev    = ajfill + ainew[row] + shift + nz + 1;
720       nnz     = ainew[row+1] - ainew[row] - nz - 1;
721       fm      = row;
722       while (nnz-- > 0) {
723         idx = *xi++ + shift;
724         if (*flev + incrlev > levels) {
725           flev++;
726           continue;
727         }
728         do {
729           m  = fm;
730           fm = fill[m];
731         } while (fm < idx);
732         if (fm != idx) {
733           im[idx]   = *flev + incrlev;
734           fill[m]   = idx;
735           fill[idx] = fm;
736           fm        = idx;
737           nzf++;
738         } else {
739           if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev;
740         }
741         flev++;
742       }
743       row = fill[row];
744       nzi++;
745     }
746     /* copy new filled row into permanent storage */
747     ainew[prow+1] = ainew[prow] + nzf;
748     if (ainew[prow+1] > jmax-shift) {
749 
750       /* estimate how much additional space we will need */
751       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
752       /* just double the memory each time */
753       /*  maxadd = (int) ((f*(ai[n]+!shift)*(n-prow+5))/n); */
754       int maxadd = jmax;
755       if (maxadd < nzf) maxadd = (n-prow)*(nzf+1);
756       jmax += maxadd;
757 
758       /* allocate a longer ajnew and ajfill */
759       xi     = (int *) PetscMalloc( jmax*sizeof(int) );CHKPTRQ(xi);
760       ierr   = PetscMemcpy(xi,ajnew,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr);
761       ierr = PetscFree(ajnew);CHKERRQ(ierr);
762       ajnew  = xi;
763       xi     = (int *) PetscMalloc( jmax*sizeof(int) );CHKPTRQ(xi);
764       ierr   = PetscMemcpy(xi,ajfill,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr);
765       ierr = PetscFree(ajfill);CHKERRQ(ierr);
766       ajfill = xi;
767       realloc++; /* count how many times we realloc */
768     }
769     xi          = ajnew + ainew[prow] + shift;
770     flev        = ajfill + ainew[prow] + shift;
771     dloc[prow]  = nzi;
772     fm          = fill[n];
773     while (nzf--) {
774       *xi++   = fm - shift;
775       *flev++ = im[fm];
776       fm      = fill[fm];
777     }
778     /* make sure row has diagonal entry */
779     if (ajnew[ainew[prow]+shift+dloc[prow]]+shift != prow) {
780       SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,1,"Row %d has missing diagonal in factored matrix\n\
781     try running with -pc_ilu_nonzeros_along_diagonal or -pc_ilu_diagonal_fill",prow);
782     }
783   }
784   ierr = PetscFree(ajfill); CHKERRQ(ierr);
785   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
786   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
787   ierr = PetscFree(fill);CHKERRQ(ierr);
788   ierr = PetscFree(im);CHKERRQ(ierr);
789 
790   {
791     double af = ((double)ainew[n])/((double)ai[n]);
792     PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af);
793     PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Run with -pc_ilu_fill %g or use \n",af);
794     PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:PCILUSetFill(pc,%g);\n",af);
795     PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:for best performance.\n");
796     if (diagonal_fill) {
797       PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Detected and replace %d missing diagonals",dcount);
798     }
799   }
800 
801   /* put together the new matrix */
802   ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,fact);CHKERRQ(ierr);
803   PLogObjectParent(*fact,isicol);
804   b = (Mat_SeqAIJ *) (*fact)->data;
805   ierr = PetscFree(b->imax);CHKERRQ(ierr);
806   b->singlemalloc = 0;
807   len = (ainew[n] + shift)*sizeof(Scalar);
808   /* the next line frees the default space generated by the Create() */
809   ierr = PetscFree(b->a);CHKERRQ(ierr);
810   ierr = PetscFree(b->ilen);CHKERRQ(ierr);
811   b->a          = (Scalar *) PetscMalloc( len+1 );CHKPTRQ(b->a);
812   b->j          = ajnew;
813   b->i          = ainew;
814   for ( i=0; i<n; i++ ) dloc[i] += ainew[i];
815   b->diag       = dloc;
816   b->ilen       = 0;
817   b->imax       = 0;
818   b->row        = isrow;
819   b->col        = iscol;
820   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
821   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
822   b->icol       = isicol;
823   b->solve_work = (Scalar *) PetscMalloc( (n+1)*sizeof(Scalar));CHKPTRQ(b->solve_work);
824   /* In b structure:  Free imax, ilen, old a, old j.
825      Allocate dloc, solve_work, new a, new j */
826   PLogObjectMemory(*fact,(ainew[n]+shift-n) * (sizeof(int)+sizeof(Scalar)));
827   b->maxnz          = b->nz = ainew[n] + shift;
828   (*fact)->factor   = FACTOR_LU;
829 
830   (*fact)->info.factor_mallocs    = realloc;
831   (*fact)->info.fill_ratio_given  = f;
832   (*fact)->info.fill_ratio_needed = ((double)ainew[n])/((double)ai[prow]);
833   (*fact)->factor                 =  FACTOR_LU;;
834 
835   PetscFunctionReturn(0);
836 }
837 
838 
839 
840 
841