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