Actual source code: baijfact2.c
petsc-3.7.3 2016-08-01
2: /*
3: Factorization code for BAIJ format.
4: */
6: #include <../src/mat/impls/baij/seq/baij.h>
7: #include <petsc/private/kernels/blockinvert.h>
8: #include <petscbt.h>
9: #include <../src/mat/utils/freespace.h>
11: /* ----------------------------------------------------------------*/
12: extern PetscErrorCode MatDuplicateNoCreate_SeqBAIJ(Mat,Mat,MatDuplicateOption,PetscBool);
16: /*
17: This is not much faster than MatLUFactorNumeric_SeqBAIJ_N() but the solve is faster at least sometimes
18: */
19: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_15_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info)
20: {
21: Mat C =B;
22: Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data;
23: PetscErrorCode ierr;
24: PetscInt i,j,k,ipvt[15];
25: const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*ajtmp,*bjtmp,*bdiag=b->diag,*pj;
26: PetscInt nz,nzL,row;
27: MatScalar *rtmp,*pc,*mwork,*pv,*vv,work[225];
28: const MatScalar *v,*aa=a->a;
29: PetscInt bs2 = a->bs2,bs=A->rmap->bs,flg;
30: PetscInt sol_ver;
31: PetscBool allowzeropivot,zeropivotdetected;
34: allowzeropivot = PetscNot(A->erroriffailure);
35: PetscOptionsGetInt(NULL,((PetscObject)A)->prefix,"-sol_ver",&sol_ver,NULL);
37: /* generate work space needed by the factorization */
38: PetscMalloc2(bs2*n,&rtmp,bs2,&mwork);
39: PetscMemzero(rtmp,bs2*n*sizeof(MatScalar));
41: for (i=0; i<n; i++) {
42: /* zero rtmp */
43: /* L part */
44: nz = bi[i+1] - bi[i];
45: bjtmp = bj + bi[i];
46: for (j=0; j<nz; j++) {
47: PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));
48: }
50: /* U part */
51: nz = bdiag[i] - bdiag[i+1];
52: bjtmp = bj + bdiag[i+1]+1;
53: for (j=0; j<nz; j++) {
54: PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));
55: }
57: /* load in initial (unfactored row) */
58: nz = ai[i+1] - ai[i];
59: ajtmp = aj + ai[i];
60: v = aa + bs2*ai[i];
61: for (j=0; j<nz; j++) {
62: PetscMemcpy(rtmp+bs2*ajtmp[j],v+bs2*j,bs2*sizeof(MatScalar));
63: }
65: /* elimination */
66: bjtmp = bj + bi[i];
67: nzL = bi[i+1] - bi[i];
68: for (k=0; k < nzL; k++) {
69: row = bjtmp[k];
70: pc = rtmp + bs2*row;
71: for (flg=0,j=0; j<bs2; j++) {
72: if (pc[j]!=0.0) {
73: flg = 1;
74: break;
75: }
76: }
77: if (flg) {
78: pv = b->a + bs2*bdiag[row];
79: PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork);
80: /*PetscKernel_A_gets_A_times_B_15(pc,pv,mwork);*/
81: pj = b->j + bdiag[row+1]+1; /* begining of U(row,:) */
82: pv = b->a + bs2*(bdiag[row+1]+1);
83: nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */
84: for (j=0; j<nz; j++) {
85: vv = rtmp + bs2*pj[j];
86: PetscKernel_A_gets_A_minus_B_times_C(bs,vv,pc,pv);
87: /* PetscKernel_A_gets_A_minus_B_times_C_15(vv,pc,pv); */
88: pv += bs2;
89: }
90: PetscLogFlops(2*bs2*bs*(nz+1)-bs2); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
91: }
92: }
94: /* finished row so stick it into b->a */
95: /* L part */
96: pv = b->a + bs2*bi[i];
97: pj = b->j + bi[i];
98: nz = bi[i+1] - bi[i];
99: for (j=0; j<nz; j++) {
100: PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));
101: }
103: /* Mark diagonal and invert diagonal for simplier triangular solves */
104: pv = b->a + bs2*bdiag[i];
105: pj = b->j + bdiag[i];
106: PetscMemcpy(pv,rtmp+bs2*pj[0],bs2*sizeof(MatScalar));
107: PetscKernel_A_gets_inverse_A_15(pv,ipvt,work,info->shiftamount,allowzeropivot,&zeropivotdetected);
108: if (zeropivotdetected) C->errortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
110: /* U part */
111: pv = b->a + bs2*(bdiag[i+1]+1);
112: pj = b->j + bdiag[i+1]+1;
113: nz = bdiag[i] - bdiag[i+1] - 1;
114: for (j=0; j<nz; j++) {
115: PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));
116: }
117: }
119: PetscFree2(rtmp,mwork);
121: C->ops->solve = MatSolve_SeqBAIJ_15_NaturalOrdering_ver1;
122: C->ops->solvetranspose = MatSolve_SeqBAIJ_N_NaturalOrdering;
123: C->assembled = PETSC_TRUE;
125: PetscLogFlops(1.333333333333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
126: return(0);
127: }
131: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_N(Mat B,Mat A,const MatFactorInfo *info)
132: {
133: Mat C =B;
134: Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data;
135: IS isrow = b->row,isicol = b->icol;
137: const PetscInt *r,*ic;
138: PetscInt i,j,k,n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
139: PetscInt *ajtmp,*bjtmp,nz,nzL,row,*bdiag=b->diag,*pj;
140: MatScalar *rtmp,*pc,*mwork,*v,*pv,*aa=a->a;
141: PetscInt bs=A->rmap->bs,bs2 = a->bs2,*v_pivots,flg;
142: MatScalar *v_work;
143: PetscBool col_identity,row_identity,both_identity;
144: PetscBool allowzeropivot,zeropivotdetected;
147: ISGetIndices(isrow,&r);
148: ISGetIndices(isicol,&ic);
149: allowzeropivot = PetscNot(A->erroriffailure);
151: PetscMalloc1(bs2*n,&rtmp);
152: PetscMemzero(rtmp,bs2*n*sizeof(MatScalar));
154: /* generate work space needed by dense LU factorization */
155: PetscMalloc3(bs,&v_work,bs2,&mwork,bs,&v_pivots);
157: for (i=0; i<n; i++) {
158: /* zero rtmp */
159: /* L part */
160: nz = bi[i+1] - bi[i];
161: bjtmp = bj + bi[i];
162: for (j=0; j<nz; j++) {
163: PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));
164: }
166: /* U part */
167: nz = bdiag[i] - bdiag[i+1];
168: bjtmp = bj + bdiag[i+1]+1;
169: for (j=0; j<nz; j++) {
170: PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));
171: }
173: /* load in initial (unfactored row) */
174: nz = ai[r[i]+1] - ai[r[i]];
175: ajtmp = aj + ai[r[i]];
176: v = aa + bs2*ai[r[i]];
177: for (j=0; j<nz; j++) {
178: PetscMemcpy(rtmp+bs2*ic[ajtmp[j]],v+bs2*j,bs2*sizeof(MatScalar));
179: }
181: /* elimination */
182: bjtmp = bj + bi[i];
183: nzL = bi[i+1] - bi[i];
184: for (k=0; k < nzL; k++) {
185: row = bjtmp[k];
186: pc = rtmp + bs2*row;
187: for (flg=0,j=0; j<bs2; j++) {
188: if (pc[j]!=0.0) {
189: flg = 1;
190: break;
191: }
192: }
193: if (flg) {
194: pv = b->a + bs2*bdiag[row];
195: PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); /* *pc = *pc * (*pv); */
196: pj = b->j + bdiag[row+1]+1; /* begining of U(row,:) */
197: pv = b->a + bs2*(bdiag[row+1]+1);
198: nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */
199: for (j=0; j<nz; j++) {
200: PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j);
201: }
202: PetscLogFlops(2*bs2*bs*(nz+1)-bs2); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
203: }
204: }
206: /* finished row so stick it into b->a */
207: /* L part */
208: pv = b->a + bs2*bi[i];
209: pj = b->j + bi[i];
210: nz = bi[i+1] - bi[i];
211: for (j=0; j<nz; j++) {
212: PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));
213: }
215: /* Mark diagonal and invert diagonal for simplier triangular solves */
216: pv = b->a + bs2*bdiag[i];
217: pj = b->j + bdiag[i];
218: /* if (*pj != i)SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"row %d != *pj %d",i,*pj); */
219: PetscMemcpy(pv,rtmp+bs2*pj[0],bs2*sizeof(MatScalar));
221: PetscKernel_A_gets_inverse_A(bs,pv,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
222: if (zeropivotdetected) B->errortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
224: /* U part */
225: pv = b->a + bs2*(bdiag[i+1]+1);
226: pj = b->j + bdiag[i+1]+1;
227: nz = bdiag[i] - bdiag[i+1] - 1;
228: for (j=0; j<nz; j++) {
229: PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));
230: }
231: }
233: PetscFree(rtmp);
234: PetscFree3(v_work,mwork,v_pivots);
235: ISRestoreIndices(isicol,&ic);
236: ISRestoreIndices(isrow,&r);
238: ISIdentity(isrow,&row_identity);
239: ISIdentity(isicol,&col_identity);
241: both_identity = (PetscBool) (row_identity && col_identity);
242: if (both_identity) {
243: C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
244: } else {
245: C->ops->solve = MatSolve_SeqBAIJ_N;
246: }
247: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N;
249: C->assembled = PETSC_TRUE;
251: PetscLogFlops(1.333333333333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
252: return(0);
253: }
255: /*
256: ilu(0) with natural ordering under new data structure.
257: See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed description
258: because this code is almost identical to MatILUFactorSymbolic_SeqAIJ_ilu0_inplace().
259: */
263: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_ilu0(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
264: {
266: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b;
268: PetscInt n=a->mbs,*ai=a->i,*aj,*adiag=a->diag,bs2 = a->bs2;
269: PetscInt i,j,nz,*bi,*bj,*bdiag,bi_temp;
272: MatDuplicateNoCreate_SeqBAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_FALSE);
273: b = (Mat_SeqBAIJ*)(fact)->data;
275: /* allocate matrix arrays for new data structure */
276: PetscMalloc3(bs2*ai[n]+1,&b->a,ai[n]+1,&b->j,n+1,&b->i);
277: PetscLogObjectMemory((PetscObject)fact,ai[n]*(bs2*sizeof(PetscScalar)+sizeof(PetscInt))+(n+1)*sizeof(PetscInt));
279: b->singlemalloc = PETSC_TRUE;
280: b->free_a = PETSC_TRUE;
281: b->free_ij = PETSC_TRUE;
282: fact->preallocated = PETSC_TRUE;
283: fact->assembled = PETSC_TRUE;
284: if (!b->diag) {
285: PetscMalloc1(n+1,&b->diag);
286: PetscLogObjectMemory((PetscObject)fact,(n+1)*sizeof(PetscInt));
287: }
288: bdiag = b->diag;
290: if (n > 0) {
291: PetscMemzero(b->a,bs2*ai[n]*sizeof(MatScalar));
292: }
294: /* set bi and bj with new data structure */
295: bi = b->i;
296: bj = b->j;
298: /* L part */
299: bi[0] = 0;
300: for (i=0; i<n; i++) {
301: nz = adiag[i] - ai[i];
302: bi[i+1] = bi[i] + nz;
303: aj = a->j + ai[i];
304: for (j=0; j<nz; j++) {
305: *bj = aj[j]; bj++;
306: }
307: }
309: /* U part */
310: bi_temp = bi[n];
311: bdiag[n] = bi[n]-1;
312: for (i=n-1; i>=0; i--) {
313: nz = ai[i+1] - adiag[i] - 1;
314: bi_temp = bi_temp + nz + 1;
315: aj = a->j + adiag[i] + 1;
316: for (j=0; j<nz; j++) {
317: *bj = aj[j]; bj++;
318: }
319: /* diag[i] */
320: *bj = i; bj++;
321: bdiag[i] = bi_temp - 1;
322: }
323: return(0);
324: }
328: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
329: {
330: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b;
331: IS isicol;
332: PetscErrorCode ierr;
333: const PetscInt *r,*ic;
334: PetscInt n=a->mbs,*ai=a->i,*aj=a->j,d;
335: PetscInt *bi,*cols,nnz,*cols_lvl;
336: PetscInt *bdiag,prow,fm,nzbd,reallocs=0,dcount=0;
337: PetscInt i,levels,diagonal_fill;
338: PetscBool col_identity,row_identity,both_identity;
339: PetscReal f;
340: PetscInt nlnk,*lnk,*lnk_lvl=NULL;
341: PetscBT lnkbt;
342: PetscInt nzi,*bj,**bj_ptr,**bjlvl_ptr;
343: PetscFreeSpaceList free_space =NULL,current_space=NULL;
344: PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL;
345: PetscBool missing;
346: PetscInt bs=A->rmap->bs,bs2=a->bs2;
349: if (A->rmap->n != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %D columns %D",A->rmap->n,A->cmap->n);
350: if (bs>1) { /* check shifttype */
351: if (info->shifttype == MAT_SHIFT_NONZERO || info->shifttype == MAT_SHIFT_POSITIVE_DEFINITE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only MAT_SHIFT_NONE and MAT_SHIFT_INBLOCKS are supported for BAIJ matrix");
352: }
354: MatMissingDiagonal(A,&missing,&d);
355: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
357: f = info->fill;
358: levels = (PetscInt)info->levels;
359: diagonal_fill = (PetscInt)info->diagonal_fill;
361: ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
363: ISIdentity(isrow,&row_identity);
364: ISIdentity(iscol,&col_identity);
366: both_identity = (PetscBool) (row_identity && col_identity);
368: if (!levels && both_identity) {
369: /* special case: ilu(0) with natural ordering */
370: MatILUFactorSymbolic_SeqBAIJ_ilu0(fact,A,isrow,iscol,info);
371: MatSeqBAIJSetNumericFactorization(fact,both_identity);
373: fact->factortype = MAT_FACTOR_ILU;
374: (fact)->info.factor_mallocs = 0;
375: (fact)->info.fill_ratio_given = info->fill;
376: (fact)->info.fill_ratio_needed = 1.0;
378: b = (Mat_SeqBAIJ*)(fact)->data;
379: b->row = isrow;
380: b->col = iscol;
381: b->icol = isicol;
382: PetscObjectReference((PetscObject)isrow);
383: PetscObjectReference((PetscObject)iscol);
384: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
386: PetscMalloc1((n+1)*bs,&b->solve_work);
387: return(0);
388: }
390: ISGetIndices(isrow,&r);
391: ISGetIndices(isicol,&ic);
393: /* get new row pointers */
394: PetscMalloc1(n+1,&bi);
395: bi[0] = 0;
396: /* bdiag is location of diagonal in factor */
397: PetscMalloc1(n+1,&bdiag);
398: bdiag[0] = 0;
400: PetscMalloc2(n,&bj_ptr,n,&bjlvl_ptr);
402: /* create a linked list for storing column indices of the active row */
403: nlnk = n + 1;
404: PetscIncompleteLLCreate(n,n,nlnk,lnk,lnk_lvl,lnkbt);
406: /* initial FreeSpace size is f*(ai[n]+1) */
407: PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
408: current_space = free_space;
409: PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space_lvl);
410: current_space_lvl = free_space_lvl;
412: for (i=0; i<n; i++) {
413: nzi = 0;
414: /* copy current row into linked list */
415: nnz = ai[r[i]+1] - ai[r[i]];
416: if (!nnz) 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);
417: cols = aj + ai[r[i]];
418: lnk[i] = -1; /* marker to indicate if diagonal exists */
419: PetscIncompleteLLInit(nnz,cols,n,ic,nlnk,lnk,lnk_lvl,lnkbt);
420: nzi += nlnk;
422: /* make sure diagonal entry is included */
423: if (diagonal_fill && lnk[i] == -1) {
424: fm = n;
425: while (lnk[fm] < i) fm = lnk[fm];
426: lnk[i] = lnk[fm]; /* insert diagonal into linked list */
427: lnk[fm] = i;
428: lnk_lvl[i] = 0;
429: nzi++; dcount++;
430: }
432: /* add pivot rows into the active row */
433: nzbd = 0;
434: prow = lnk[n];
435: while (prow < i) {
436: nnz = bdiag[prow];
437: cols = bj_ptr[prow] + nnz + 1;
438: cols_lvl = bjlvl_ptr[prow] + nnz + 1;
439: nnz = bi[prow+1] - bi[prow] - nnz - 1;
441: PetscILULLAddSorted(nnz,cols,levels,cols_lvl,prow,nlnk,lnk,lnk_lvl,lnkbt,prow);
442: nzi += nlnk;
443: prow = lnk[prow];
444: nzbd++;
445: }
446: bdiag[i] = nzbd;
447: bi[i+1] = bi[i] + nzi;
449: /* if free space is not available, make more free space */
450: if (current_space->local_remaining<nzi) {
451: nnz = PetscIntMultTruncate(2,PetscIntMultTruncate(nzi,(n - i))); /* estimated and max additional space needed */
452: PetscFreeSpaceGet(nnz,¤t_space);
453: PetscFreeSpaceGet(nnz,¤t_space_lvl);
454: reallocs++;
455: }
457: /* copy data into free_space and free_space_lvl, then initialize lnk */
458: PetscIncompleteLLClean(n,n,nzi,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);
460: bj_ptr[i] = current_space->array;
461: bjlvl_ptr[i] = current_space_lvl->array;
463: /* make sure the active row i has diagonal entry */
464: if (*(bj_ptr[i]+bdiag[i]) != i) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Row %D has missing diagonal in factored matrix\ntry running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill",i);
466: current_space->array += nzi;
467: current_space->local_used += nzi;
468: current_space->local_remaining -= nzi;
470: current_space_lvl->array += nzi;
471: current_space_lvl->local_used += nzi;
472: current_space_lvl->local_remaining -= nzi;
473: }
475: ISRestoreIndices(isrow,&r);
476: ISRestoreIndices(isicol,&ic);
478: /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
479: PetscMalloc1(bi[n]+1,&bj);
480: PetscFreeSpaceContiguous_LU(&free_space,bj,n,bi,bdiag);
482: PetscIncompleteLLDestroy(lnk,lnkbt);
483: PetscFreeSpaceDestroy(free_space_lvl);
484: PetscFree2(bj_ptr,bjlvl_ptr);
486: #if defined(PETSC_USE_INFO)
487: {
488: PetscReal af = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]);
489: PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);
490: PetscInfo1(A,"Run with -[sub_]pc_factor_fill %g or use \n",(double)af);
491: PetscInfo1(A,"PCFactorSetFill([sub]pc,%g);\n",(double)af);
492: PetscInfo(A,"for best performance.\n");
493: if (diagonal_fill) {
494: PetscInfo1(A,"Detected and replaced %D missing diagonals\n",dcount);
495: }
496: }
497: #endif
499: /* put together the new matrix */
500: MatSeqBAIJSetPreallocation_SeqBAIJ(fact,bs,MAT_SKIP_ALLOCATION,NULL);
501: PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol);
503: b = (Mat_SeqBAIJ*)(fact)->data;
504: b->free_a = PETSC_TRUE;
505: b->free_ij = PETSC_TRUE;
506: b->singlemalloc = PETSC_FALSE;
508: PetscMalloc1(bs2*(bdiag[0]+1),&b->a);
510: b->j = bj;
511: b->i = bi;
512: b->diag = bdiag;
513: b->free_diag = PETSC_TRUE;
514: b->ilen = 0;
515: b->imax = 0;
516: b->row = isrow;
517: b->col = iscol;
518: PetscObjectReference((PetscObject)isrow);
519: PetscObjectReference((PetscObject)iscol);
520: b->icol = isicol;
522: PetscMalloc1(bs*n+bs,&b->solve_work);
523: /* In b structure: Free imax, ilen, old a, old j.
524: Allocate bdiag, solve_work, new a, new j */
525: PetscLogObjectMemory((PetscObject)fact,(bdiag[0]+1) * (sizeof(PetscInt)+bs2*sizeof(PetscScalar)));
526: b->maxnz = b->nz = bdiag[0]+1;
528: fact->info.factor_mallocs = reallocs;
529: fact->info.fill_ratio_given = f;
530: fact->info.fill_ratio_needed = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]);
532: MatSeqBAIJSetNumericFactorization(fact,both_identity);
533: return(0);
534: }
536: /*
537: This code is virtually identical to MatILUFactorSymbolic_SeqAIJ
538: except that the data structure of Mat_SeqAIJ is slightly different.
539: Not a good example of code reuse.
540: */
543: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_inplace(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
544: {
545: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b;
546: IS isicol;
548: const PetscInt *r,*ic,*ai = a->i,*aj = a->j,*xi;
549: PetscInt prow,n = a->mbs,*ainew,*ajnew,jmax,*fill,nz,*im,*ajfill,*flev,*xitmp;
550: PetscInt *dloc,idx,row,m,fm,nzf,nzi,reallocate = 0,dcount = 0;
551: PetscInt incrlev,nnz,i,bs = A->rmap->bs,bs2 = a->bs2,levels,diagonal_fill,dd;
552: PetscBool col_identity,row_identity,both_identity,flg;
553: PetscReal f;
556: MatMissingDiagonal_SeqBAIJ(A,&flg,&dd);
557: if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix A is missing diagonal entry in row %D",dd);
559: f = info->fill;
560: levels = (PetscInt)info->levels;
561: diagonal_fill = (PetscInt)info->diagonal_fill;
563: ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
565: ISIdentity(isrow,&row_identity);
566: ISIdentity(iscol,&col_identity);
567: both_identity = (PetscBool) (row_identity && col_identity);
569: if (!levels && both_identity) { /* special case copy the nonzero structure */
570: MatDuplicateNoCreate_SeqBAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_TRUE);
571: MatSeqBAIJSetNumericFactorization_inplace(fact,both_identity);
573: fact->factortype = MAT_FACTOR_ILU;
574: b = (Mat_SeqBAIJ*)fact->data;
575: b->row = isrow;
576: b->col = iscol;
577: PetscObjectReference((PetscObject)isrow);
578: PetscObjectReference((PetscObject)iscol);
579: b->icol = isicol;
580: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
582: PetscMalloc1((n+1)*bs,&b->solve_work);
583: return(0);
584: }
586: /* general case perform the symbolic factorization */
587: ISGetIndices(isrow,&r);
588: ISGetIndices(isicol,&ic);
590: /* get new row pointers */
591: PetscMalloc1(n+1,&ainew);
592: ainew[0] = 0;
593: /* don't know how many column pointers are needed so estimate */
594: jmax = (PetscInt)(f*ai[n] + 1);
595: PetscMalloc1(jmax,&ajnew);
596: /* ajfill is level of fill for each fill entry */
597: PetscMalloc1(jmax,&ajfill);
598: /* fill is a linked list of nonzeros in active row */
599: PetscMalloc1(n+1,&fill);
600: /* im is level for each filled value */
601: PetscMalloc1(n+1,&im);
602: /* dloc is location of diagonal in factor */
603: PetscMalloc1(n+1,&dloc);
604: dloc[0] = 0;
605: for (prow=0; prow<n; prow++) {
607: /* copy prow into linked list */
608: nzf = nz = ai[r[prow]+1] - ai[r[prow]];
609: if (!nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",r[prow],prow);
610: xi = aj + ai[r[prow]];
611: fill[n] = n;
612: fill[prow] = -1; /* marker for diagonal entry */
613: while (nz--) {
614: fm = n;
615: idx = ic[*xi++];
616: do {
617: m = fm;
618: fm = fill[m];
619: } while (fm < idx);
620: fill[m] = idx;
621: fill[idx] = fm;
622: im[idx] = 0;
623: }
625: /* make sure diagonal entry is included */
626: if (diagonal_fill && fill[prow] == -1) {
627: fm = n;
628: while (fill[fm] < prow) fm = fill[fm];
629: fill[prow] = fill[fm]; /* insert diagonal into linked list */
630: fill[fm] = prow;
631: im[prow] = 0;
632: nzf++;
633: dcount++;
634: }
636: nzi = 0;
637: row = fill[n];
638: while (row < prow) {
639: incrlev = im[row] + 1;
640: nz = dloc[row];
641: xi = ajnew + ainew[row] + nz + 1;
642: flev = ajfill + ainew[row] + nz + 1;
643: nnz = ainew[row+1] - ainew[row] - nz - 1;
644: fm = row;
645: while (nnz-- > 0) {
646: idx = *xi++;
647: if (*flev + incrlev > levels) {
648: flev++;
649: continue;
650: }
651: do {
652: m = fm;
653: fm = fill[m];
654: } while (fm < idx);
655: if (fm != idx) {
656: im[idx] = *flev + incrlev;
657: fill[m] = idx;
658: fill[idx] = fm;
659: fm = idx;
660: nzf++;
661: } else if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev;
662: flev++;
663: }
664: row = fill[row];
665: nzi++;
666: }
667: /* copy new filled row into permanent storage */
668: ainew[prow+1] = ainew[prow] + nzf;
669: if (ainew[prow+1] > jmax) {
671: /* estimate how much additional space we will need */
672: /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
673: /* just double the memory each time */
674: PetscInt maxadd = jmax;
675: /* maxadd = (int)(((f*ai[n]+1)*(n-prow+5))/n); */
676: if (maxadd < nzf) maxadd = (n-prow)*(nzf+1);
677: jmax += maxadd;
679: /* allocate a longer ajnew and ajfill */
680: PetscMalloc1(jmax,&xitmp);
681: PetscMemcpy(xitmp,ajnew,ainew[prow]*sizeof(PetscInt));
682: PetscFree(ajnew);
683: ajnew = xitmp;
684: PetscMalloc1(jmax,&xitmp);
685: PetscMemcpy(xitmp,ajfill,ainew[prow]*sizeof(PetscInt));
686: PetscFree(ajfill);
687: ajfill = xitmp;
688: reallocate++; /* count how many reallocations are needed */
689: }
690: xitmp = ajnew + ainew[prow];
691: flev = ajfill + ainew[prow];
692: dloc[prow] = nzi;
693: fm = fill[n];
694: while (nzf--) {
695: *xitmp++ = fm;
696: *flev++ = im[fm];
697: fm = fill[fm];
698: }
699: /* make sure row has diagonal entry */
700: if (ajnew[ainew[prow]+dloc[prow]] != prow) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Row %D has missing diagonal in factored matrix\n\
701: try running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill",prow);
702: }
703: PetscFree(ajfill);
704: ISRestoreIndices(isrow,&r);
705: ISRestoreIndices(isicol,&ic);
706: PetscFree(fill);
707: PetscFree(im);
709: #if defined(PETSC_USE_INFO)
710: {
711: PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
712: PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocate,(double)f,(double)af);
713: PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
714: PetscInfo1(A,"PCFactorSetFill(pc,%g);\n",(double)af);
715: PetscInfo(A,"for best performance.\n");
716: if (diagonal_fill) {
717: PetscInfo1(A,"Detected and replaced %D missing diagonals\n",dcount);
718: }
719: }
720: #endif
722: /* put together the new matrix */
723: MatSeqBAIJSetPreallocation_SeqBAIJ(fact,bs,MAT_SKIP_ALLOCATION,NULL);
724: PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol);
725: b = (Mat_SeqBAIJ*)fact->data;
727: b->free_a = PETSC_TRUE;
728: b->free_ij = PETSC_TRUE;
729: b->singlemalloc = PETSC_FALSE;
731: PetscMalloc1(bs2*ainew[n],&b->a);
733: b->j = ajnew;
734: b->i = ainew;
735: for (i=0; i<n; i++) dloc[i] += ainew[i];
736: b->diag = dloc;
737: b->free_diag = PETSC_TRUE;
738: b->ilen = 0;
739: b->imax = 0;
740: b->row = isrow;
741: b->col = iscol;
742: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
744: PetscObjectReference((PetscObject)isrow);
745: PetscObjectReference((PetscObject)iscol);
746: b->icol = isicol;
747: PetscMalloc1(bs*n+bs,&b->solve_work);
748: /* In b structure: Free imax, ilen, old a, old j.
749: Allocate dloc, solve_work, new a, new j */
750: PetscLogObjectMemory((PetscObject)fact,(ainew[n]-n)*(sizeof(PetscInt))+bs2*ainew[n]*sizeof(PetscScalar));
751: b->maxnz = b->nz = ainew[n];
753: fact->info.factor_mallocs = reallocate;
754: fact->info.fill_ratio_given = f;
755: fact->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[prow]);
757: MatSeqBAIJSetNumericFactorization_inplace(fact,both_identity);
758: return(0);
759: }
763: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE(Mat A)
764: {
765: /* Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; */
766: /* int i,*AJ=a->j,nz=a->nz; */
769: /* Undo Column scaling */
770: /* while (nz--) { */
771: /* AJ[i] = AJ[i]/4; */
772: /* } */
773: /* This should really invoke a push/pop logic, but we don't have that yet. */
774: A->ops->setunfactored = NULL;
775: return(0);
776: }
780: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE_usj(Mat A)
781: {
782: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
783: PetscInt *AJ=a->j,nz=a->nz;
784: unsigned short *aj=(unsigned short*)AJ;
787: /* Is this really necessary? */
788: while (nz--) {
789: AJ[nz] = (int)((unsigned int)aj[nz]); /* First extend, then convert to signed. */
790: }
791: A->ops->setunfactored = NULL;
792: return(0);
793: }