1 
2 /*
3     Factorization code for BAIJ format.
4 */
5 #include <../src/mat/impls/baij/seq/baij.h>
6 #include <petsc/private/kernels/blockinvert.h>
7 
MatLUFactorNumeric_SeqBAIJ_2(Mat B,Mat A,const MatFactorInfo * info)8 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2(Mat B,Mat A,const MatFactorInfo *info)
9 {
10   Mat            C     =B;
11   Mat_SeqBAIJ    *a    =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data;
12   IS             isrow = b->row,isicol = b->icol;
13   PetscErrorCode ierr;
14   const PetscInt *r,*ic;
15   PetscInt       i,j,k,nz,nzL,row,*pj;
16   const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,bs2=a->bs2;
17   const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag;
18   MatScalar      *rtmp,*pc,*mwork,*pv;
19   MatScalar      *aa=a->a,*v;
20   PetscInt       flg;
21   PetscReal      shift = info->shiftamount;
22   PetscBool      allowzeropivot,zeropivotdetected;
23 
24   PetscFunctionBegin;
25   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
26   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
27   allowzeropivot = PetscNot(A->erroriffailure);
28 
29   /* generate work space needed by the factorization */
30   ierr = PetscMalloc2(bs2*n,&rtmp,bs2,&mwork);CHKERRQ(ierr);
31   ierr = PetscArrayzero(rtmp,bs2*n);CHKERRQ(ierr);
32 
33   for (i=0; i<n; i++) {
34     /* zero rtmp */
35     /* L part */
36     nz    = bi[i+1] - bi[i];
37     bjtmp = bj + bi[i];
38     for  (j=0; j<nz; j++) {
39       ierr = PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);CHKERRQ(ierr);
40     }
41 
42     /* U part */
43     nz    = bdiag[i] - bdiag[i+1];
44     bjtmp = bj + bdiag[i+1]+1;
45     for  (j=0; j<nz; j++) {
46       ierr = PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);CHKERRQ(ierr);
47     }
48 
49     /* load in initial (unfactored row) */
50     nz    = ai[r[i]+1] - ai[r[i]];
51     ajtmp = aj + ai[r[i]];
52     v     = aa + bs2*ai[r[i]];
53     for (j=0; j<nz; j++) {
54       ierr = PetscArraycpy(rtmp+bs2*ic[ajtmp[j]],v+bs2*j,bs2);CHKERRQ(ierr);
55     }
56 
57     /* elimination */
58     bjtmp = bj + bi[i];
59     nzL   = bi[i+1] - bi[i];
60     for (k=0; k < nzL; k++) {
61       row = bjtmp[k];
62       pc  = rtmp + bs2*row;
63       for (flg=0,j=0; j<bs2; j++) {
64         if (pc[j] != (PetscScalar)0.0) {
65           flg = 1;
66           break;
67         }
68       }
69       if (flg) {
70         pv = b->a + bs2*bdiag[row];
71         /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
72         ierr = PetscKernel_A_gets_A_times_B_2(pc,pv,mwork);CHKERRQ(ierr);
73 
74         pj = b->j + bdiag[row+1]+1; /* begining of U(row,:) */
75         pv = b->a + bs2*(bdiag[row+1]+1);
76         nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */
77         for (j=0; j<nz; j++) {
78           /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
79           /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
80           v    = rtmp + 4*pj[j];
81           ierr = PetscKernel_A_gets_A_minus_B_times_C_2(v,pc,pv);CHKERRQ(ierr);
82           pv  += 4;
83         }
84         ierr = PetscLogFlops(16.0*nz+12);CHKERRQ(ierr); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
85       }
86     }
87 
88     /* finished row so stick it into b->a */
89     /* L part */
90     pv = b->a + bs2*bi[i];
91     pj = b->j + bi[i];
92     nz = bi[i+1] - bi[i];
93     for (j=0; j<nz; j++) {
94       ierr = PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);CHKERRQ(ierr);
95     }
96 
97     /* Mark diagonal and invert diagonal for simplier triangular solves */
98     pv   = b->a + bs2*bdiag[i];
99     pj   = b->j + bdiag[i];
100     ierr = PetscArraycpy(pv,rtmp+bs2*pj[0],bs2);CHKERRQ(ierr);
101     ierr = PetscKernel_A_gets_inverse_A_2(pv,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr);
102     if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
103 
104     /* U part */
105     pv = b->a + bs2*(bdiag[i+1]+1);
106     pj = b->j + bdiag[i+1]+1;
107     nz = bdiag[i] - bdiag[i+1] - 1;
108     for (j=0; j<nz; j++) {
109       ierr = PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);CHKERRQ(ierr);
110     }
111   }
112 
113   ierr = PetscFree2(rtmp,mwork);CHKERRQ(ierr);
114   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
115   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
116 
117   C->ops->solve          = MatSolve_SeqBAIJ_2;
118   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2;
119   C->assembled           = PETSC_TRUE;
120 
121   ierr = PetscLogFlops(1.333333333333*2*2*2*n);CHKERRQ(ierr); /* from inverting diagonal blocks */
122   PetscFunctionReturn(0);
123 }
124 
MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering(Mat B,Mat A,const MatFactorInfo * info)125 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info)
126 {
127   Mat            C =B;
128   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data;
129   PetscErrorCode ierr;
130   PetscInt       i,j,k,nz,nzL,row,*pj;
131   const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,bs2=a->bs2;
132   const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag;
133   MatScalar      *rtmp,*pc,*mwork,*pv;
134   MatScalar      *aa=a->a,*v;
135   PetscInt       flg;
136   PetscReal      shift = info->shiftamount;
137   PetscBool      allowzeropivot,zeropivotdetected;
138 
139   PetscFunctionBegin;
140   allowzeropivot = PetscNot(A->erroriffailure);
141 
142   /* generate work space needed by the factorization */
143   ierr = PetscMalloc2(bs2*n,&rtmp,bs2,&mwork);CHKERRQ(ierr);
144   ierr = PetscArrayzero(rtmp,bs2*n);CHKERRQ(ierr);
145 
146   for (i=0; i<n; i++) {
147     /* zero rtmp */
148     /* L part */
149     nz    = bi[i+1] - bi[i];
150     bjtmp = bj + bi[i];
151     for  (j=0; j<nz; j++) {
152       ierr = PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);CHKERRQ(ierr);
153     }
154 
155     /* U part */
156     nz    = bdiag[i] - bdiag[i+1];
157     bjtmp = bj + bdiag[i+1]+1;
158     for  (j=0; j<nz; j++) {
159       ierr = PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);CHKERRQ(ierr);
160     }
161 
162     /* load in initial (unfactored row) */
163     nz    = ai[i+1] - ai[i];
164     ajtmp = aj + ai[i];
165     v     = aa + bs2*ai[i];
166     for (j=0; j<nz; j++) {
167       ierr = PetscArraycpy(rtmp+bs2*ajtmp[j],v+bs2*j,bs2);CHKERRQ(ierr);
168     }
169 
170     /* elimination */
171     bjtmp = bj + bi[i];
172     nzL   = bi[i+1] - bi[i];
173     for (k=0; k < nzL; k++) {
174       row = bjtmp[k];
175       pc  = rtmp + bs2*row;
176       for (flg=0,j=0; j<bs2; j++) {
177         if (pc[j]!=(PetscScalar)0.0) {
178           flg = 1;
179           break;
180         }
181       }
182       if (flg) {
183         pv = b->a + bs2*bdiag[row];
184         /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
185         ierr = PetscKernel_A_gets_A_times_B_2(pc,pv,mwork);CHKERRQ(ierr);
186 
187         pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
188         pv = b->a + bs2*(bdiag[row+1]+1);
189         nz = bdiag[row]-bdiag[row+1] - 1; /* num of entries in U(row,:) excluding diag */
190         for (j=0; j<nz; j++) {
191           /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
192           /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
193           v    = rtmp + 4*pj[j];
194           ierr = PetscKernel_A_gets_A_minus_B_times_C_2(v,pc,pv);CHKERRQ(ierr);
195           pv  += 4;
196         }
197         ierr = PetscLogFlops(16.0*nz+12);CHKERRQ(ierr); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
198       }
199     }
200 
201     /* finished row so stick it into b->a */
202     /* L part */
203     pv = b->a + bs2*bi[i];
204     pj = b->j + bi[i];
205     nz = bi[i+1] - bi[i];
206     for (j=0; j<nz; j++) {
207       ierr = PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);CHKERRQ(ierr);
208     }
209 
210     /* Mark diagonal and invert diagonal for simplier triangular solves */
211     pv   = b->a + bs2*bdiag[i];
212     pj   = b->j + bdiag[i];
213     ierr = PetscArraycpy(pv,rtmp+bs2*pj[0],bs2);CHKERRQ(ierr);
214     ierr = PetscKernel_A_gets_inverse_A_2(pv,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr);
215     if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
216 
217     /* U part */
218     /*
219     pv = b->a + bs2*bi[2*n-i];
220     pj = b->j + bi[2*n-i];
221     nz = bi[2*n-i+1] - bi[2*n-i] - 1;
222     */
223     pv = b->a + bs2*(bdiag[i+1]+1);
224     pj = b->j + bdiag[i+1]+1;
225     nz = bdiag[i] - bdiag[i+1] - 1;
226     for (j=0; j<nz; j++) {
227       ierr = PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);CHKERRQ(ierr);
228     }
229   }
230   ierr = PetscFree2(rtmp,mwork);CHKERRQ(ierr);
231 
232   C->ops->solve          = MatSolve_SeqBAIJ_2_NaturalOrdering;
233   C->ops->forwardsolve   = MatForwardSolve_SeqBAIJ_2_NaturalOrdering;
234   C->ops->backwardsolve  = MatBackwardSolve_SeqBAIJ_2_NaturalOrdering;
235   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering;
236   C->assembled           = PETSC_TRUE;
237 
238   ierr = PetscLogFlops(1.333333333333*2*2*2*n);CHKERRQ(ierr); /* from inverting diagonal blocks */
239   PetscFunctionReturn(0);
240 }
241 
MatLUFactorNumeric_SeqBAIJ_2_inplace(Mat B,Mat A,const MatFactorInfo * info)242 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_inplace(Mat B,Mat A,const MatFactorInfo *info)
243 {
244   Mat            C     = B;
245   Mat_SeqBAIJ    *a    = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data;
246   IS             isrow = b->row,isicol = b->icol;
247   PetscErrorCode ierr;
248   const PetscInt *r,*ic;
249   PetscInt       i,j,n = a->mbs,*bi = b->i,*bj = b->j;
250   PetscInt       *ajtmpold,*ajtmp,nz,row;
251   PetscInt       *diag_offset=b->diag,idx,*ai=a->i,*aj=a->j,*pj;
252   MatScalar      *pv,*v,*rtmp,m1,m2,m3,m4,*pc,*w,*x,x1,x2,x3,x4;
253   MatScalar      p1,p2,p3,p4;
254   MatScalar      *ba   = b->a,*aa = a->a;
255   PetscReal      shift = info->shiftamount;
256   PetscBool      allowzeropivot,zeropivotdetected;
257 
258   PetscFunctionBegin;
259   allowzeropivot = PetscNot(A->erroriffailure);
260   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
261   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
262   ierr = PetscMalloc1(4*(n+1),&rtmp);CHKERRQ(ierr);
263 
264   for (i=0; i<n; i++) {
265     nz    = bi[i+1] - bi[i];
266     ajtmp = bj + bi[i];
267     for  (j=0; j<nz; j++) {
268       x = rtmp+4*ajtmp[j]; x[0] = x[1] = x[2] = x[3] = 0.0;
269     }
270     /* load in initial (unfactored row) */
271     idx      = r[i];
272     nz       = ai[idx+1] - ai[idx];
273     ajtmpold = aj + ai[idx];
274     v        = aa + 4*ai[idx];
275     for (j=0; j<nz; j++) {
276       x    = rtmp+4*ic[ajtmpold[j]];
277       x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3];
278       v   += 4;
279     }
280     row = *ajtmp++;
281     while (row < i) {
282       pc = rtmp + 4*row;
283       p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3];
284       if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) {
285         pv    = ba + 4*diag_offset[row];
286         pj    = bj + diag_offset[row] + 1;
287         x1    = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3];
288         pc[0] = m1 = p1*x1 + p3*x2;
289         pc[1] = m2 = p2*x1 + p4*x2;
290         pc[2] = m3 = p1*x3 + p3*x4;
291         pc[3] = m4 = p2*x3 + p4*x4;
292         nz    = bi[row+1] - diag_offset[row] - 1;
293         pv   += 4;
294         for (j=0; j<nz; j++) {
295           x1    = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3];
296           x     = rtmp + 4*pj[j];
297           x[0] -= m1*x1 + m3*x2;
298           x[1] -= m2*x1 + m4*x2;
299           x[2] -= m1*x3 + m3*x4;
300           x[3] -= m2*x3 + m4*x4;
301           pv   += 4;
302         }
303         ierr = PetscLogFlops(16.0*nz+12.0);CHKERRQ(ierr);
304       }
305       row = *ajtmp++;
306     }
307     /* finished row so stick it into b->a */
308     pv = ba + 4*bi[i];
309     pj = bj + bi[i];
310     nz = bi[i+1] - bi[i];
311     for (j=0; j<nz; j++) {
312       x     = rtmp+4*pj[j];
313       pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3];
314       pv   += 4;
315     }
316     /* invert diagonal block */
317     w    = ba + 4*diag_offset[i];
318     ierr = PetscKernel_A_gets_inverse_A_2(w,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr);
319     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
320   }
321 
322   ierr = PetscFree(rtmp);CHKERRQ(ierr);
323   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
324   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
325 
326   C->ops->solve          = MatSolve_SeqBAIJ_2_inplace;
327   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_inplace;
328   C->assembled           = PETSC_TRUE;
329 
330   ierr = PetscLogFlops(1.333333333333*8*b->mbs);CHKERRQ(ierr); /* from inverting diagonal blocks */
331   PetscFunctionReturn(0);
332 }
333 /*
334       Version for when blocks are 2 by 2 Using natural ordering
335 */
MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo * info)336 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo *info)
337 {
338   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data;
339   PetscErrorCode ierr;
340   PetscInt       i,j,n = a->mbs,*bi = b->i,*bj = b->j;
341   PetscInt       *ajtmpold,*ajtmp,nz,row;
342   PetscInt       *diag_offset = b->diag,*ai=a->i,*aj=a->j,*pj;
343   MatScalar      *pv,*v,*rtmp,*pc,*w,*x;
344   MatScalar      p1,p2,p3,p4,m1,m2,m3,m4,x1,x2,x3,x4;
345   MatScalar      *ba   = b->a,*aa = a->a;
346   PetscReal      shift = info->shiftamount;
347   PetscBool      allowzeropivot,zeropivotdetected;
348 
349   PetscFunctionBegin;
350   allowzeropivot = PetscNot(A->erroriffailure);
351   ierr = PetscMalloc1(4*(n+1),&rtmp);CHKERRQ(ierr);
352   for (i=0; i<n; i++) {
353     nz    = bi[i+1] - bi[i];
354     ajtmp = bj + bi[i];
355     for  (j=0; j<nz; j++) {
356       x    = rtmp+4*ajtmp[j];
357       x[0] = x[1]  = x[2]  = x[3]  = 0.0;
358     }
359     /* load in initial (unfactored row) */
360     nz       = ai[i+1] - ai[i];
361     ajtmpold = aj + ai[i];
362     v        = aa + 4*ai[i];
363     for (j=0; j<nz; j++) {
364       x    = rtmp+4*ajtmpold[j];
365       x[0] = v[0];  x[1]  = v[1];  x[2]  = v[2];  x[3]  = v[3];
366       v   += 4;
367     }
368     row = *ajtmp++;
369     while (row < i) {
370       pc = rtmp + 4*row;
371       p1 = pc[0];  p2  = pc[1];  p3  = pc[2];  p4  = pc[3];
372       if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) {
373         pv    = ba + 4*diag_offset[row];
374         pj    = bj + diag_offset[row] + 1;
375         x1    = pv[0];  x2  = pv[1];  x3  = pv[2];  x4  = pv[3];
376         pc[0] = m1 = p1*x1 + p3*x2;
377         pc[1] = m2 = p2*x1 + p4*x2;
378         pc[2] = m3 = p1*x3 + p3*x4;
379         pc[3] = m4 = p2*x3 + p4*x4;
380         nz    = bi[row+1] - diag_offset[row] - 1;
381         pv   += 4;
382         for (j=0; j<nz; j++) {
383           x1    = pv[0];  x2  = pv[1];   x3 = pv[2];  x4  = pv[3];
384           x     = rtmp + 4*pj[j];
385           x[0] -= m1*x1 + m3*x2;
386           x[1] -= m2*x1 + m4*x2;
387           x[2] -= m1*x3 + m3*x4;
388           x[3] -= m2*x3 + m4*x4;
389           pv   += 4;
390         }
391         ierr = PetscLogFlops(16.0*nz+12.0);CHKERRQ(ierr);
392       }
393       row = *ajtmp++;
394     }
395     /* finished row so stick it into b->a */
396     pv = ba + 4*bi[i];
397     pj = bj + bi[i];
398     nz = bi[i+1] - bi[i];
399     for (j=0; j<nz; j++) {
400       x     = rtmp+4*pj[j];
401       pv[0] = x[0];  pv[1]  = x[1];  pv[2]  = x[2];  pv[3]  = x[3];
402       /*
403       printf(" col %d:",pj[j]);
404       PetscInt j1;
405       for (j1=0; j1<4; j1++) printf(" %g,",*(pv+j1));
406       printf("\n");
407       */
408       pv += 4;
409     }
410     /* invert diagonal block */
411     w = ba + 4*diag_offset[i];
412     ierr = PetscKernel_A_gets_inverse_A_2(w,shift, allowzeropivot,&zeropivotdetected);CHKERRQ(ierr);
413     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
414   }
415 
416   ierr = PetscFree(rtmp);CHKERRQ(ierr);
417 
418   C->ops->solve          = MatSolve_SeqBAIJ_2_NaturalOrdering_inplace;
419   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering_inplace;
420   C->assembled           = PETSC_TRUE;
421 
422   ierr = PetscLogFlops(1.333333333333*8*b->mbs);CHKERRQ(ierr); /* from inverting diagonal blocks */
423   PetscFunctionReturn(0);
424 }
425 
426 /* ----------------------------------------------------------- */
427 /*
428      Version for when blocks are 1 by 1.
429 */
MatLUFactorNumeric_SeqBAIJ_1(Mat B,Mat A,const MatFactorInfo * info)430 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1(Mat B,Mat A,const MatFactorInfo *info)
431 {
432   Mat             C     =B;
433   Mat_SeqBAIJ     *a    =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data;
434   IS              isrow = b->row,isicol = b->icol;
435   PetscErrorCode  ierr;
436   const PetscInt  *r,*ic,*ics;
437   const PetscInt  n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bdiag=b->diag;
438   PetscInt        i,j,k,nz,nzL,row,*pj;
439   const PetscInt  *ajtmp,*bjtmp;
440   MatScalar       *rtmp,*pc,multiplier,*pv;
441   const MatScalar *aa=a->a,*v;
442   PetscBool       row_identity,col_identity;
443   FactorShiftCtx  sctx;
444   const PetscInt  *ddiag;
445   PetscReal       rs;
446   MatScalar       d;
447 
448   PetscFunctionBegin;
449   /* MatPivotSetUp(): initialize shift context sctx */
450   ierr = PetscMemzero(&sctx,sizeof(FactorShiftCtx));CHKERRQ(ierr);
451 
452   if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
453     ddiag          = a->diag;
454     sctx.shift_top = info->zeropivot;
455     for (i=0; i<n; i++) {
456       /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
457       d  = (aa)[ddiag[i]];
458       rs = -PetscAbsScalar(d) - PetscRealPart(d);
459       v  = aa+ai[i];
460       nz = ai[i+1] - ai[i];
461       for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
462       if (rs>sctx.shift_top) sctx.shift_top = rs;
463     }
464     sctx.shift_top *= 1.1;
465     sctx.nshift_max = 5;
466     sctx.shift_lo   = 0.;
467     sctx.shift_hi   = 1.;
468   }
469 
470   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
471   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
472   ierr = PetscMalloc1(n+1,&rtmp);CHKERRQ(ierr);
473   ics  = ic;
474 
475   do {
476     sctx.newshift = PETSC_FALSE;
477     for (i=0; i<n; i++) {
478       /* zero rtmp */
479       /* L part */
480       nz    = bi[i+1] - bi[i];
481       bjtmp = bj + bi[i];
482       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;
483 
484       /* U part */
485       nz    = bdiag[i]-bdiag[i+1];
486       bjtmp = bj + bdiag[i+1]+1;
487       for  (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;
488 
489       /* load in initial (unfactored row) */
490       nz    = ai[r[i]+1] - ai[r[i]];
491       ajtmp = aj + ai[r[i]];
492       v     = aa + ai[r[i]];
493       for (j=0; j<nz; j++) rtmp[ics[ajtmp[j]]] = v[j];
494 
495       /* ZeropivotApply() */
496       rtmp[i] += sctx.shift_amount;  /* shift the diagonal of the matrix */
497 
498       /* elimination */
499       bjtmp = bj + bi[i];
500       row   = *bjtmp++;
501       nzL   = bi[i+1] - bi[i];
502       for (k=0; k < nzL; k++) {
503         pc = rtmp + row;
504         if (*pc != (PetscScalar)0.0) {
505           pv         = b->a + bdiag[row];
506           multiplier = *pc * (*pv);
507           *pc        = multiplier;
508 
509           pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
510           pv = b->a + bdiag[row+1]+1;
511           nz = bdiag[row]-bdiag[row+1]-1; /* num of entries in U(row,:) excluding diag */
512           for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
513           ierr = PetscLogFlops(2.0*nz);CHKERRQ(ierr);
514         }
515         row = *bjtmp++;
516       }
517 
518       /* finished row so stick it into b->a */
519       rs = 0.0;
520       /* L part */
521       pv = b->a + bi[i];
522       pj = b->j + bi[i];
523       nz = bi[i+1] - bi[i];
524       for (j=0; j<nz; j++) {
525         pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
526       }
527 
528       /* U part */
529       pv = b->a + bdiag[i+1]+1;
530       pj = b->j + bdiag[i+1]+1;
531       nz = bdiag[i] - bdiag[i+1]-1;
532       for (j=0; j<nz; j++) {
533         pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
534       }
535 
536       sctx.rs = rs;
537       sctx.pv = rtmp[i];
538       ierr    = MatPivotCheck(B,A,info,&sctx,i);CHKERRQ(ierr);
539       if (sctx.newshift) break; /* break for-loop */
540       rtmp[i] = sctx.pv; /* sctx.pv might be updated in the case of MAT_SHIFT_INBLOCKS */
541 
542       /* Mark diagonal and invert diagonal for simplier triangular solves */
543       pv  = b->a + bdiag[i];
544       *pv = (PetscScalar)1.0/rtmp[i];
545 
546     } /* endof for (i=0; i<n; i++) { */
547 
548     /* MatPivotRefine() */
549     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
550       /*
551        * if no shift in this attempt & shifting & started shifting & can refine,
552        * then try lower shift
553        */
554       sctx.shift_hi       = sctx.shift_fraction;
555       sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
556       sctx.shift_amount   = sctx.shift_fraction * sctx.shift_top;
557       sctx.newshift       = PETSC_TRUE;
558       sctx.nshift++;
559     }
560   } while (sctx.newshift);
561 
562   ierr = PetscFree(rtmp);CHKERRQ(ierr);
563   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
564   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
565 
566   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
567   ierr = ISIdentity(isicol,&col_identity);CHKERRQ(ierr);
568   if (row_identity && col_identity) {
569     C->ops->solve          = MatSolve_SeqBAIJ_1_NaturalOrdering;
570     C->ops->forwardsolve   = MatForwardSolve_SeqBAIJ_1_NaturalOrdering;
571     C->ops->backwardsolve  = MatBackwardSolve_SeqBAIJ_1_NaturalOrdering;
572     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering;
573   } else {
574     C->ops->solve          = MatSolve_SeqBAIJ_1;
575     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1;
576   }
577   C->assembled = PETSC_TRUE;
578   ierr         = PetscLogFlops(C->cmap->n);CHKERRQ(ierr);
579 
580   /* MatShiftView(A,info,&sctx) */
581   if (sctx.nshift) {
582     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
583       ierr = PetscInfo4(A,"number of shift_pd tries %D, shift_amount %g, diagonal shifted up by %e fraction top_value %e\n",sctx.nshift,(double)sctx.shift_amount,(double)sctx.shift_fraction,(double)sctx.shift_top);CHKERRQ(ierr);
584     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
585       ierr = PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);CHKERRQ(ierr);
586     } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
587       ierr = PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %g\n",sctx.nshift,(double)info->shiftamount);CHKERRQ(ierr);
588     }
589   }
590   PetscFunctionReturn(0);
591 }
592 
MatLUFactorNumeric_SeqBAIJ_1_inplace(Mat C,Mat A,const MatFactorInfo * info)593 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1_inplace(Mat C,Mat A,const MatFactorInfo *info)
594 {
595   Mat_SeqBAIJ    *a    = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data;
596   IS             isrow = b->row,isicol = b->icol;
597   PetscErrorCode ierr;
598   const PetscInt *r,*ic;
599   PetscInt       i,j,n = a->mbs,*bi = b->i,*bj = b->j;
600   PetscInt       *ajtmpold,*ajtmp,nz,row,*ai = a->i,*aj = a->j;
601   PetscInt       *diag_offset = b->diag,diag,*pj;
602   MatScalar      *pv,*v,*rtmp,multiplier,*pc;
603   MatScalar      *ba = b->a,*aa = a->a;
604   PetscBool      row_identity, col_identity;
605 
606   PetscFunctionBegin;
607   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
608   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
609   ierr = PetscMalloc1(n+1,&rtmp);CHKERRQ(ierr);
610 
611   for (i=0; i<n; i++) {
612     nz    = bi[i+1] - bi[i];
613     ajtmp = bj + bi[i];
614     for  (j=0; j<nz; j++) rtmp[ajtmp[j]] = 0.0;
615 
616     /* load in initial (unfactored row) */
617     nz       = ai[r[i]+1] - ai[r[i]];
618     ajtmpold = aj + ai[r[i]];
619     v        = aa + ai[r[i]];
620     for (j=0; j<nz; j++) rtmp[ic[ajtmpold[j]]] =  v[j];
621 
622     row = *ajtmp++;
623     while (row < i) {
624       pc = rtmp + row;
625       if (*pc != 0.0) {
626         pv         = ba + diag_offset[row];
627         pj         = bj + diag_offset[row] + 1;
628         multiplier = *pc * *pv++;
629         *pc        = multiplier;
630         nz         = bi[row+1] - diag_offset[row] - 1;
631         for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
632         ierr = PetscLogFlops(1.0+2.0*nz);CHKERRQ(ierr);
633       }
634       row = *ajtmp++;
635     }
636     /* finished row so stick it into b->a */
637     pv = ba + bi[i];
638     pj = bj + bi[i];
639     nz = bi[i+1] - bi[i];
640     for (j=0; j<nz; j++) pv[j] = rtmp[pj[j]];
641     diag = diag_offset[i] - bi[i];
642     /* check pivot entry for current row */
643     if (pv[diag] == 0.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot: row in original ordering %D in permuted ordering %D",r[i],i);
644     pv[diag] = 1.0/pv[diag];
645   }
646 
647   ierr = PetscFree(rtmp);CHKERRQ(ierr);
648   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
649   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
650   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
651   ierr = ISIdentity(isicol,&col_identity);CHKERRQ(ierr);
652   if (row_identity && col_identity) {
653     C->ops->solve          = MatSolve_SeqBAIJ_1_NaturalOrdering_inplace;
654     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering_inplace;
655   } else {
656     C->ops->solve          = MatSolve_SeqBAIJ_1_inplace;
657     C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_inplace;
658   }
659   C->assembled = PETSC_TRUE;
660   ierr         = PetscLogFlops(C->cmap->n);CHKERRQ(ierr);
661   PetscFunctionReturn(0);
662 }
663 
MatGetFactor_seqbaij_petsc(Mat A,MatFactorType ftype,Mat * B)664 PETSC_INTERN PetscErrorCode MatGetFactor_seqbaij_petsc(Mat A,MatFactorType ftype,Mat *B)
665 {
666   PetscInt       n = A->rmap->n;
667   PetscErrorCode ierr;
668 
669   PetscFunctionBegin;
670 #if defined(PETSC_USE_COMPLEX)
671   if (A->hermitian && (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported");
672 #endif
673   ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr);
674   ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr);
675   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU || ftype == MAT_FACTOR_ILUDT) {
676     ierr = MatSetType(*B,MATSEQBAIJ);CHKERRQ(ierr);
677 
678     (*B)->ops->lufactorsymbolic  = MatLUFactorSymbolic_SeqBAIJ;
679     (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqBAIJ;
680   } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
681     ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr);
682     ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr);
683 
684     (*B)->ops->iccfactorsymbolic      = MatICCFactorSymbolic_SeqBAIJ;
685     (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqBAIJ;
686   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
687   (*B)->factortype = ftype;
688   (*B)->useordering = PETSC_TRUE;
689 
690   ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr);
691   ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr);
692   PetscFunctionReturn(0);
693 }
694 
695 /* ----------------------------------------------------------- */
MatLUFactor_SeqBAIJ(Mat A,IS row,IS col,const MatFactorInfo * info)696 PetscErrorCode MatLUFactor_SeqBAIJ(Mat A,IS row,IS col,const MatFactorInfo *info)
697 {
698   PetscErrorCode ierr;
699   Mat            C;
700 
701   PetscFunctionBegin;
702   ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&C);CHKERRQ(ierr);
703   ierr = MatLUFactorSymbolic(C,A,row,col,info);CHKERRQ(ierr);
704   ierr = MatLUFactorNumeric(C,A,info);CHKERRQ(ierr);
705 
706   A->ops->solve          = C->ops->solve;
707   A->ops->solvetranspose = C->ops->solvetranspose;
708 
709   ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr);
710   ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)((Mat_SeqBAIJ*)(A->data))->icol);CHKERRQ(ierr);
711   PetscFunctionReturn(0);
712 }
713 
714 #include <../src/mat/impls/sbaij/seq/sbaij.h>
MatCholeskyFactorNumeric_SeqBAIJ_N(Mat C,Mat A,const MatFactorInfo * info)715 PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N(Mat C,Mat A,const MatFactorInfo *info)
716 {
717   PetscErrorCode ierr;
718   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data;
719   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
720   IS             ip=b->row;
721   const PetscInt *rip;
722   PetscInt       i,j,mbs=a->mbs,bs=A->rmap->bs,*bi=b->i,*bj=b->j,*bcol;
723   PetscInt       *ai=a->i,*aj=a->j;
724   PetscInt       k,jmin,jmax,*jl,*il,col,nexti,ili,nz;
725   MatScalar      *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
726   PetscReal      rs;
727   FactorShiftCtx sctx;
728 
729   PetscFunctionBegin;
730   if (bs > 1) { /* convert A to a SBAIJ matrix and apply Cholesky factorization from it */
731     if (!a->sbaijMat) {
732       ierr = MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);CHKERRQ(ierr);
733     }
734     ierr = (a->sbaijMat)->ops->choleskyfactornumeric(C,a->sbaijMat,info);CHKERRQ(ierr);
735     ierr = MatDestroy(&a->sbaijMat);CHKERRQ(ierr);
736     PetscFunctionReturn(0);
737   }
738 
739   /* MatPivotSetUp(): initialize shift context sctx */
740   ierr = PetscMemzero(&sctx,sizeof(FactorShiftCtx));CHKERRQ(ierr);
741 
742   ierr = ISGetIndices(ip,&rip);CHKERRQ(ierr);
743   ierr = PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&jl);CHKERRQ(ierr);
744 
745   sctx.shift_amount = 0.;
746   sctx.nshift       = 0;
747   do {
748     sctx.newshift = PETSC_FALSE;
749     for (i=0; i<mbs; i++) {
750       rtmp[i] = 0.0; jl[i] = mbs; il[0] = 0;
751     }
752 
753     for (k = 0; k<mbs; k++) {
754       bval = ba + bi[k];
755       /* initialize k-th row by the perm[k]-th row of A */
756       jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
757       for (j = jmin; j < jmax; j++) {
758         col = rip[aj[j]];
759         if (col >= k) { /* only take upper triangular entry */
760           rtmp[col] = aa[j];
761           *bval++   = 0.0; /* for in-place factorization */
762         }
763       }
764 
765       /* shift the diagonal of the matrix */
766       if (sctx.nshift) rtmp[k] += sctx.shift_amount;
767 
768       /* modify k-th row by adding in those rows i with U(i,k)!=0 */
769       dk = rtmp[k];
770       i  = jl[k]; /* first row to be added to k_th row  */
771 
772       while (i < k) {
773         nexti = jl[i]; /* next row to be added to k_th row */
774 
775         /* compute multiplier, update diag(k) and U(i,k) */
776         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
777         uikdi   = -ba[ili]*ba[bi[i]]; /* diagonal(k) */
778         dk     += uikdi*ba[ili];
779         ba[ili] = uikdi; /* -U(i,k) */
780 
781         /* add multiple of row i to k-th row */
782         jmin = ili + 1; jmax = bi[i+1];
783         if (jmin < jmax) {
784           for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
785           /* update il and jl for row i */
786           il[i] = jmin;
787           j     = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
788         }
789         i = nexti;
790       }
791 
792       /* shift the diagonals when zero pivot is detected */
793       /* compute rs=sum of abs(off-diagonal) */
794       rs   = 0.0;
795       jmin = bi[k]+1;
796       nz   = bi[k+1] - jmin;
797       if (nz) {
798         bcol = bj + jmin;
799         while (nz--) {
800           rs += PetscAbsScalar(rtmp[*bcol]);
801           bcol++;
802         }
803       }
804 
805       sctx.rs = rs;
806       sctx.pv = dk;
807       ierr    = MatPivotCheck(C,A,info,&sctx,k);CHKERRQ(ierr);
808       if (sctx.newshift) break;
809       dk = sctx.pv;
810 
811       /* copy data into U(k,:) */
812       ba[bi[k]] = 1.0/dk; /* U(k,k) */
813       jmin      = bi[k]+1; jmax = bi[k+1];
814       if (jmin < jmax) {
815         for (j=jmin; j<jmax; j++) {
816           col = bj[j]; ba[j] = rtmp[col]; rtmp[col] = 0.0;
817         }
818         /* add the k-th row into il and jl */
819         il[k] = jmin;
820         i     = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
821       }
822     }
823   } while (sctx.newshift);
824   ierr = PetscFree3(rtmp,il,jl);CHKERRQ(ierr);
825 
826   ierr = ISRestoreIndices(ip,&rip);CHKERRQ(ierr);
827 
828   C->assembled    = PETSC_TRUE;
829   C->preallocated = PETSC_TRUE;
830 
831   ierr = PetscLogFlops(C->rmap->N);CHKERRQ(ierr);
832   if (sctx.nshift) {
833     if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
834       ierr = PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);CHKERRQ(ierr);
835     } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
836       ierr = PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);CHKERRQ(ierr);
837     }
838   }
839   PetscFunctionReturn(0);
840 }
841 
MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering(Mat C,Mat A,const MatFactorInfo * info)842 PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info)
843 {
844   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data;
845   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
846   PetscErrorCode ierr;
847   PetscInt       i,j,am=a->mbs;
848   PetscInt       *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
849   PetscInt       k,jmin,*jl,*il,nexti,ili,*acol,*bcol,nz;
850   MatScalar      *rtmp,*ba=b->a,*aa=a->a,dk,uikdi,*aval,*bval;
851   PetscReal      rs;
852   FactorShiftCtx sctx;
853 
854   PetscFunctionBegin;
855   /* MatPivotSetUp(): initialize shift context sctx */
856   ierr = PetscMemzero(&sctx,sizeof(FactorShiftCtx));CHKERRQ(ierr);
857 
858   ierr = PetscMalloc3(am,&rtmp,am,&il,am,&jl);CHKERRQ(ierr);
859 
860   do {
861     sctx.newshift = PETSC_FALSE;
862     for (i=0; i<am; i++) {
863       rtmp[i] = 0.0; jl[i] = am; il[0] = 0;
864     }
865 
866     for (k = 0; k<am; k++) {
867       /* initialize k-th row with elements nonzero in row perm(k) of A */
868       nz   = ai[k+1] - ai[k];
869       acol = aj + ai[k];
870       aval = aa + ai[k];
871       bval = ba + bi[k];
872       while (nz--) {
873         if (*acol < k) { /* skip lower triangular entries */
874           acol++; aval++;
875         } else {
876           rtmp[*acol++] = *aval++;
877           *bval++       = 0.0; /* for in-place factorization */
878         }
879       }
880 
881       /* shift the diagonal of the matrix */
882       if (sctx.nshift) rtmp[k] += sctx.shift_amount;
883 
884       /* modify k-th row by adding in those rows i with U(i,k)!=0 */
885       dk = rtmp[k];
886       i  = jl[k]; /* first row to be added to k_th row  */
887 
888       while (i < k) {
889         nexti = jl[i]; /* next row to be added to k_th row */
890         /* compute multiplier, update D(k) and U(i,k) */
891         ili     = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
892         uikdi   = -ba[ili]*ba[bi[i]];
893         dk     += uikdi*ba[ili];
894         ba[ili] = uikdi; /* -U(i,k) */
895 
896         /* add multiple of row i to k-th row ... */
897         jmin = ili + 1;
898         nz   = bi[i+1] - jmin;
899         if (nz > 0) {
900           bcol = bj + jmin;
901           bval = ba + jmin;
902           while (nz--) rtmp[*bcol++] += uikdi*(*bval++);
903           /* update il and jl for i-th row */
904           il[i] = jmin;
905           j     = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
906         }
907         i = nexti;
908       }
909 
910       /* shift the diagonals when zero pivot is detected */
911       /* compute rs=sum of abs(off-diagonal) */
912       rs   = 0.0;
913       jmin = bi[k]+1;
914       nz   = bi[k+1] - jmin;
915       if (nz) {
916         bcol = bj + jmin;
917         while (nz--) {
918           rs += PetscAbsScalar(rtmp[*bcol]);
919           bcol++;
920         }
921       }
922 
923       sctx.rs = rs;
924       sctx.pv = dk;
925       ierr    = MatPivotCheck(C,A,info,&sctx,k);CHKERRQ(ierr);
926       if (sctx.newshift) break;    /* sctx.shift_amount is updated */
927       dk = sctx.pv;
928 
929       /* copy data into U(k,:) */
930       ba[bi[k]] = 1.0/dk;
931       jmin      = bi[k]+1;
932       nz        = bi[k+1] - jmin;
933       if (nz) {
934         bcol = bj + jmin;
935         bval = ba + jmin;
936         while (nz--) {
937           *bval++       = rtmp[*bcol];
938           rtmp[*bcol++] = 0.0;
939         }
940         /* add k-th row into il and jl */
941         il[k] = jmin;
942         i     = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
943       }
944     }
945   } while (sctx.newshift);
946   ierr = PetscFree3(rtmp,il,jl);CHKERRQ(ierr);
947 
948   C->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
949   C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
950   C->assembled           = PETSC_TRUE;
951   C->preallocated        = PETSC_TRUE;
952 
953   ierr = PetscLogFlops(C->rmap->N);CHKERRQ(ierr);
954   if (sctx.nshift) {
955     if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
956       ierr = PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);CHKERRQ(ierr);
957     } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
958       ierr = PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);CHKERRQ(ierr);
959     }
960   }
961   PetscFunctionReturn(0);
962 }
963 
964 #include <petscbt.h>
965 #include <../src/mat/utils/freespace.h>
MatICCFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo * info)966 PetscErrorCode MatICCFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
967 {
968   Mat_SeqBAIJ        *a = (Mat_SeqBAIJ*)A->data;
969   Mat_SeqSBAIJ       *b;
970   Mat                B;
971   PetscErrorCode     ierr;
972   PetscBool          perm_identity,missing;
973   PetscInt           reallocs=0,i,*ai=a->i,*aj=a->j,am=a->mbs,bs=A->rmap->bs,*ui;
974   const PetscInt     *rip;
975   PetscInt           jmin,jmax,nzk,k,j,*jl,prow,*il,nextprow;
976   PetscInt           nlnk,*lnk,*lnk_lvl=NULL,ncols,ncols_upper,*cols,*cols_lvl,*uj,**uj_ptr,**uj_lvl_ptr;
977   PetscReal          fill          =info->fill,levels=info->levels;
978   PetscFreeSpaceList free_space    =NULL,current_space=NULL;
979   PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL;
980   PetscBT            lnkbt;
981 
982   PetscFunctionBegin;
983   ierr = MatMissingDiagonal(A,&missing,&i);CHKERRQ(ierr);
984   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);
985 
986   if (bs > 1) {
987     if (!a->sbaijMat) {
988       ierr = MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);CHKERRQ(ierr);
989     }
990     (fact)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ;  /* undue the change made in MatGetFactor_seqbaij_petsc */
991 
992     ierr = MatICCFactorSymbolic(fact,a->sbaijMat,perm,info);CHKERRQ(ierr);
993     PetscFunctionReturn(0);
994   }
995 
996   ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr);
997   ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr);
998 
999   /* special case that simply copies fill pattern */
1000   if (!levels && perm_identity) {
1001     ierr = PetscMalloc1(am+1,&ui);CHKERRQ(ierr);
1002     for (i=0; i<am; i++) ui[i] = ai[i+1] - a->diag[i]; /* ui: rowlengths - changes when !perm_identity */
1003     B    = fact;
1004     ierr = MatSeqSBAIJSetPreallocation(B,1,0,ui);CHKERRQ(ierr);
1005 
1006 
1007     b  = (Mat_SeqSBAIJ*)B->data;
1008     uj = b->j;
1009     for (i=0; i<am; i++) {
1010       aj = a->j + a->diag[i];
1011       for (j=0; j<ui[i]; j++) *uj++ = *aj++;
1012       b->ilen[i] = ui[i];
1013     }
1014     ierr = PetscFree(ui);CHKERRQ(ierr);
1015 
1016     B->factortype = MAT_FACTOR_NONE;
1017 
1018     ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1019     ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1020     B->factortype = MAT_FACTOR_ICC;
1021 
1022     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1023     PetscFunctionReturn(0);
1024   }
1025 
1026   /* initialization */
1027   ierr  = PetscMalloc1(am+1,&ui);CHKERRQ(ierr);
1028   ui[0] = 0;
1029   ierr  = PetscMalloc1(2*am+1,&cols_lvl);CHKERRQ(ierr);
1030 
1031   /* jl: linked list for storing indices of the pivot rows
1032      il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
1033   ierr = PetscMalloc4(am,&uj_ptr,am,&uj_lvl_ptr,am,&il,am,&jl);CHKERRQ(ierr);
1034   for (i=0; i<am; i++) {
1035     jl[i] = am; il[i] = 0;
1036   }
1037 
1038   /* create and initialize a linked list for storing column indices of the active row k */
1039   nlnk = am + 1;
1040   ierr = PetscIncompleteLLCreate(am,am,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr);
1041 
1042   /* initial FreeSpace size is fill*(ai[am]+am)/2 */
1043   ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am]/2,am/2)),&free_space);CHKERRQ(ierr);
1044 
1045   current_space = free_space;
1046 
1047   ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am]/2,am/2)),&free_space_lvl);CHKERRQ(ierr);
1048   current_space_lvl = free_space_lvl;
1049 
1050   for (k=0; k<am; k++) {  /* for each active row k */
1051     /* initialize lnk by the column indices of row rip[k] of A */
1052     nzk         = 0;
1053     ncols       = ai[rip[k]+1] - ai[rip[k]];
1054     ncols_upper = 0;
1055     cols        = cols_lvl + am;
1056     for (j=0; j<ncols; j++) {
1057       i = rip[*(aj + ai[rip[k]] + j)];
1058       if (i >= k) { /* only take upper triangular entry */
1059         cols[ncols_upper]     = i;
1060         cols_lvl[ncols_upper] = -1;  /* initialize level for nonzero entries */
1061         ncols_upper++;
1062       }
1063     }
1064     ierr = PetscIncompleteLLAdd(ncols_upper,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr);
1065     nzk += nlnk;
1066 
1067     /* update lnk by computing fill-in for each pivot row to be merged in */
1068     prow = jl[k]; /* 1st pivot row */
1069 
1070     while (prow < k) {
1071       nextprow = jl[prow];
1072 
1073       /* merge prow into k-th row */
1074       jmin  = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
1075       jmax  = ui[prow+1];
1076       ncols = jmax-jmin;
1077       i     = jmin - ui[prow];
1078       cols  = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
1079       for (j=0; j<ncols; j++) cols_lvl[j] = *(uj_lvl_ptr[prow] + i + j);
1080       ierr = PetscIncompleteLLAddSorted(ncols,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr);
1081       nzk += nlnk;
1082 
1083       /* update il and jl for prow */
1084       if (jmin < jmax) {
1085         il[prow] = jmin;
1086 
1087         j = *cols; jl[prow] = jl[j]; jl[j] = prow;
1088       }
1089       prow = nextprow;
1090     }
1091 
1092     /* if free space is not available, make more free space */
1093     if (current_space->local_remaining<nzk) {
1094       i    = am - k + 1; /* num of unfactored rows */
1095       i    = PetscMin(PetscIntMultTruncate(i,nzk), PetscIntMultTruncate(i,i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1096       ierr = PetscFreeSpaceGet(i,&current_space);CHKERRQ(ierr);
1097       ierr = PetscFreeSpaceGet(i,&current_space_lvl);CHKERRQ(ierr);
1098       reallocs++;
1099     }
1100 
1101     /* copy data into free_space and free_space_lvl, then initialize lnk */
1102     ierr = PetscIncompleteLLClean(am,am,nzk,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);CHKERRQ(ierr);
1103 
1104     /* add the k-th row into il and jl */
1105     if (nzk-1 > 0) {
1106       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
1107       jl[k] = jl[i]; jl[i] = k;
1108       il[k] = ui[k] + 1;
1109     }
1110     uj_ptr[k]     = current_space->array;
1111     uj_lvl_ptr[k] = current_space_lvl->array;
1112 
1113     current_space->array           += nzk;
1114     current_space->local_used      += nzk;
1115     current_space->local_remaining -= nzk;
1116 
1117     current_space_lvl->array           += nzk;
1118     current_space_lvl->local_used      += nzk;
1119     current_space_lvl->local_remaining -= nzk;
1120 
1121     ui[k+1] = ui[k] + nzk;
1122   }
1123 
1124   ierr = ISRestoreIndices(perm,&rip);CHKERRQ(ierr);
1125   ierr = PetscFree4(uj_ptr,uj_lvl_ptr,il,jl);CHKERRQ(ierr);
1126   ierr = PetscFree(cols_lvl);CHKERRQ(ierr);
1127 
1128   /* copy free_space into uj and free free_space; set uj in new datastructure; */
1129   ierr = PetscMalloc1(ui[am]+1,&uj);CHKERRQ(ierr);
1130   ierr = PetscFreeSpaceContiguous(&free_space,uj);CHKERRQ(ierr);
1131   ierr = PetscIncompleteLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
1132   ierr = PetscFreeSpaceDestroy(free_space_lvl);CHKERRQ(ierr);
1133 
1134   /* put together the new matrix in MATSEQSBAIJ format */
1135   B    = fact;
1136   ierr = MatSeqSBAIJSetPreallocation(B,1,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr);
1137 
1138   b                = (Mat_SeqSBAIJ*)B->data;
1139   b->singlemalloc  = PETSC_FALSE;
1140   b->free_a        = PETSC_TRUE;
1141   b->free_ij       = PETSC_TRUE;
1142 
1143   ierr = PetscMalloc1(ui[am]+1,&b->a);CHKERRQ(ierr);
1144 
1145   b->j             = uj;
1146   b->i             = ui;
1147   b->diag          = NULL;
1148   b->ilen          = NULL;
1149   b->imax          = NULL;
1150   b->row           = perm;
1151   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
1152 
1153   ierr    = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr);
1154 
1155   b->icol = perm;
1156 
1157   ierr    = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr);
1158   ierr    = PetscMalloc1(am+1,&b->solve_work);CHKERRQ(ierr);
1159   ierr    = PetscLogObjectMemory((PetscObject)B,(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr);
1160 
1161   b->maxnz = b->nz = ui[am];
1162 
1163   B->info.factor_mallocs   = reallocs;
1164   B->info.fill_ratio_given = fill;
1165   if (ai[am] != 0.) {
1166     /* nonzeros in lower triangular part of A (includign diagonals)= (ai[am]+am)/2 */
1167     B->info.fill_ratio_needed = ((PetscReal)2*ui[am])/(ai[am]+am);
1168   } else {
1169     B->info.fill_ratio_needed = 0.0;
1170   }
1171 #if defined(PETSC_USE_INFO)
1172   if (ai[am] != 0) {
1173     PetscReal af = B->info.fill_ratio_needed;
1174     ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);CHKERRQ(ierr);
1175     ierr = PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);CHKERRQ(ierr);
1176     ierr = PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);CHKERRQ(ierr);
1177   } else {
1178     ierr = PetscInfo(A,"Empty matrix\n");CHKERRQ(ierr);
1179   }
1180 #endif
1181   if (perm_identity) {
1182     B->ops->solve                 = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1183     B->ops->solvetranspose        = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1184     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1185   } else {
1186     (fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
1187   }
1188   PetscFunctionReturn(0);
1189 }
1190 
MatCholeskyFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo * info)1191 PetscErrorCode MatCholeskyFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
1192 {
1193   Mat_SeqBAIJ        *a = (Mat_SeqBAIJ*)A->data;
1194   Mat_SeqSBAIJ       *b;
1195   Mat                B;
1196   PetscErrorCode     ierr;
1197   PetscBool          perm_identity,missing;
1198   PetscReal          fill = info->fill;
1199   const PetscInt     *rip;
1200   PetscInt           i,mbs=a->mbs,bs=A->rmap->bs,*ai=a->i,*aj=a->j,reallocs=0,prow;
1201   PetscInt           *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
1202   PetscInt           nlnk,*lnk,ncols,ncols_upper,*cols,*uj,**ui_ptr,*uj_ptr;
1203   PetscFreeSpaceList free_space=NULL,current_space=NULL;
1204   PetscBT            lnkbt;
1205 
1206   PetscFunctionBegin;
1207   if (bs > 1) { /* convert to seqsbaij */
1208     if (!a->sbaijMat) {
1209       ierr = MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);CHKERRQ(ierr);
1210     }
1211     (fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */
1212 
1213     ierr = MatCholeskyFactorSymbolic(fact,a->sbaijMat,perm,info);CHKERRQ(ierr);
1214     PetscFunctionReturn(0);
1215   }
1216 
1217   ierr = MatMissingDiagonal(A,&missing,&i);CHKERRQ(ierr);
1218   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);
1219 
1220   /* check whether perm is the identity mapping */
1221   ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr);
1222   if (!perm_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported");
1223   ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr);
1224 
1225   /* initialization */
1226   ierr  = PetscMalloc1(mbs+1,&ui);CHKERRQ(ierr);
1227   ui[0] = 0;
1228 
1229   /* jl: linked list for storing indices of the pivot rows
1230      il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
1231   ierr = PetscMalloc4(mbs,&ui_ptr,mbs,&il,mbs,&jl,mbs,&cols);CHKERRQ(ierr);
1232   for (i=0; i<mbs; i++) {
1233     jl[i] = mbs; il[i] = 0;
1234   }
1235 
1236   /* create and initialize a linked list for storing column indices of the active row k */
1237   nlnk = mbs + 1;
1238   ierr = PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr);
1239 
1240   /* initial FreeSpace size is fill* (ai[mbs]+mbs)/2 */
1241   ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[mbs]/2,mbs/2)),&free_space);CHKERRQ(ierr);
1242 
1243   current_space = free_space;
1244 
1245   for (k=0; k<mbs; k++) {  /* for each active row k */
1246     /* initialize lnk by the column indices of row rip[k] of A */
1247     nzk         = 0;
1248     ncols       = ai[rip[k]+1] - ai[rip[k]];
1249     ncols_upper = 0;
1250     for (j=0; j<ncols; j++) {
1251       i = rip[*(aj + ai[rip[k]] + j)];
1252       if (i >= k) { /* only take upper triangular entry */
1253         cols[ncols_upper] = i;
1254         ncols_upper++;
1255       }
1256     }
1257     ierr = PetscLLAdd(ncols_upper,cols,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr);
1258     nzk += nlnk;
1259 
1260     /* update lnk by computing fill-in for each pivot row to be merged in */
1261     prow = jl[k]; /* 1st pivot row */
1262 
1263     while (prow < k) {
1264       nextprow = jl[prow];
1265       /* merge prow into k-th row */
1266       jmin   = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
1267       jmax   = ui[prow+1];
1268       ncols  = jmax-jmin;
1269       uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
1270       ierr   = PetscLLAddSorted(ncols,uj_ptr,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr);
1271       nzk   += nlnk;
1272 
1273       /* update il and jl for prow */
1274       if (jmin < jmax) {
1275         il[prow] = jmin;
1276         j        = *uj_ptr;
1277         jl[prow] = jl[j];
1278         jl[j]    = prow;
1279       }
1280       prow = nextprow;
1281     }
1282 
1283     /* if free space is not available, make more free space */
1284     if (current_space->local_remaining<nzk) {
1285       i    = mbs - k + 1; /* num of unfactored rows */
1286       i    = PetscMin(PetscIntMultTruncate(i,nzk), PetscIntMultTruncate(i,i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1287       ierr = PetscFreeSpaceGet(i,&current_space);CHKERRQ(ierr);
1288       reallocs++;
1289     }
1290 
1291     /* copy data into free space, then initialize lnk */
1292     ierr = PetscLLClean(mbs,mbs,nzk,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
1293 
1294     /* add the k-th row into il and jl */
1295     if (nzk-1 > 0) {
1296       i     = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
1297       jl[k] = jl[i]; jl[i] = k;
1298       il[k] = ui[k] + 1;
1299     }
1300     ui_ptr[k]                       = current_space->array;
1301     current_space->array           += nzk;
1302     current_space->local_used      += nzk;
1303     current_space->local_remaining -= nzk;
1304 
1305     ui[k+1] = ui[k] + nzk;
1306   }
1307 
1308   ierr = ISRestoreIndices(perm,&rip);CHKERRQ(ierr);
1309   ierr = PetscFree4(ui_ptr,il,jl,cols);CHKERRQ(ierr);
1310 
1311   /* copy free_space into uj and free free_space; set uj in new datastructure; */
1312   ierr = PetscMalloc1(ui[mbs]+1,&uj);CHKERRQ(ierr);
1313   ierr = PetscFreeSpaceContiguous(&free_space,uj);CHKERRQ(ierr);
1314   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
1315 
1316   /* put together the new matrix in MATSEQSBAIJ format */
1317   B    = fact;
1318   ierr = MatSeqSBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr);
1319 
1320   b               = (Mat_SeqSBAIJ*)B->data;
1321   b->singlemalloc = PETSC_FALSE;
1322   b->free_a       = PETSC_TRUE;
1323   b->free_ij      = PETSC_TRUE;
1324 
1325   ierr = PetscMalloc1(ui[mbs]+1,&b->a);CHKERRQ(ierr);
1326 
1327   b->j             = uj;
1328   b->i             = ui;
1329   b->diag          = NULL;
1330   b->ilen          = NULL;
1331   b->imax          = NULL;
1332   b->row           = perm;
1333   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
1334 
1335   ierr     = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr);
1336   b->icol  = perm;
1337   ierr     = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr);
1338   ierr     = PetscMalloc1(mbs+1,&b->solve_work);CHKERRQ(ierr);
1339   ierr     = PetscLogObjectMemory((PetscObject)B,(ui[mbs]-mbs)*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr);
1340   b->maxnz = b->nz = ui[mbs];
1341 
1342   B->info.factor_mallocs   = reallocs;
1343   B->info.fill_ratio_given = fill;
1344   if (ai[mbs] != 0.) {
1345     /* nonzeros in lower triangular part of A = (ai[mbs]+mbs)/2 */
1346     B->info.fill_ratio_needed = ((PetscReal)2*ui[mbs])/(ai[mbs]+mbs);
1347   } else {
1348     B->info.fill_ratio_needed = 0.0;
1349   }
1350 #if defined(PETSC_USE_INFO)
1351   if (ai[mbs] != 0.) {
1352     PetscReal af = B->info.fill_ratio_needed;
1353     ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);CHKERRQ(ierr);
1354     ierr = PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);CHKERRQ(ierr);
1355     ierr = PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);CHKERRQ(ierr);
1356   } else {
1357     ierr = PetscInfo(A,"Empty matrix\n");CHKERRQ(ierr);
1358   }
1359 #endif
1360   if (perm_identity) {
1361     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1362   } else {
1363     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
1364   }
1365   PetscFunctionReturn(0);
1366 }
1367 
MatSolve_SeqBAIJ_N_NaturalOrdering(Mat A,Vec bb,Vec xx)1368 PetscErrorCode MatSolve_SeqBAIJ_N_NaturalOrdering(Mat A,Vec bb,Vec xx)
1369 {
1370   Mat_SeqBAIJ       *a=(Mat_SeqBAIJ*)A->data;
1371   PetscErrorCode    ierr;
1372   const PetscInt    *ai=a->i,*aj=a->j,*adiag=a->diag,*vi;
1373   PetscInt          i,k,n=a->mbs;
1374   PetscInt          nz,bs=A->rmap->bs,bs2=a->bs2;
1375   const MatScalar   *aa=a->a,*v;
1376   PetscScalar       *x,*s,*t,*ls;
1377   const PetscScalar *b;
1378 
1379   PetscFunctionBegin;
1380   ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr);
1381   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
1382   t    = a->solve_work;
1383 
1384   /* forward solve the lower triangular */
1385   ierr = PetscArraycpy(t,b,bs);CHKERRQ(ierr); /* copy 1st block of b to t */
1386 
1387   for (i=1; i<n; i++) {
1388     v    = aa + bs2*ai[i];
1389     vi   = aj + ai[i];
1390     nz   = ai[i+1] - ai[i];
1391     s    = t + bs*i;
1392     ierr = PetscArraycpy(s,b+bs*i,bs);CHKERRQ(ierr); /* copy i_th block of b to t */
1393     for (k=0;k<nz;k++) {
1394       PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[k]);
1395       v += bs2;
1396     }
1397   }
1398 
1399   /* backward solve the upper triangular */
1400   ls = a->solve_work + A->cmap->n;
1401   for (i=n-1; i>=0; i--) {
1402     v    = aa + bs2*(adiag[i+1]+1);
1403     vi   = aj + adiag[i+1]+1;
1404     nz   = adiag[i] - adiag[i+1]-1;
1405     ierr = PetscArraycpy(ls,t+i*bs,bs);CHKERRQ(ierr);
1406     for (k=0; k<nz; k++) {
1407       PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[k]);
1408       v += bs2;
1409     }
1410     PetscKernel_w_gets_A_times_v(bs,ls,aa+bs2*adiag[i],t+i*bs); /* *inv(diagonal[i]) */
1411     ierr = PetscArraycpy(x+i*bs,t+i*bs,bs);CHKERRQ(ierr);
1412   }
1413 
1414   ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr);
1415   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
1416   ierr = PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);CHKERRQ(ierr);
1417   PetscFunctionReturn(0);
1418 }
1419 
MatSolve_SeqBAIJ_N(Mat A,Vec bb,Vec xx)1420 PetscErrorCode MatSolve_SeqBAIJ_N(Mat A,Vec bb,Vec xx)
1421 {
1422   Mat_SeqBAIJ        *a   =(Mat_SeqBAIJ*)A->data;
1423   IS                 iscol=a->col,isrow=a->row;
1424   PetscErrorCode     ierr;
1425   const PetscInt     *r,*c,*rout,*cout,*ai=a->i,*aj=a->j,*adiag=a->diag,*vi;
1426   PetscInt           i,m,n=a->mbs;
1427   PetscInt           nz,bs=A->rmap->bs,bs2=a->bs2;
1428   const MatScalar    *aa=a->a,*v;
1429   PetscScalar        *x,*s,*t,*ls;
1430   const PetscScalar  *b;
1431 
1432   PetscFunctionBegin;
1433   ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr);
1434   ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
1435   t    = a->solve_work;
1436 
1437   ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout;
1438   ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout;
1439 
1440   /* forward solve the lower triangular */
1441   ierr = PetscArraycpy(t,b+bs*r[0],bs);CHKERRQ(ierr);
1442   for (i=1; i<n; i++) {
1443     v    = aa + bs2*ai[i];
1444     vi   = aj + ai[i];
1445     nz   = ai[i+1] - ai[i];
1446     s    = t + bs*i;
1447     ierr = PetscArraycpy(s,b+bs*r[i],bs);CHKERRQ(ierr);
1448     for (m=0; m<nz; m++) {
1449       PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[m]);
1450       v += bs2;
1451     }
1452   }
1453 
1454   /* backward solve the upper triangular */
1455   ls = a->solve_work + A->cmap->n;
1456   for (i=n-1; i>=0; i--) {
1457     v    = aa + bs2*(adiag[i+1]+1);
1458     vi   = aj + adiag[i+1]+1;
1459     nz   = adiag[i] - adiag[i+1] - 1;
1460     ierr = PetscArraycpy(ls,t+i*bs,bs);CHKERRQ(ierr);
1461     for (m=0; m<nz; m++) {
1462       PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[m]);
1463       v += bs2;
1464     }
1465     PetscKernel_w_gets_A_times_v(bs,ls,v,t+i*bs); /* *inv(diagonal[i]) */
1466     ierr = PetscArraycpy(x + bs*c[i],t+i*bs,bs);CHKERRQ(ierr);
1467   }
1468   ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr);
1469   ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr);
1470   ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr);
1471   ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
1472   ierr = PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);CHKERRQ(ierr);
1473   PetscFunctionReturn(0);
1474 }
1475 
1476 /*
1477     For each block in an block array saves the largest absolute value in the block into another array
1478 */
MatBlockAbs_private(PetscInt nbs,PetscInt bs2,PetscScalar * blockarray,PetscReal * absarray)1479 static PetscErrorCode MatBlockAbs_private(PetscInt nbs,PetscInt bs2,PetscScalar *blockarray,PetscReal *absarray)
1480 {
1481   PetscErrorCode ierr;
1482   PetscInt       i,j;
1483 
1484   PetscFunctionBegin;
1485   ierr = PetscArrayzero(absarray,nbs+1);CHKERRQ(ierr);
1486   for (i=0; i<nbs; i++) {
1487     for (j=0; j<bs2; j++) {
1488       if (absarray[i] < PetscAbsScalar(blockarray[i*nbs+j])) absarray[i] = PetscAbsScalar(blockarray[i*nbs+j]);
1489     }
1490   }
1491   PetscFunctionReturn(0);
1492 }
1493 
1494 /*
1495      This needs to be renamed and called by the regular MatILUFactor_SeqBAIJ when drop tolerance is used
1496 */
MatILUDTFactor_SeqBAIJ(Mat A,IS isrow,IS iscol,const MatFactorInfo * info,Mat * fact)1497 PetscErrorCode MatILUDTFactor_SeqBAIJ(Mat A,IS isrow,IS iscol,const MatFactorInfo *info,Mat *fact)
1498 {
1499   Mat            B = *fact;
1500   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data,*b;
1501   IS             isicol;
1502   PetscErrorCode ierr;
1503   const PetscInt *r,*ic;
1504   PetscInt       i,mbs=a->mbs,bs=A->rmap->bs,bs2=a->bs2,*ai=a->i,*aj=a->j,*ajtmp,*adiag;
1505   PetscInt       *bi,*bj,*bdiag;
1506 
1507   PetscInt  row,nzi,nzi_bl,nzi_bu,*im,dtcount,nzi_al,nzi_au;
1508   PetscInt  nlnk,*lnk;
1509   PetscBT   lnkbt;
1510   PetscBool row_identity,icol_identity;
1511   MatScalar *aatmp,*pv,*batmp,*ba,*rtmp,*pc,*multiplier,*vtmp;
1512   PetscInt  j,nz,*pj,*bjtmp,k,ncut,*jtmp;
1513 
1514   PetscReal dt=info->dt;          /* shift=info->shiftamount; */
1515   PetscInt  nnz_max;
1516   PetscBool missing;
1517   PetscReal *vtmp_abs;
1518   MatScalar *v_work;
1519   PetscInt  *v_pivots;
1520   PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE;
1521 
1522   PetscFunctionBegin;
1523   /* ------- symbolic factorization, can be reused ---------*/
1524   ierr = MatMissingDiagonal(A,&missing,&i);CHKERRQ(ierr);
1525   if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);
1526   adiag=a->diag;
1527 
1528   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
1529 
1530   /* bdiag is location of diagonal in factor */
1531   ierr = PetscMalloc1(mbs+1,&bdiag);CHKERRQ(ierr);
1532 
1533   /* allocate row pointers bi */
1534   ierr = PetscMalloc1(2*mbs+2,&bi);CHKERRQ(ierr);
1535 
1536   /* allocate bj and ba; max num of nonzero entries is (ai[n]+2*n*dtcount+2) */
1537   dtcount = (PetscInt)info->dtcount;
1538   if (dtcount > mbs-1) dtcount = mbs-1;
1539   nnz_max = ai[mbs]+2*mbs*dtcount +2;
1540   /* printf("MatILUDTFactor_SeqBAIJ, bs %d, ai[mbs] %d, nnz_max  %d, dtcount %d\n",bs,ai[mbs],nnz_max,dtcount); */
1541   ierr    = PetscMalloc1(nnz_max,&bj);CHKERRQ(ierr);
1542   nnz_max = nnz_max*bs2;
1543   ierr    = PetscMalloc1(nnz_max,&ba);CHKERRQ(ierr);
1544 
1545   /* put together the new matrix */
1546   ierr = MatSeqBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr);
1547   ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);CHKERRQ(ierr);
1548 
1549   b               = (Mat_SeqBAIJ*)(B)->data;
1550   b->free_a       = PETSC_TRUE;
1551   b->free_ij      = PETSC_TRUE;
1552   b->singlemalloc = PETSC_FALSE;
1553 
1554   b->a    = ba;
1555   b->j    = bj;
1556   b->i    = bi;
1557   b->diag = bdiag;
1558   b->ilen = NULL;
1559   b->imax = NULL;
1560   b->row  = isrow;
1561   b->col  = iscol;
1562 
1563   ierr = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
1564   ierr = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
1565 
1566   b->icol  = isicol;
1567   ierr     = PetscMalloc1(bs*(mbs+1),&b->solve_work);CHKERRQ(ierr);
1568   ierr     = PetscLogObjectMemory((PetscObject)B,nnz_max*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr);
1569   b->maxnz = nnz_max/bs2;
1570 
1571   (B)->factortype            = MAT_FACTOR_ILUDT;
1572   (B)->info.factor_mallocs   = 0;
1573   (B)->info.fill_ratio_given = ((PetscReal)nnz_max)/((PetscReal)(ai[mbs]*bs2));
1574   /* ------- end of symbolic factorization ---------*/
1575   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
1576   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
1577 
1578   /* linked list for storing column indices of the active row */
1579   nlnk = mbs + 1;
1580   ierr = PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr);
1581 
1582   /* im: used by PetscLLAddSortedLU(); jtmp: working array for column indices of active row */
1583   ierr = PetscMalloc2(mbs,&im,mbs,&jtmp);CHKERRQ(ierr);
1584   /* rtmp, vtmp: working arrays for sparse and contiguous row entries of active row */
1585   ierr = PetscMalloc2(mbs*bs2,&rtmp,mbs*bs2,&vtmp);CHKERRQ(ierr);
1586   ierr = PetscMalloc1(mbs+1,&vtmp_abs);CHKERRQ(ierr);
1587   ierr = PetscMalloc3(bs,&v_work,bs2,&multiplier,bs,&v_pivots);CHKERRQ(ierr);
1588 
1589   allowzeropivot = PetscNot(A->erroriffailure);
1590   bi[0]       = 0;
1591   bdiag[0]    = (nnz_max/bs2)-1; /* location of diagonal in factor B */
1592   bi[2*mbs+1] = bdiag[0]+1; /* endof bj and ba array */
1593   for (i=0; i<mbs; i++) {
1594     /* copy initial fill into linked list */
1595     nzi = ai[r[i]+1] - ai[r[i]];
1596     if (!nzi) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",r[i],i);
1597     nzi_al = adiag[r[i]] - ai[r[i]];
1598     nzi_au = ai[r[i]+1] - adiag[r[i]] -1;
1599 
1600     /* load in initial unfactored row */
1601     ajtmp = aj + ai[r[i]];
1602     ierr  = PetscLLAddPerm(nzi,ajtmp,ic,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr);
1603     ierr  = PetscArrayzero(rtmp,mbs*bs2);CHKERRQ(ierr);
1604     aatmp = a->a + bs2*ai[r[i]];
1605     for (j=0; j<nzi; j++) {
1606       ierr = PetscArraycpy(rtmp+bs2*ic[ajtmp[j]],aatmp+bs2*j,bs2);CHKERRQ(ierr);
1607     }
1608 
1609     /* add pivot rows into linked list */
1610     row = lnk[mbs];
1611     while (row < i) {
1612       nzi_bl = bi[row+1] - bi[row] + 1;
1613       bjtmp  = bj + bdiag[row+1]+1; /* points to 1st column next to the diagonal in U */
1614       ierr   = PetscLLAddSortedLU(bjtmp,row,nlnk,lnk,lnkbt,i,nzi_bl,im);CHKERRQ(ierr);
1615       nzi   += nlnk;
1616       row    = lnk[row];
1617     }
1618 
1619     /* copy data from lnk into jtmp, then initialize lnk */
1620     ierr = PetscLLClean(mbs,mbs,nzi,lnk,jtmp,lnkbt);CHKERRQ(ierr);
1621 
1622     /* numerical factorization */
1623     bjtmp = jtmp;
1624     row   = *bjtmp++; /* 1st pivot row */
1625 
1626     while  (row < i) {
1627       pc = rtmp + bs2*row;
1628       pv = ba + bs2*bdiag[row]; /* inv(diag) of the pivot row */
1629       PetscKernel_A_gets_A_times_B(bs,pc,pv,multiplier); /* pc= multiplier = pc*inv(diag[row]) */
1630       ierr = MatBlockAbs_private(1,bs2,pc,vtmp_abs);CHKERRQ(ierr);
1631       if (vtmp_abs[0] > dt) { /* apply tolerance dropping rule */
1632         pj = bj + bdiag[row+1] + 1;         /* point to 1st entry of U(row,:) */
1633         pv = ba + bs2*(bdiag[row+1] + 1);
1634         nz = bdiag[row] - bdiag[row+1] - 1;         /* num of entries in U(row,:), excluding diagonal */
1635         for (j=0; j<nz; j++) {
1636           PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j);
1637         }
1638         /* ierr = PetscLogFlops(bslog*(nz+1.0)-bs);CHKERRQ(ierr); */
1639       }
1640       row = *bjtmp++;
1641     }
1642 
1643     /* copy sparse rtmp into contiguous vtmp; separate L and U part */
1644     nzi_bl = 0; j = 0;
1645     while (jtmp[j] < i) { /* L-part. Note: jtmp is sorted */
1646       ierr = PetscArraycpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2);CHKERRQ(ierr);
1647       nzi_bl++; j++;
1648     }
1649     nzi_bu = nzi - nzi_bl -1;
1650 
1651     while (j < nzi) { /* U-part */
1652       ierr = PetscArraycpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2);CHKERRQ(ierr);
1653       j++;
1654     }
1655 
1656     ierr = MatBlockAbs_private(nzi,bs2,vtmp,vtmp_abs);CHKERRQ(ierr);
1657 
1658     bjtmp = bj + bi[i];
1659     batmp = ba + bs2*bi[i];
1660     /* apply level dropping rule to L part */
1661     ncut = nzi_al + dtcount;
1662     if (ncut < nzi_bl) {
1663       ierr = PetscSortSplitReal(ncut,nzi_bl,vtmp_abs,jtmp);CHKERRQ(ierr);
1664       ierr = PetscSortIntWithScalarArray(ncut,jtmp,vtmp);CHKERRQ(ierr);
1665     } else {
1666       ncut = nzi_bl;
1667     }
1668     for (j=0; j<ncut; j++) {
1669       bjtmp[j] = jtmp[j];
1670       ierr     = PetscArraycpy(batmp+bs2*j,rtmp+bs2*bjtmp[j],bs2);CHKERRQ(ierr);
1671     }
1672     bi[i+1] = bi[i] + ncut;
1673     nzi     = ncut + 1;
1674 
1675     /* apply level dropping rule to U part */
1676     ncut = nzi_au + dtcount;
1677     if (ncut < nzi_bu) {
1678       ierr = PetscSortSplitReal(ncut,nzi_bu,vtmp_abs+nzi_bl+1,jtmp+nzi_bl+1);CHKERRQ(ierr);
1679       ierr = PetscSortIntWithScalarArray(ncut,jtmp+nzi_bl+1,vtmp+nzi_bl+1);CHKERRQ(ierr);
1680     } else {
1681       ncut = nzi_bu;
1682     }
1683     nzi += ncut;
1684 
1685     /* mark bdiagonal */
1686     bdiag[i+1]    = bdiag[i] - (ncut + 1);
1687     bi[2*mbs - i] = bi[2*mbs - i +1] - (ncut + 1);
1688 
1689     bjtmp  = bj + bdiag[i];
1690     batmp  = ba + bs2*bdiag[i];
1691     ierr   = PetscArraycpy(batmp,rtmp+bs2*i,bs2);CHKERRQ(ierr);
1692     *bjtmp = i;
1693 
1694     bjtmp = bj + bdiag[i+1]+1;
1695     batmp = ba + (bdiag[i+1]+1)*bs2;
1696 
1697     for (k=0; k<ncut; k++) {
1698       bjtmp[k] = jtmp[nzi_bl+1+k];
1699       ierr     = PetscArraycpy(batmp+bs2*k,rtmp+bs2*bjtmp[k],bs2);CHKERRQ(ierr);
1700     }
1701 
1702     im[i] = nzi; /* used by PetscLLAddSortedLU() */
1703 
1704     /* invert diagonal block for simplier triangular solves - add shift??? */
1705     batmp = ba + bs2*bdiag[i];
1706 
1707     ierr  = PetscKernel_A_gets_inverse_A(bs,batmp,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr);
1708     if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1709   } /* for (i=0; i<mbs; i++) */
1710   ierr = PetscFree3(v_work,multiplier,v_pivots);CHKERRQ(ierr);
1711 
1712   /* printf("end of L %d, beginning of U %d\n",bi[mbs],bdiag[mbs]); */
1713   if (bi[mbs] >= bdiag[mbs]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"end of L array %d cannot >= the beginning of U array %d",bi[mbs],bdiag[mbs]);
1714 
1715   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
1716   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
1717 
1718   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
1719 
1720   ierr = PetscFree2(im,jtmp);CHKERRQ(ierr);
1721   ierr = PetscFree2(rtmp,vtmp);CHKERRQ(ierr);
1722 
1723   ierr     = PetscLogFlops(bs2*B->cmap->n);CHKERRQ(ierr);
1724   b->maxnz = b->nz = bi[mbs] + bdiag[0] - bdiag[mbs];
1725 
1726   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
1727   ierr = ISIdentity(isicol,&icol_identity);CHKERRQ(ierr);
1728   if (row_identity && icol_identity) {
1729     B->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
1730   } else {
1731     B->ops->solve = MatSolve_SeqBAIJ_N;
1732   }
1733 
1734   B->ops->solveadd          = NULL;
1735   B->ops->solvetranspose    = NULL;
1736   B->ops->solvetransposeadd = NULL;
1737   B->ops->matsolve          = NULL;
1738   B->assembled              = PETSC_TRUE;
1739   B->preallocated           = PETSC_TRUE;
1740   PetscFunctionReturn(0);
1741 }
1742