Actual source code: baijfact.c
petsc-3.4.5 2014-06-29
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,MatScalar,&rtmp,bs2,MatScalar,&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,MatScalar,&rtmp,bs2,MatScalar,&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->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering;
233: C->assembled = PETSC_TRUE;
235: PetscLogFlops(1.333333333333*2*2*2*n); /* from inverting diagonal blocks */
236: return(0);
237: }
241: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_inplace(Mat B,Mat A,const MatFactorInfo *info)
242: {
243: Mat C = B;
244: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data;
245: IS isrow = b->row,isicol = b->icol;
247: const PetscInt *r,*ic;
248: PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j;
249: PetscInt *ajtmpold,*ajtmp,nz,row;
250: PetscInt *diag_offset=b->diag,idx,*ai=a->i,*aj=a->j,*pj;
251: MatScalar *pv,*v,*rtmp,m1,m2,m3,m4,*pc,*w,*x,x1,x2,x3,x4;
252: MatScalar p1,p2,p3,p4;
253: MatScalar *ba = b->a,*aa = a->a;
254: PetscReal shift = info->shiftamount;
257: ISGetIndices(isrow,&r);
258: ISGetIndices(isicol,&ic);
259: PetscMalloc(4*(n+1)*sizeof(MatScalar),&rtmp);
261: for (i=0; i<n; i++) {
262: nz = bi[i+1] - bi[i];
263: ajtmp = bj + bi[i];
264: for (j=0; j<nz; j++) {
265: x = rtmp+4*ajtmp[j]; x[0] = x[1] = x[2] = x[3] = 0.0;
266: }
267: /* load in initial (unfactored row) */
268: idx = r[i];
269: nz = ai[idx+1] - ai[idx];
270: ajtmpold = aj + ai[idx];
271: v = aa + 4*ai[idx];
272: for (j=0; j<nz; j++) {
273: x = rtmp+4*ic[ajtmpold[j]];
274: x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3];
275: v += 4;
276: }
277: row = *ajtmp++;
278: while (row < i) {
279: pc = rtmp + 4*row;
280: p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3];
281: if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) {
282: pv = ba + 4*diag_offset[row];
283: pj = bj + diag_offset[row] + 1;
284: x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3];
285: pc[0] = m1 = p1*x1 + p3*x2;
286: pc[1] = m2 = p2*x1 + p4*x2;
287: pc[2] = m3 = p1*x3 + p3*x4;
288: pc[3] = m4 = p2*x3 + p4*x4;
289: nz = bi[row+1] - diag_offset[row] - 1;
290: pv += 4;
291: for (j=0; j<nz; j++) {
292: x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3];
293: x = rtmp + 4*pj[j];
294: x[0] -= m1*x1 + m3*x2;
295: x[1] -= m2*x1 + m4*x2;
296: x[2] -= m1*x3 + m3*x4;
297: x[3] -= m2*x3 + m4*x4;
298: pv += 4;
299: }
300: PetscLogFlops(16.0*nz+12.0);
301: }
302: row = *ajtmp++;
303: }
304: /* finished row so stick it into b->a */
305: pv = ba + 4*bi[i];
306: pj = bj + bi[i];
307: nz = bi[i+1] - bi[i];
308: for (j=0; j<nz; j++) {
309: x = rtmp+4*pj[j];
310: pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3];
311: pv += 4;
312: }
313: /* invert diagonal block */
314: w = ba + 4*diag_offset[i];
315: PetscKernel_A_gets_inverse_A_2(w,shift);
316: }
318: PetscFree(rtmp);
319: ISRestoreIndices(isicol,&ic);
320: ISRestoreIndices(isrow,&r);
322: C->ops->solve = MatSolve_SeqBAIJ_2_inplace;
323: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_inplace;
324: C->assembled = PETSC_TRUE;
326: PetscLogFlops(1.333333333333*8*b->mbs); /* from inverting diagonal blocks */
327: return(0);
328: }
329: /*
330: Version for when blocks are 2 by 2 Using natural ordering
331: */
334: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo *info)
335: {
336: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data;
338: PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j;
339: PetscInt *ajtmpold,*ajtmp,nz,row;
340: PetscInt *diag_offset = b->diag,*ai=a->i,*aj=a->j,*pj;
341: MatScalar *pv,*v,*rtmp,*pc,*w,*x;
342: MatScalar p1,p2,p3,p4,m1,m2,m3,m4,x1,x2,x3,x4;
343: MatScalar *ba = b->a,*aa = a->a;
344: PetscReal shift = info->shiftamount;
347: PetscMalloc(4*(n+1)*sizeof(MatScalar),&rtmp);
348: for (i=0; i<n; i++) {
349: nz = bi[i+1] - bi[i];
350: ajtmp = bj + bi[i];
351: for (j=0; j<nz; j++) {
352: x = rtmp+4*ajtmp[j];
353: x[0] = x[1] = x[2] = x[3] = 0.0;
354: }
355: /* load in initial (unfactored row) */
356: nz = ai[i+1] - ai[i];
357: ajtmpold = aj + ai[i];
358: v = aa + 4*ai[i];
359: for (j=0; j<nz; j++) {
360: x = rtmp+4*ajtmpold[j];
361: x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3];
362: v += 4;
363: }
364: row = *ajtmp++;
365: while (row < i) {
366: pc = rtmp + 4*row;
367: p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3];
368: if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) {
369: pv = ba + 4*diag_offset[row];
370: pj = bj + diag_offset[row] + 1;
371: x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3];
372: pc[0] = m1 = p1*x1 + p3*x2;
373: pc[1] = m2 = p2*x1 + p4*x2;
374: pc[2] = m3 = p1*x3 + p3*x4;
375: pc[3] = m4 = p2*x3 + p4*x4;
376: nz = bi[row+1] - diag_offset[row] - 1;
377: pv += 4;
378: for (j=0; j<nz; j++) {
379: x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3];
380: x = rtmp + 4*pj[j];
381: x[0] -= m1*x1 + m3*x2;
382: x[1] -= m2*x1 + m4*x2;
383: x[2] -= m1*x3 + m3*x4;
384: x[3] -= m2*x3 + m4*x4;
385: pv += 4;
386: }
387: PetscLogFlops(16.0*nz+12.0);
388: }
389: row = *ajtmp++;
390: }
391: /* finished row so stick it into b->a */
392: pv = ba + 4*bi[i];
393: pj = bj + bi[i];
394: nz = bi[i+1] - bi[i];
395: for (j=0; j<nz; j++) {
396: x = rtmp+4*pj[j];
397: pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3];
398: /*
399: printf(" col %d:",pj[j]);
400: PetscInt j1;
401: for (j1=0; j1<4; j1++) printf(" %g,",*(pv+j1));
402: printf("\n");
403: */
404: pv += 4;
405: }
406: /* invert diagonal block */
407: w = ba + 4*diag_offset[i];
408: /*
409: printf(" \n%d -th: diag: ",i);
410: for (j=0; j<4; j++) {
411: printf(" %g,",w[j]);
412: }
413: printf("\n----------------------------\n");
414: */
415: PetscKernel_A_gets_inverse_A_2(w,shift);
416: }
418: PetscFree(rtmp);
420: C->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering_inplace;
421: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering_inplace;
422: C->assembled = PETSC_TRUE;
424: PetscLogFlops(1.333333333333*8*b->mbs); /* from inverting diagonal blocks */
425: return(0);
426: }
428: /* ----------------------------------------------------------- */
429: /*
430: Version for when blocks are 1 by 1.
431: */
434: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1(Mat B,Mat A,const MatFactorInfo *info)
435: {
436: Mat C =B;
437: Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data;
438: IS isrow = b->row,isicol = b->icol;
439: PetscErrorCode ierr;
440: const PetscInt *r,*ic,*ics;
441: const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bdiag=b->diag;
442: PetscInt i,j,k,nz,nzL,row,*pj;
443: const PetscInt *ajtmp,*bjtmp;
444: MatScalar *rtmp,*pc,multiplier,*pv;
445: const MatScalar *aa=a->a,*v;
446: PetscBool row_identity,col_identity;
447: FactorShiftCtx sctx;
448: const PetscInt *ddiag;
449: PetscReal rs;
450: MatScalar d;
453: /* MatPivotSetUp(): initialize shift context sctx */
454: PetscMemzero(&sctx,sizeof(FactorShiftCtx));
456: if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
457: ddiag = a->diag;
458: sctx.shift_top = info->zeropivot;
459: for (i=0; i<n; i++) {
460: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
461: d = (aa)[ddiag[i]];
462: rs = -PetscAbsScalar(d) - PetscRealPart(d);
463: v = aa+ai[i];
464: nz = ai[i+1] - ai[i];
465: for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
466: if (rs>sctx.shift_top) sctx.shift_top = rs;
467: }
468: sctx.shift_top *= 1.1;
469: sctx.nshift_max = 5;
470: sctx.shift_lo = 0.;
471: sctx.shift_hi = 1.;
472: }
474: ISGetIndices(isrow,&r);
475: ISGetIndices(isicol,&ic);
476: PetscMalloc((n+1)*sizeof(MatScalar),&rtmp);
477: ics = ic;
479: do {
480: sctx.newshift = PETSC_FALSE;
481: for (i=0; i<n; i++) {
482: /* zero rtmp */
483: /* L part */
484: nz = bi[i+1] - bi[i];
485: bjtmp = bj + bi[i];
486: for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;
488: /* U part */
489: nz = bdiag[i]-bdiag[i+1];
490: bjtmp = bj + bdiag[i+1]+1;
491: for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0;
493: /* load in initial (unfactored row) */
494: nz = ai[r[i]+1] - ai[r[i]];
495: ajtmp = aj + ai[r[i]];
496: v = aa + ai[r[i]];
497: for (j=0; j<nz; j++) rtmp[ics[ajtmp[j]]] = v[j];
499: /* ZeropivotApply() */
500: rtmp[i] += sctx.shift_amount; /* shift the diagonal of the matrix */
502: /* elimination */
503: bjtmp = bj + bi[i];
504: row = *bjtmp++;
505: nzL = bi[i+1] - bi[i];
506: for (k=0; k < nzL; k++) {
507: pc = rtmp + row;
508: if (*pc != (PetscScalar)0.0) {
509: pv = b->a + bdiag[row];
510: multiplier = *pc * (*pv);
511: *pc = multiplier;
513: pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
514: pv = b->a + bdiag[row+1]+1;
515: nz = bdiag[row]-bdiag[row+1]-1; /* num of entries in U(row,:) excluding diag */
516: for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
517: PetscLogFlops(2.0*nz);
518: }
519: row = *bjtmp++;
520: }
522: /* finished row so stick it into b->a */
523: rs = 0.0;
524: /* L part */
525: pv = b->a + bi[i];
526: pj = b->j + bi[i];
527: nz = bi[i+1] - bi[i];
528: for (j=0; j<nz; j++) {
529: pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
530: }
532: /* U part */
533: pv = b->a + bdiag[i+1]+1;
534: pj = b->j + bdiag[i+1]+1;
535: nz = bdiag[i] - bdiag[i+1]-1;
536: for (j=0; j<nz; j++) {
537: pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]);
538: }
540: sctx.rs = rs;
541: sctx.pv = rtmp[i];
542: MatPivotCheck(A,info,&sctx,i);
543: if (sctx.newshift) break; /* break for-loop */
544: rtmp[i] = sctx.pv; /* sctx.pv might be updated in the case of MAT_SHIFT_INBLOCKS */
546: /* Mark diagonal and invert diagonal for simplier triangular solves */
547: pv = b->a + bdiag[i];
548: *pv = (PetscScalar)1.0/rtmp[i];
550: } /* endof for (i=0; i<n; i++) { */
552: /* MatPivotRefine() */
553: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
554: /*
555: * if no shift in this attempt & shifting & started shifting & can refine,
556: * then try lower shift
557: */
558: sctx.shift_hi = sctx.shift_fraction;
559: sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
560: sctx.shift_amount = sctx.shift_fraction * sctx.shift_top;
561: sctx.newshift = PETSC_TRUE;
562: sctx.nshift++;
563: }
564: } while (sctx.newshift);
566: PetscFree(rtmp);
567: ISRestoreIndices(isicol,&ic);
568: ISRestoreIndices(isrow,&r);
570: ISIdentity(isrow,&row_identity);
571: ISIdentity(isicol,&col_identity);
572: if (row_identity && col_identity) {
573: C->ops->solve = MatSolve_SeqBAIJ_1_NaturalOrdering;
574: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering;
575: } else {
576: C->ops->solve = MatSolve_SeqBAIJ_1;
577: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1;
578: }
579: C->assembled = PETSC_TRUE;
580: PetscLogFlops(C->cmap->n);
582: /* MatShiftView(A,info,&sctx) */
583: if (sctx.nshift) {
584: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
585: PetscInfo4(A,"number of shift_pd tries %D, shift_amount %G, diagonal shifted up by %e fraction top_value %e\n",sctx.nshift,sctx.shift_amount,sctx.shift_fraction,sctx.shift_top);
586: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
587: PetscInfo2(A,"number of shift_nz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);
588: } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
589: PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %G\n",sctx.nshift,info->shiftamount);
590: }
591: }
592: return(0);
593: }
597: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1_inplace(Mat C,Mat A,const MatFactorInfo *info)
598: {
599: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data;
600: IS isrow = b->row,isicol = b->icol;
602: const PetscInt *r,*ic;
603: PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j;
604: PetscInt *ajtmpold,*ajtmp,nz,row,*ai = a->i,*aj = a->j;
605: PetscInt *diag_offset = b->diag,diag,*pj;
606: MatScalar *pv,*v,*rtmp,multiplier,*pc;
607: MatScalar *ba = b->a,*aa = a->a;
608: PetscBool row_identity, col_identity;
611: ISGetIndices(isrow,&r);
612: ISGetIndices(isicol,&ic);
613: PetscMalloc((n+1)*sizeof(MatScalar),&rtmp);
615: for (i=0; i<n; i++) {
616: nz = bi[i+1] - bi[i];
617: ajtmp = bj + bi[i];
618: for (j=0; j<nz; j++) rtmp[ajtmp[j]] = 0.0;
620: /* load in initial (unfactored row) */
621: nz = ai[r[i]+1] - ai[r[i]];
622: ajtmpold = aj + ai[r[i]];
623: v = aa + ai[r[i]];
624: for (j=0; j<nz; j++) rtmp[ic[ajtmpold[j]]] = v[j];
626: row = *ajtmp++;
627: while (row < i) {
628: pc = rtmp + row;
629: if (*pc != 0.0) {
630: pv = ba + diag_offset[row];
631: pj = bj + diag_offset[row] + 1;
632: multiplier = *pc * *pv++;
633: *pc = multiplier;
634: nz = bi[row+1] - diag_offset[row] - 1;
635: for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
636: PetscLogFlops(1.0+2.0*nz);
637: }
638: row = *ajtmp++;
639: }
640: /* finished row so stick it into b->a */
641: pv = ba + bi[i];
642: pj = bj + bi[i];
643: nz = bi[i+1] - bi[i];
644: for (j=0; j<nz; j++) pv[j] = rtmp[pj[j]];
645: diag = diag_offset[i] - bi[i];
646: /* check pivot entry for current row */
647: 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);
648: pv[diag] = 1.0/pv[diag];
649: }
651: PetscFree(rtmp);
652: ISRestoreIndices(isicol,&ic);
653: ISRestoreIndices(isrow,&r);
654: ISIdentity(isrow,&row_identity);
655: ISIdentity(isicol,&col_identity);
656: if (row_identity && col_identity) {
657: C->ops->solve = MatSolve_SeqBAIJ_1_NaturalOrdering_inplace;
658: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering_inplace;
659: } else {
660: C->ops->solve = MatSolve_SeqBAIJ_1_inplace;
661: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_inplace;
662: }
663: C->assembled = PETSC_TRUE;
664: PetscLogFlops(C->cmap->n);
665: return(0);
666: }
670: PETSC_EXTERN PetscErrorCode MatGetFactor_seqbaij_petsc(Mat A,MatFactorType ftype,Mat *B)
671: {
672: PetscInt n = A->rmap->n;
676: #if defined(PETSC_USE_COMPLEX)
677: if (A->hermitian && (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported");
678: #endif
679: MatCreate(PetscObjectComm((PetscObject)A),B);
680: MatSetSizes(*B,n,n,n,n);
681: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU || ftype == MAT_FACTOR_ILUDT) {
682: MatSetType(*B,MATSEQBAIJ);
684: (*B)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqBAIJ;
685: (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqBAIJ;
686: } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
687: MatSetType(*B,MATSEQSBAIJ);
688: MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);
690: (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqBAIJ;
691: (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqBAIJ;
692: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
693: (*B)->factortype = ftype;
694: return(0);
695: }
699: PetscErrorCode MatGetFactorAvailable_seqbaij_petsc(Mat A,MatFactorType ftype,PetscBool *flg)
700: {
702: *flg = PETSC_TRUE;
703: return(0);
704: }
706: /* ----------------------------------------------------------- */
709: PetscErrorCode MatLUFactor_SeqBAIJ(Mat A,IS row,IS col,const MatFactorInfo *info)
710: {
712: Mat C;
715: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&C);
716: MatLUFactorSymbolic(C,A,row,col,info);
717: MatLUFactorNumeric(C,A,info);
719: A->ops->solve = C->ops->solve;
720: A->ops->solvetranspose = C->ops->solvetranspose;
722: MatHeaderMerge(A,C);
723: PetscLogObjectParent(A,((Mat_SeqBAIJ*)(A->data))->icol);
724: return(0);
725: }
727: #include <../src/mat/impls/sbaij/seq/sbaij.h>
730: PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N(Mat C,Mat A,const MatFactorInfo *info)
731: {
733: Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data;
734: Mat_SeqSBAIJ *b=(Mat_SeqSBAIJ*)C->data;
735: IS ip=b->row;
736: const PetscInt *rip;
737: PetscInt i,j,mbs=a->mbs,bs=A->rmap->bs,*bi=b->i,*bj=b->j,*bcol;
738: PetscInt *ai=a->i,*aj=a->j;
739: PetscInt k,jmin,jmax,*jl,*il,col,nexti,ili,nz;
740: MatScalar *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
741: PetscReal rs;
742: FactorShiftCtx sctx;
745: if (bs > 1) { /* convert A to a SBAIJ matrix and apply Cholesky factorization from it */
746: if (!a->sbaijMat) {
747: MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
748: }
749: (a->sbaijMat)->ops->choleskyfactornumeric(C,a->sbaijMat,info);
750: MatDestroy(&a->sbaijMat);
751: return(0);
752: }
754: /* MatPivotSetUp(): initialize shift context sctx */
755: PetscMemzero(&sctx,sizeof(FactorShiftCtx));
757: ISGetIndices(ip,&rip);
758: PetscMalloc3(mbs,MatScalar,&rtmp,mbs,PetscInt,&il,mbs,PetscInt,&jl);
760: sctx.shift_amount = 0.;
761: sctx.nshift = 0;
762: do {
763: sctx.newshift = PETSC_FALSE;
764: for (i=0; i<mbs; i++) {
765: rtmp[i] = 0.0; jl[i] = mbs; il[0] = 0;
766: }
768: for (k = 0; k<mbs; k++) {
769: bval = ba + bi[k];
770: /* initialize k-th row by the perm[k]-th row of A */
771: jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
772: for (j = jmin; j < jmax; j++) {
773: col = rip[aj[j]];
774: if (col >= k) { /* only take upper triangular entry */
775: rtmp[col] = aa[j];
776: *bval++ = 0.0; /* for in-place factorization */
777: }
778: }
780: /* shift the diagonal of the matrix */
781: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
783: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
784: dk = rtmp[k];
785: i = jl[k]; /* first row to be added to k_th row */
787: while (i < k) {
788: nexti = jl[i]; /* next row to be added to k_th row */
790: /* compute multiplier, update diag(k) and U(i,k) */
791: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
792: uikdi = -ba[ili]*ba[bi[i]]; /* diagonal(k) */
793: dk += uikdi*ba[ili];
794: ba[ili] = uikdi; /* -U(i,k) */
796: /* add multiple of row i to k-th row */
797: jmin = ili + 1; jmax = bi[i+1];
798: if (jmin < jmax) {
799: for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
800: /* update il and jl for row i */
801: il[i] = jmin;
802: j = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
803: }
804: i = nexti;
805: }
807: /* shift the diagonals when zero pivot is detected */
808: /* compute rs=sum of abs(off-diagonal) */
809: rs = 0.0;
810: jmin = bi[k]+1;
811: nz = bi[k+1] - jmin;
812: if (nz) {
813: bcol = bj + jmin;
814: while (nz--) {
815: rs += PetscAbsScalar(rtmp[*bcol]);
816: bcol++;
817: }
818: }
820: sctx.rs = rs;
821: sctx.pv = dk;
822: MatPivotCheck(A,info,&sctx,k);
823: if (sctx.newshift) break;
824: dk = sctx.pv;
826: /* copy data into U(k,:) */
827: ba[bi[k]] = 1.0/dk; /* U(k,k) */
828: jmin = bi[k]+1; jmax = bi[k+1];
829: if (jmin < jmax) {
830: for (j=jmin; j<jmax; j++) {
831: col = bj[j]; ba[j] = rtmp[col]; rtmp[col] = 0.0;
832: }
833: /* add the k-th row into il and jl */
834: il[k] = jmin;
835: i = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
836: }
837: }
838: } while (sctx.newshift);
839: PetscFree3(rtmp,il,jl);
841: ISRestoreIndices(ip,&rip);
843: C->assembled = PETSC_TRUE;
844: C->preallocated = PETSC_TRUE;
846: PetscLogFlops(C->rmap->N);
847: if (sctx.nshift) {
848: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
849: PetscInfo2(A,"number of shiftpd tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);
850: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
851: PetscInfo2(A,"number of shiftnz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);
852: }
853: }
854: return(0);
855: }
859: PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info)
860: {
861: Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data;
862: Mat_SeqSBAIJ *b=(Mat_SeqSBAIJ*)C->data;
864: PetscInt i,j,am=a->mbs;
865: PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
866: PetscInt k,jmin,*jl,*il,nexti,ili,*acol,*bcol,nz;
867: MatScalar *rtmp,*ba=b->a,*aa=a->a,dk,uikdi,*aval,*bval;
868: PetscReal rs;
869: FactorShiftCtx sctx;
872: /* MatPivotSetUp(): initialize shift context sctx */
873: PetscMemzero(&sctx,sizeof(FactorShiftCtx));
875: PetscMalloc3(am,MatScalar,&rtmp,am,PetscInt,&il,am,PetscInt,&jl);
877: do {
878: sctx.newshift = PETSC_FALSE;
879: for (i=0; i<am; i++) {
880: rtmp[i] = 0.0; jl[i] = am; il[0] = 0;
881: }
883: for (k = 0; k<am; k++) {
884: /* initialize k-th row with elements nonzero in row perm(k) of A */
885: nz = ai[k+1] - ai[k];
886: acol = aj + ai[k];
887: aval = aa + ai[k];
888: bval = ba + bi[k];
889: while (nz--) {
890: if (*acol < k) { /* skip lower triangular entries */
891: acol++; aval++;
892: } else {
893: rtmp[*acol++] = *aval++;
894: *bval++ = 0.0; /* for in-place factorization */
895: }
896: }
898: /* shift the diagonal of the matrix */
899: if (sctx.nshift) rtmp[k] += sctx.shift_amount;
901: /* modify k-th row by adding in those rows i with U(i,k)!=0 */
902: dk = rtmp[k];
903: i = jl[k]; /* first row to be added to k_th row */
905: while (i < k) {
906: nexti = jl[i]; /* next row to be added to k_th row */
907: /* compute multiplier, update D(k) and U(i,k) */
908: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
909: uikdi = -ba[ili]*ba[bi[i]];
910: dk += uikdi*ba[ili];
911: ba[ili] = uikdi; /* -U(i,k) */
913: /* add multiple of row i to k-th row ... */
914: jmin = ili + 1;
915: nz = bi[i+1] - jmin;
916: if (nz > 0) {
917: bcol = bj + jmin;
918: bval = ba + jmin;
919: while (nz--) rtmp[*bcol++] += uikdi*(*bval++);
920: /* update il and jl for i-th row */
921: il[i] = jmin;
922: j = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
923: }
924: i = nexti;
925: }
927: /* shift the diagonals when zero pivot is detected */
928: /* compute rs=sum of abs(off-diagonal) */
929: rs = 0.0;
930: jmin = bi[k]+1;
931: nz = bi[k+1] - jmin;
932: if (nz) {
933: bcol = bj + jmin;
934: while (nz--) {
935: rs += PetscAbsScalar(rtmp[*bcol]);
936: bcol++;
937: }
938: }
940: sctx.rs = rs;
941: sctx.pv = dk;
942: MatPivotCheck(A,info,&sctx,k);
943: if (sctx.newshift) break; /* sctx.shift_amount is updated */
944: dk = sctx.pv;
946: /* copy data into U(k,:) */
947: ba[bi[k]] = 1.0/dk;
948: jmin = bi[k]+1;
949: nz = bi[k+1] - jmin;
950: if (nz) {
951: bcol = bj + jmin;
952: bval = ba + jmin;
953: while (nz--) {
954: *bval++ = rtmp[*bcol];
955: rtmp[*bcol++] = 0.0;
956: }
957: /* add k-th row into il and jl */
958: il[k] = jmin;
959: i = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
960: }
961: }
962: } while (sctx.newshift);
963: PetscFree3(rtmp,il,jl);
965: C->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
966: C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
967: C->assembled = PETSC_TRUE;
968: C->preallocated = PETSC_TRUE;
970: PetscLogFlops(C->rmap->N);
971: if (sctx.nshift) {
972: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
973: PetscInfo2(A,"number of shiftnz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);
974: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
975: PetscInfo2(A,"number of shiftpd tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);
976: }
977: }
978: return(0);
979: }
981: #include <petscbt.h>
982: #include <../src/mat/utils/freespace.h>
985: PetscErrorCode MatICCFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
986: {
987: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
988: Mat_SeqSBAIJ *b;
989: Mat B;
990: PetscErrorCode ierr;
991: PetscBool perm_identity;
992: PetscInt reallocs=0,i,*ai=a->i,*aj=a->j,am=a->mbs,bs=A->rmap->bs,*ui;
993: const PetscInt *rip;
994: PetscInt jmin,jmax,nzk,k,j,*jl,prow,*il,nextprow;
995: PetscInt nlnk,*lnk,*lnk_lvl=NULL,ncols,ncols_upper,*cols,*cols_lvl,*uj,**uj_ptr,**uj_lvl_ptr;
996: PetscReal fill =info->fill,levels=info->levels;
997: PetscFreeSpaceList free_space =NULL,current_space=NULL;
998: PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL;
999: PetscBT lnkbt;
1002: if (bs > 1) {
1003: if (!a->sbaijMat) {
1004: MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
1005: }
1006: (fact)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */
1008: MatICCFactorSymbolic(fact,a->sbaijMat,perm,info);
1009: return(0);
1010: }
1012: ISIdentity(perm,&perm_identity);
1013: ISGetIndices(perm,&rip);
1015: /* special case that simply copies fill pattern */
1016: if (!levels && perm_identity) {
1017: MatMarkDiagonal_SeqBAIJ(A);
1018: PetscMalloc((am+1)*sizeof(PetscInt),&ui);
1019: for (i=0; i<am; i++) ui[i] = ai[i+1] - a->diag[i]; /* ui: rowlengths - changes when !perm_identity */
1020: B = fact;
1021: MatSeqSBAIJSetPreallocation(B,1,0,ui);
1024: b = (Mat_SeqSBAIJ*)B->data;
1025: uj = b->j;
1026: for (i=0; i<am; i++) {
1027: aj = a->j + a->diag[i];
1028: for (j=0; j<ui[i]; j++) *uj++ = *aj++;
1029: b->ilen[i] = ui[i];
1030: }
1031: PetscFree(ui);
1033: B->factortype = MAT_FACTOR_NONE;
1035: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1036: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1037: B->factortype = MAT_FACTOR_ICC;
1039: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1040: return(0);
1041: }
1043: /* initialization */
1044: PetscMalloc((am+1)*sizeof(PetscInt),&ui);
1045: ui[0] = 0;
1046: PetscMalloc((2*am+1)*sizeof(PetscInt),&cols_lvl);
1048: /* jl: linked list for storing indices of the pivot rows
1049: il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
1050: PetscMalloc4(am,PetscInt*,&uj_ptr,am,PetscInt*,&uj_lvl_ptr,am,PetscInt,&il,am,PetscInt,&jl);
1051: for (i=0; i<am; i++) {
1052: jl[i] = am; il[i] = 0;
1053: }
1055: /* create and initialize a linked list for storing column indices of the active row k */
1056: nlnk = am + 1;
1057: PetscIncompleteLLCreate(am,am,nlnk,lnk,lnk_lvl,lnkbt);
1059: /* initial FreeSpace size is fill*(ai[am]+am)/2 */
1060: PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+am)/2),&free_space);
1062: current_space = free_space;
1064: PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+am)/2),&free_space_lvl);
1065: current_space_lvl = free_space_lvl;
1067: for (k=0; k<am; k++) { /* for each active row k */
1068: /* initialize lnk by the column indices of row rip[k] of A */
1069: nzk = 0;
1070: ncols = ai[rip[k]+1] - ai[rip[k]];
1071: ncols_upper = 0;
1072: cols = cols_lvl + am;
1073: for (j=0; j<ncols; j++) {
1074: i = rip[*(aj + ai[rip[k]] + j)];
1075: if (i >= k) { /* only take upper triangular entry */
1076: cols[ncols_upper] = i;
1077: cols_lvl[ncols_upper] = -1; /* initialize level for nonzero entries */
1078: ncols_upper++;
1079: }
1080: }
1081: PetscIncompleteLLAdd(ncols_upper,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);
1082: nzk += nlnk;
1084: /* update lnk by computing fill-in for each pivot row to be merged in */
1085: prow = jl[k]; /* 1st pivot row */
1087: while (prow < k) {
1088: nextprow = jl[prow];
1090: /* merge prow into k-th row */
1091: jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */
1092: jmax = ui[prow+1];
1093: ncols = jmax-jmin;
1094: i = jmin - ui[prow];
1095: cols = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
1096: for (j=0; j<ncols; j++) cols_lvl[j] = *(uj_lvl_ptr[prow] + i + j);
1097: PetscIncompleteLLAddSorted(ncols,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);
1098: nzk += nlnk;
1100: /* update il and jl for prow */
1101: if (jmin < jmax) {
1102: il[prow] = jmin;
1104: j = *cols; jl[prow] = jl[j]; jl[j] = prow;
1105: }
1106: prow = nextprow;
1107: }
1109: /* if free space is not available, make more free space */
1110: if (current_space->local_remaining<nzk) {
1111: i = am - k + 1; /* num of unfactored rows */
1112: i = PetscMin(i*nzk, i*(i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
1113: PetscFreeSpaceGet(i,¤t_space);
1114: PetscFreeSpaceGet(i,¤t_space_lvl);
1115: reallocs++;
1116: }
1118: /* copy data into free_space and free_space_lvl, then initialize lnk */
1119: PetscIncompleteLLClean(am,am,nzk,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);
1121: /* add the k-th row into il and jl */
1122: if (nzk-1 > 0) {
1123: i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
1124: jl[k] = jl[i]; jl[i] = k;
1125: il[k] = ui[k] + 1;
1126: }
1127: uj_ptr[k] = current_space->array;
1128: uj_lvl_ptr[k] = current_space_lvl->array;
1130: current_space->array += nzk;
1131: current_space->local_used += nzk;
1132: current_space->local_remaining -= nzk;
1134: current_space_lvl->array += nzk;
1135: current_space_lvl->local_used += nzk;
1136: current_space_lvl->local_remaining -= nzk;
1138: ui[k+1] = ui[k] + nzk;
1139: }
1141: ISRestoreIndices(perm,&rip);
1142: PetscFree4(uj_ptr,uj_lvl_ptr,il,jl);
1143: PetscFree(cols_lvl);
1145: /* copy free_space into uj and free free_space; set uj in new datastructure; */
1146: PetscMalloc((ui[am]+1)*sizeof(PetscInt),&uj);
1147: PetscFreeSpaceContiguous(&free_space,uj);
1148: PetscIncompleteLLDestroy(lnk,lnkbt);
1149: PetscFreeSpaceDestroy(free_space_lvl);
1151: /* put together the new matrix in MATSEQSBAIJ format */
1152: B = fact;
1153: MatSeqSBAIJSetPreallocation(B,1,MAT_SKIP_ALLOCATION,NULL);
1155: b = (Mat_SeqSBAIJ*)B->data;
1156: b->singlemalloc = PETSC_FALSE;
1157: b->free_a = PETSC_TRUE;
1158: b->free_ij = PETSC_TRUE;
1160: PetscMalloc((ui[am]+1)*sizeof(MatScalar),&b->a);
1162: b->j = uj;
1163: b->i = ui;
1164: b->diag = 0;
1165: b->ilen = 0;
1166: b->imax = 0;
1167: b->row = perm;
1168: b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
1170: PetscObjectReference((PetscObject)perm);
1172: b->icol = perm;
1174: PetscObjectReference((PetscObject)perm);
1175: PetscMalloc((am+1)*sizeof(PetscScalar),&b->solve_work);
1176: PetscLogObjectMemory(B,(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));
1178: b->maxnz = b->nz = ui[am];
1180: B->info.factor_mallocs = reallocs;
1181: B->info.fill_ratio_given = fill;
1182: if (ai[am] != 0.) {
1183: /* nonzeros in lower triangular part of A (includign diagonals)= (ai[am]+am)/2 */
1184: B->info.fill_ratio_needed = ((PetscReal)2*ui[am])/(ai[am]+am);
1185: } else {
1186: B->info.fill_ratio_needed = 0.0;
1187: }
1188: #if defined(PETSC_USE_INFO)
1189: if (ai[am] != 0) {
1190: PetscReal af = B->info.fill_ratio_needed;
1191: PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,fill,af);
1192: PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);
1193: PetscInfo1(A,"PCFactorSetFill(pc,%G) for best performance.\n",af);
1194: } else {
1195: PetscInfo(A,"Empty matrix.\n");
1196: }
1197: #endif
1198: if (perm_identity) {
1199: B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1200: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
1201: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
1202: } else {
1203: (fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
1204: }
1205: return(0);
1206: }
1210: PetscErrorCode MatCholeskyFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
1211: {
1212: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
1213: Mat_SeqSBAIJ *b;
1214: Mat B;
1215: PetscErrorCode ierr;
1216: PetscBool perm_identity;
1217: PetscReal fill = info->fill;
1218: const PetscInt *rip;
1219: PetscInt i,mbs=a->mbs,bs=A->rmap->bs,*ai=a->i,*aj=a->j,reallocs=0,prow;
1220: PetscInt *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
1221: PetscInt nlnk,*lnk,ncols,ncols_upper,*cols,*uj,**ui_ptr,*uj_ptr;
1222: PetscFreeSpaceList free_space=NULL,current_space=NULL;
1223: PetscBT lnkbt;
1226: if (bs > 1) { /* convert to seqsbaij */
1227: if (!a->sbaijMat) {
1228: MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
1229: }
1230: (fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */
1232: MatCholeskyFactorSymbolic(fact,a->sbaijMat,perm,info);
1233: return(0);
1234: }
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: PetscMalloc((mbs+1)*sizeof(PetscInt),&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,PetscInt*,&ui_ptr,mbs,PetscInt,&il,mbs,PetscInt,&jl,mbs,PetscInt,&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: PetscMalloc((ui[mbs]+1)*sizeof(PetscInt),&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: PetscMalloc((ui[mbs]+1)*sizeof(MatScalar),&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: PetscMalloc((mbs+1)*sizeof(PetscScalar),&b->solve_work);
1355: PetscLogObjectMemory(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,fill,af);
1370: PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);
1371: PetscInfo1(A,"PCFactorSetFill(pc,%G) for best performance.\n",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;
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: MatScalar *aa=a->a,*v;
1394: PetscScalar *x,*b,*s,*t,*ls;
1397: VecGetArray(bb,&b);
1398: VecGetArray(xx,&x);
1399: t = a->solve_work;
1401: /* forward solve the lower triangular */
1402: PetscMemcpy(t,b,bs*sizeof(PetscScalar)); /* copy 1st block of b to t */
1404: for (i=1; i<n; i++) {
1405: v = aa + bs2*ai[i];
1406: vi = aj + ai[i];
1407: nz = ai[i+1] - ai[i];
1408: s = t + bs*i;
1409: PetscMemcpy(s,b+bs*i,bs*sizeof(PetscScalar)); /* copy i_th block of b to t */
1410: for (k=0;k<nz;k++) {
1411: PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[k]);
1412: v += bs2;
1413: }
1414: }
1416: /* backward solve the upper triangular */
1417: ls = a->solve_work + A->cmap->n;
1418: for (i=n-1; i>=0; i--) {
1419: v = aa + bs2*(adiag[i+1]+1);
1420: vi = aj + adiag[i+1]+1;
1421: nz = adiag[i] - adiag[i+1]-1;
1422: PetscMemcpy(ls,t+i*bs,bs*sizeof(PetscScalar));
1423: for (k=0; k<nz; k++) {
1424: PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[k]);
1425: v += bs2;
1426: }
1427: PetscKernel_w_gets_A_times_v(bs,ls,aa+bs2*adiag[i],t+i*bs); /* *inv(diagonal[i]) */
1428: PetscMemcpy(x+i*bs,t+i*bs,bs*sizeof(PetscScalar));
1429: }
1431: VecRestoreArray(bb,&b);
1432: VecRestoreArray(xx,&x);
1433: PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);
1434: return(0);
1435: }
1439: PetscErrorCode MatSolve_SeqBAIJ_N(Mat A,Vec bb,Vec xx)
1440: {
1441: Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)A->data;
1442: IS iscol=a->col,isrow=a->row;
1444: const PetscInt *r,*c,*rout,*cout,*ai=a->i,*aj=a->j,*adiag=a->diag,*vi;
1445: PetscInt i,m,n=a->mbs;
1446: PetscInt nz,bs=A->rmap->bs,bs2=a->bs2;
1447: MatScalar *aa=a->a,*v;
1448: PetscScalar *x,*b,*s,*t,*ls;
1451: VecGetArray(bb,&b);
1452: VecGetArray(xx,&x);
1453: t = a->solve_work;
1455: ISGetIndices(isrow,&rout); r = rout;
1456: ISGetIndices(iscol,&cout); c = cout;
1458: /* forward solve the lower triangular */
1459: PetscMemcpy(t,b+bs*r[0],bs*sizeof(PetscScalar));
1460: for (i=1; i<n; i++) {
1461: v = aa + bs2*ai[i];
1462: vi = aj + ai[i];
1463: nz = ai[i+1] - ai[i];
1464: s = t + bs*i;
1465: PetscMemcpy(s,b+bs*r[i],bs*sizeof(PetscScalar));
1466: for (m=0; m<nz; m++) {
1467: PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[m]);
1468: v += bs2;
1469: }
1470: }
1472: /* backward solve the upper triangular */
1473: ls = a->solve_work + A->cmap->n;
1474: for (i=n-1; i>=0; i--) {
1475: v = aa + bs2*(adiag[i+1]+1);
1476: vi = aj + adiag[i+1]+1;
1477: nz = adiag[i] - adiag[i+1] - 1;
1478: PetscMemcpy(ls,t+i*bs,bs*sizeof(PetscScalar));
1479: for (m=0; m<nz; m++) {
1480: PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[m]);
1481: v += bs2;
1482: }
1483: PetscKernel_w_gets_A_times_v(bs,ls,v,t+i*bs); /* *inv(diagonal[i]) */
1484: PetscMemcpy(x + bs*c[i],t+i*bs,bs*sizeof(PetscScalar));
1485: }
1486: ISRestoreIndices(isrow,&rout);
1487: ISRestoreIndices(iscol,&cout);
1488: VecRestoreArray(bb,&b);
1489: VecRestoreArray(xx,&x);
1490: PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);
1491: return(0);
1492: }
1496: /*
1497: For each block in an block array saves the largest absolute value in the block into another array
1498: */
1499: static PetscErrorCode MatBlockAbs_private(PetscInt nbs,PetscInt bs2,PetscScalar *blockarray,PetscReal *absarray)
1500: {
1502: PetscInt i,j;
1505: PetscMemzero(absarray,(nbs+1)*sizeof(PetscReal));
1506: for (i=0; i<nbs; i++) {
1507: for (j=0; j<bs2; j++) {
1508: if (absarray[i] < PetscAbsScalar(blockarray[i*nbs+j])) absarray[i] = PetscAbsScalar(blockarray[i*nbs+j]);
1509: }
1510: }
1511: return(0);
1512: }
1516: /*
1517: This needs to be renamed and called by the regular MatILUFactor_SeqBAIJ when drop tolerance is used
1518: */
1519: PetscErrorCode MatILUDTFactor_SeqBAIJ(Mat A,IS isrow,IS iscol,const MatFactorInfo *info,Mat *fact)
1520: {
1521: Mat B = *fact;
1522: Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b;
1523: IS isicol;
1525: const PetscInt *r,*ic;
1526: PetscInt i,mbs=a->mbs,bs=A->rmap->bs,bs2=a->bs2,*ai=a->i,*aj=a->j,*ajtmp,*adiag;
1527: PetscInt *bi,*bj,*bdiag;
1529: PetscInt row,nzi,nzi_bl,nzi_bu,*im,dtcount,nzi_al,nzi_au;
1530: PetscInt nlnk,*lnk;
1531: PetscBT lnkbt;
1532: PetscBool row_identity,icol_identity;
1533: MatScalar *aatmp,*pv,*batmp,*ba,*rtmp,*pc,*multiplier,*vtmp;
1534: PetscInt j,nz,*pj,*bjtmp,k,ncut,*jtmp;
1536: PetscReal dt=info->dt; /* shift=info->shiftamount; */
1537: PetscInt nnz_max;
1538: PetscBool missing;
1539: PetscReal *vtmp_abs;
1540: MatScalar *v_work;
1541: PetscInt *v_pivots;
1544: /* ------- symbolic factorization, can be reused ---------*/
1545: MatMissingDiagonal(A,&missing,&i);
1546: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i);
1547: adiag=a->diag;
1549: ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
1551: /* bdiag is location of diagonal in factor */
1552: PetscMalloc((mbs+1)*sizeof(PetscInt),&bdiag);
1554: /* allocate row pointers bi */
1555: PetscMalloc((2*mbs+2)*sizeof(PetscInt),&bi);
1557: /* allocate bj and ba; max num of nonzero entries is (ai[n]+2*n*dtcount+2) */
1558: dtcount = (PetscInt)info->dtcount;
1559: if (dtcount > mbs-1) dtcount = mbs-1;
1560: nnz_max = ai[mbs]+2*mbs*dtcount +2;
1561: /* printf("MatILUDTFactor_SeqBAIJ, bs %d, ai[mbs] %d, nnz_max %d, dtcount %d\n",bs,ai[mbs],nnz_max,dtcount); */
1562: PetscMalloc(nnz_max*sizeof(PetscInt),&bj);
1563: nnz_max = nnz_max*bs2;
1564: PetscMalloc(nnz_max*sizeof(MatScalar),&ba);
1566: /* put together the new matrix */
1567: MatSeqBAIJSetPreallocation_SeqBAIJ(B,bs,MAT_SKIP_ALLOCATION,NULL);
1568: PetscLogObjectParent(B,isicol);
1570: b = (Mat_SeqBAIJ*)(B)->data;
1571: b->free_a = PETSC_TRUE;
1572: b->free_ij = PETSC_TRUE;
1573: b->singlemalloc = PETSC_FALSE;
1575: b->a = ba;
1576: b->j = bj;
1577: b->i = bi;
1578: b->diag = bdiag;
1579: b->ilen = 0;
1580: b->imax = 0;
1581: b->row = isrow;
1582: b->col = iscol;
1584: PetscObjectReference((PetscObject)isrow);
1585: PetscObjectReference((PetscObject)iscol);
1587: b->icol = isicol;
1588: PetscMalloc((bs*(mbs+1))*sizeof(PetscScalar),&b->solve_work);
1589: PetscLogObjectMemory(B,nnz_max*(sizeof(PetscInt)+sizeof(MatScalar)));
1590: b->maxnz = nnz_max/bs2;
1592: (B)->factortype = MAT_FACTOR_ILUDT;
1593: (B)->info.factor_mallocs = 0;
1594: (B)->info.fill_ratio_given = ((PetscReal)nnz_max)/((PetscReal)(ai[mbs]*bs2));
1595: /* ------- end of symbolic factorization ---------*/
1596: ISGetIndices(isrow,&r);
1597: ISGetIndices(isicol,&ic);
1599: /* linked list for storing column indices of the active row */
1600: nlnk = mbs + 1;
1601: PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);
1603: /* im: used by PetscLLAddSortedLU(); jtmp: working array for column indices of active row */
1604: PetscMalloc2(mbs,PetscInt,&im,mbs,PetscInt,&jtmp);
1605: /* rtmp, vtmp: working arrays for sparse and contiguous row entries of active row */
1606: PetscMalloc2(mbs*bs2,MatScalar,&rtmp,mbs*bs2,MatScalar,&vtmp);
1607: PetscMalloc((mbs+1)*sizeof(PetscReal),&vtmp_abs);
1608: PetscMalloc3(bs,MatScalar,&v_work,bs2,MatScalar,&multiplier,bs,PetscInt,&v_pivots);
1610: bi[0] = 0;
1611: bdiag[0] = (nnz_max/bs2)-1; /* location of diagonal in factor B */
1612: bi[2*mbs+1] = bdiag[0]+1; /* endof bj and ba array */
1613: for (i=0; i<mbs; i++) {
1614: /* copy initial fill into linked list */
1615: nzi = 0; /* nonzeros for active row i */
1616: nzi = ai[r[i]+1] - ai[r[i]];
1617: 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);
1618: nzi_al = adiag[r[i]] - ai[r[i]];
1619: nzi_au = ai[r[i]+1] - adiag[r[i]] -1;
1620: /* printf("row %d, nzi_al/au %d %d\n",i,nzi_al,nzi_au); */
1622: /* load in initial unfactored row */
1623: ajtmp = aj + ai[r[i]];
1624: PetscLLAddPerm(nzi,ajtmp,ic,mbs,nlnk,lnk,lnkbt);
1625: PetscMemzero(rtmp,mbs*bs2*sizeof(PetscScalar));
1626: aatmp = a->a + bs2*ai[r[i]];
1627: for (j=0; j<nzi; j++) {
1628: PetscMemcpy(rtmp+bs2*ic[ajtmp[j]],aatmp+bs2*j,bs2*sizeof(MatScalar));
1629: }
1631: /* add pivot rows into linked list */
1632: row = lnk[mbs];
1633: while (row < i) {
1634: nzi_bl = bi[row+1] - bi[row] + 1;
1635: bjtmp = bj + bdiag[row+1]+1; /* points to 1st column next to the diagonal in U */
1636: PetscLLAddSortedLU(bjtmp,row,nlnk,lnk,lnkbt,i,nzi_bl,im);
1637: nzi += nlnk;
1638: row = lnk[row];
1639: }
1641: /* copy data from lnk into jtmp, then initialize lnk */
1642: PetscLLClean(mbs,mbs,nzi,lnk,jtmp,lnkbt);
1644: /* numerical factorization */
1645: bjtmp = jtmp;
1646: row = *bjtmp++; /* 1st pivot row */
1648: while (row < i) {
1649: pc = rtmp + bs2*row;
1650: pv = ba + bs2*bdiag[row]; /* inv(diag) of the pivot row */
1651: PetscKernel_A_gets_A_times_B(bs,pc,pv,multiplier); /* pc= multiplier = pc*inv(diag[row]) */
1652: MatBlockAbs_private(1,bs2,pc,vtmp_abs);
1653: if (vtmp_abs[0] > dt) { /* apply tolerance dropping rule */
1654: pj = bj + bdiag[row+1] + 1; /* point to 1st entry of U(row,:) */
1655: pv = ba + bs2*(bdiag[row+1] + 1);
1656: nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries in U(row,:), excluding diagonal */
1657: for (j=0; j<nz; j++) {
1658: PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j);
1659: }
1660: /* PetscLogFlops(bslog*(nz+1.0)-bs); */
1661: }
1662: row = *bjtmp++;
1663: }
1665: /* copy sparse rtmp into contiguous vtmp; separate L and U part */
1666: nzi_bl = 0; j = 0;
1667: while (jtmp[j] < i) { /* L-part. Note: jtmp is sorted */
1668: PetscMemcpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2*sizeof(MatScalar));
1669: nzi_bl++; j++;
1670: }
1671: nzi_bu = nzi - nzi_bl -1;
1672: /* printf("nzi %d, nzi_bl %d, nzi_bu %d\n",nzi,nzi_bl,nzi_bu); */
1674: while (j < nzi) { /* U-part */
1675: PetscMemcpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2*sizeof(MatScalar));
1676: /*
1677: printf(" col %d: ",jtmp[j]);
1678: for (j1=0; j1<bs2; j1++) printf(" %g",*(vtmp+bs2*j+j1));
1679: printf(" \n");
1680: */
1681: j++;
1682: }
1684: MatBlockAbs_private(nzi,bs2,vtmp,vtmp_abs);
1685: /*
1686: printf(" row %d, nzi %d, vtmp_abs\n",i,nzi);
1687: for (j1=0; j1<nzi; j1++) printf(" (%d %g),",jtmp[j1],vtmp_abs[j1]);
1688: printf(" \n");
1689: */
1690: bjtmp = bj + bi[i];
1691: batmp = ba + bs2*bi[i];
1692: /* apply level dropping rule to L part */
1693: ncut = nzi_al + dtcount;
1694: if (ncut < nzi_bl) {
1695: PetscSortSplitReal(ncut,nzi_bl,vtmp_abs,jtmp);
1696: PetscSortIntWithScalarArray(ncut,jtmp,vtmp);
1697: } else {
1698: ncut = nzi_bl;
1699: }
1700: for (j=0; j<ncut; j++) {
1701: bjtmp[j] = jtmp[j];
1702: PetscMemcpy(batmp+bs2*j,rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));
1703: /*
1704: printf(" col %d: ",bjtmp[j]);
1705: for (j1=0; j1<bs2; j1++) printf(" %g,",*(batmp+bs2*j+j1));
1706: printf("\n");
1707: */
1708: }
1709: bi[i+1] = bi[i] + ncut;
1710: nzi = ncut + 1;
1712: /* apply level dropping rule to U part */
1713: ncut = nzi_au + dtcount;
1714: if (ncut < nzi_bu) {
1715: PetscSortSplitReal(ncut,nzi_bu,vtmp_abs+nzi_bl+1,jtmp+nzi_bl+1);
1716: PetscSortIntWithScalarArray(ncut,jtmp+nzi_bl+1,vtmp+nzi_bl+1);
1717: } else {
1718: ncut = nzi_bu;
1719: }
1720: nzi += ncut;
1722: /* mark bdiagonal */
1723: bdiag[i+1] = bdiag[i] - (ncut + 1);
1724: bi[2*mbs - i] = bi[2*mbs - i +1] - (ncut + 1);
1726: bjtmp = bj + bdiag[i];
1727: batmp = ba + bs2*bdiag[i];
1728: PetscMemcpy(batmp,rtmp+bs2*i,bs2*sizeof(MatScalar));
1729: *bjtmp = i;
1730: /*
1731: printf(" diag %d: ",*bjtmp);
1732: for (j=0; j<bs2; j++) {
1733: printf(" %g,",batmp[j]);
1734: }
1735: printf("\n");
1736: */
1737: bjtmp = bj + bdiag[i+1]+1;
1738: batmp = ba + (bdiag[i+1]+1)*bs2;
1740: for (k=0; k<ncut; k++) {
1741: bjtmp[k] = jtmp[nzi_bl+1+k];
1742: PetscMemcpy(batmp+bs2*k,rtmp+bs2*bjtmp[k],bs2*sizeof(MatScalar));
1743: /*
1744: printf(" col %d:",bjtmp[k]);
1745: for (j1=0; j1<bs2; j1++) printf(" %g,",*(batmp+bs2*k+j1));
1746: printf("\n");
1747: */
1748: }
1750: im[i] = nzi; /* used by PetscLLAddSortedLU() */
1752: /* invert diagonal block for simplier triangular solves - add shift??? */
1753: batmp = ba + bs2*bdiag[i];
1754: PetscKernel_A_gets_inverse_A(bs,batmp,v_pivots,v_work);
1755: } /* for (i=0; i<mbs; i++) */
1756: PetscFree3(v_work,multiplier,v_pivots);
1758: /* printf("end of L %d, beginning of U %d\n",bi[mbs],bdiag[mbs]); */
1759: 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]);
1761: ISRestoreIndices(isrow,&r);
1762: ISRestoreIndices(isicol,&ic);
1764: PetscLLDestroy(lnk,lnkbt);
1766: PetscFree2(im,jtmp);
1767: PetscFree2(rtmp,vtmp);
1769: PetscLogFlops(bs2*B->cmap->n);
1770: b->maxnz = b->nz = bi[mbs] + bdiag[0] - bdiag[mbs];
1772: ISIdentity(isrow,&row_identity);
1773: ISIdentity(isicol,&icol_identity);
1774: if (row_identity && icol_identity) {
1775: B->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
1776: } else {
1777: B->ops->solve = MatSolve_SeqBAIJ_N;
1778: }
1780: B->ops->solveadd = 0;
1781: B->ops->solvetranspose = 0;
1782: B->ops->solvetransposeadd = 0;
1783: B->ops->matsolve = 0;
1784: B->assembled = PETSC_TRUE;
1785: B->preallocated = PETSC_TRUE;
1786: return(0);
1787: }