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