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