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