Actual source code: baijfact2.c
petsc-3.8.4 2018-03-24
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);
14: /*
15: This is not much faster than MatLUFactorNumeric_SeqBAIJ_N() but the solve is faster at least sometimes
16: */
17: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_15_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info)
18: {
19: Mat C =B;
20: Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data;
21: PetscErrorCode ierr;
22: PetscInt i,j,k,ipvt[15];
23: const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*ajtmp,*bjtmp,*bdiag=b->diag,*pj;
24: PetscInt nz,nzL,row;
25: MatScalar *rtmp,*pc,*mwork,*pv,*vv,work[225];
26: const MatScalar *v,*aa=a->a;
27: PetscInt bs2 = a->bs2,bs=A->rmap->bs,flg;
28: PetscInt sol_ver;
29: PetscBool allowzeropivot,zeropivotdetected;
32: allowzeropivot = PetscNot(A->erroriffailure);
33: PetscOptionsGetInt(NULL,((PetscObject)A)->prefix,"-sol_ver",&sol_ver,NULL);
35: /* generate work space needed by the factorization */
36: PetscMalloc2(bs2*n,&rtmp,bs2,&mwork);
37: PetscMemzero(rtmp,bs2*n*sizeof(MatScalar));
39: for (i=0; i<n; i++) {
40: /* zero rtmp */
41: /* L part */
42: nz = bi[i+1] - bi[i];
43: bjtmp = bj + bi[i];
44: for (j=0; j<nz; j++) {
45: PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));
46: }
48: /* U part */
49: nz = bdiag[i] - bdiag[i+1];
50: bjtmp = bj + bdiag[i+1]+1;
51: for (j=0; j<nz; j++) {
52: PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));
53: }
55: /* load in initial (unfactored row) */
56: nz = ai[i+1] - ai[i];
57: ajtmp = aj + ai[i];
58: v = aa + bs2*ai[i];
59: for (j=0; j<nz; j++) {
60: PetscMemcpy(rtmp+bs2*ajtmp[j],v+bs2*j,bs2*sizeof(MatScalar));
61: }
63: /* elimination */
64: bjtmp = bj + bi[i];
65: nzL = bi[i+1] - bi[i];
66: for (k=0; k < nzL; k++) {
67: row = bjtmp[k];
68: pc = rtmp + bs2*row;
69: for (flg=0,j=0; j<bs2; j++) {
70: if (pc[j]!=0.0) {
71: flg = 1;
72: break;
73: }
74: }
75: if (flg) {
76: pv = b->a + bs2*bdiag[row];
77: PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork);
78: /*PetscKernel_A_gets_A_times_B_15(pc,pv,mwork);*/
79: pj = b->j + bdiag[row+1]+1; /* begining of U(row,:) */
80: pv = b->a + bs2*(bdiag[row+1]+1);
81: nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */
82: for (j=0; j<nz; j++) {
83: vv = rtmp + bs2*pj[j];
84: PetscKernel_A_gets_A_minus_B_times_C(bs,vv,pc,pv);
85: /* PetscKernel_A_gets_A_minus_B_times_C_15(vv,pc,pv); */
86: pv += bs2;
87: }
88: PetscLogFlops(2*bs2*bs*(nz+1)-bs2); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
89: }
90: }
92: /* finished row so stick it into b->a */
93: /* L part */
94: pv = b->a + bs2*bi[i];
95: pj = b->j + bi[i];
96: nz = bi[i+1] - bi[i];
97: for (j=0; j<nz; j++) {
98: PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));
99: }
101: /* Mark diagonal and invert diagonal for simplier triangular solves */
102: pv = b->a + bs2*bdiag[i];
103: pj = b->j + bdiag[i];
104: PetscMemcpy(pv,rtmp+bs2*pj[0],bs2*sizeof(MatScalar));
105: PetscKernel_A_gets_inverse_A_15(pv,ipvt,work,info->shiftamount,allowzeropivot,&zeropivotdetected);
106: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
108: /* U part */
109: pv = b->a + bs2*(bdiag[i+1]+1);
110: pj = b->j + bdiag[i+1]+1;
111: nz = bdiag[i] - bdiag[i+1] - 1;
112: for (j=0; j<nz; j++) {
113: PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));
114: }
115: }
117: PetscFree2(rtmp,mwork);
119: C->ops->solve = MatSolve_SeqBAIJ_15_NaturalOrdering_ver1;
120: C->ops->solvetranspose = MatSolve_SeqBAIJ_N_NaturalOrdering;
121: C->assembled = PETSC_TRUE;
123: PetscLogFlops(1.333333333333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
124: return(0);
125: }
127: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_N(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;
131: IS isrow = b->row,isicol = b->icol;
133: const PetscInt *r,*ic;
134: PetscInt i,j,k,n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
135: PetscInt *ajtmp,*bjtmp,nz,nzL,row,*bdiag=b->diag,*pj;
136: MatScalar *rtmp,*pc,*mwork,*v,*pv,*aa=a->a;
137: PetscInt bs=A->rmap->bs,bs2 = a->bs2,*v_pivots,flg;
138: MatScalar *v_work;
139: PetscBool col_identity,row_identity,both_identity;
140: PetscBool allowzeropivot,zeropivotdetected;
143: ISGetIndices(isrow,&r);
144: ISGetIndices(isicol,&ic);
145: allowzeropivot = PetscNot(A->erroriffailure);
147: PetscMalloc1(bs2*n,&rtmp);
148: PetscMemzero(rtmp,bs2*n*sizeof(MatScalar));
150: /* generate work space needed by dense LU factorization */
151: PetscMalloc3(bs,&v_work,bs2,&mwork,bs,&v_pivots);
153: for (i=0; i<n; i++) {
154: /* zero rtmp */
155: /* L part */
156: nz = bi[i+1] - bi[i];
157: bjtmp = bj + bi[i];
158: for (j=0; j<nz; j++) {
159: PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));
160: }
162: /* U part */
163: nz = bdiag[i] - bdiag[i+1];
164: bjtmp = bj + bdiag[i+1]+1;
165: for (j=0; j<nz; j++) {
166: PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));
167: }
169: /* load in initial (unfactored row) */
170: nz = ai[r[i]+1] - ai[r[i]];
171: ajtmp = aj + ai[r[i]];
172: v = aa + bs2*ai[r[i]];
173: for (j=0; j<nz; j++) {
174: PetscMemcpy(rtmp+bs2*ic[ajtmp[j]],v+bs2*j,bs2*sizeof(MatScalar));
175: }
177: /* elimination */
178: bjtmp = bj + bi[i];
179: nzL = bi[i+1] - bi[i];
180: for (k=0; k < nzL; k++) {
181: row = bjtmp[k];
182: pc = rtmp + bs2*row;
183: for (flg=0,j=0; j<bs2; j++) {
184: if (pc[j]!=0.0) {
185: flg = 1;
186: break;
187: }
188: }
189: if (flg) {
190: pv = b->a + bs2*bdiag[row];
191: PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); /* *pc = *pc * (*pv); */
192: pj = b->j + bdiag[row+1]+1; /* begining of U(row,:) */
193: pv = b->a + bs2*(bdiag[row+1]+1);
194: nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */
195: for (j=0; j<nz; j++) {
196: PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j);
197: }
198: PetscLogFlops(2*bs2*bs*(nz+1)-bs2); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
199: }
200: }
202: /* finished row so stick it into b->a */
203: /* L part */
204: pv = b->a + bs2*bi[i];
205: pj = b->j + bi[i];
206: nz = bi[i+1] - bi[i];
207: for (j=0; j<nz; j++) {
208: PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));
209: }
211: /* Mark diagonal and invert diagonal for simplier triangular solves */
212: pv = b->a + bs2*bdiag[i];
213: pj = b->j + bdiag[i];
214: PetscMemcpy(pv,rtmp+bs2*pj[0],bs2*sizeof(MatScalar));
216: PetscKernel_A_gets_inverse_A(bs,pv,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
217: if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
219: /* U part */
220: pv = b->a + bs2*(bdiag[i+1]+1);
221: pj = b->j + bdiag[i+1]+1;
222: nz = bdiag[i] - bdiag[i+1] - 1;
223: for (j=0; j<nz; j++) {
224: PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));
225: }
226: }
228: PetscFree(rtmp);
229: PetscFree3(v_work,mwork,v_pivots);
230: ISRestoreIndices(isicol,&ic);
231: ISRestoreIndices(isrow,&r);
233: ISIdentity(isrow,&row_identity);
234: ISIdentity(isicol,&col_identity);
236: both_identity = (PetscBool) (row_identity && col_identity);
237: if (both_identity) {
238: switch (bs) {
239: case 11:
240: C->ops->solve = MatSolve_SeqBAIJ_11_NaturalOrdering;
241: break;
242: case 12:
243: C->ops->solve = MatSolve_SeqBAIJ_12_NaturalOrdering;
244: break;
245: case 13:
246: C->ops->solve = MatSolve_SeqBAIJ_13_NaturalOrdering;
247: break;
248: case 14:
249: C->ops->solve = MatSolve_SeqBAIJ_14_NaturalOrdering;
250: break;
251: default:
252: C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
253: break;
254: }
255: } else {
256: C->ops->solve = MatSolve_SeqBAIJ_N;
257: }
258: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N;
260: C->assembled = PETSC_TRUE;
262: PetscLogFlops(1.333333333333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
263: return(0);
264: }
266: /*
267: ilu(0) with natural ordering under new data structure.
268: See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed description
269: because this code is almost identical to MatILUFactorSymbolic_SeqAIJ_ilu0_inplace().
270: */
272: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_ilu0(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
273: {
275: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b;
277: PetscInt n=a->mbs,*ai=a->i,*aj,*adiag=a->diag,bs2 = a->bs2;
278: PetscInt i,j,nz,*bi,*bj,*bdiag,bi_temp;
281: MatDuplicateNoCreate_SeqBAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_FALSE);
282: b = (Mat_SeqBAIJ*)(fact)->data;
284: /* allocate matrix arrays for new data structure */
285: PetscMalloc3(bs2*ai[n]+1,&b->a,ai[n]+1,&b->j,n+1,&b->i);
286: PetscLogObjectMemory((PetscObject)fact,ai[n]*(bs2*sizeof(PetscScalar)+sizeof(PetscInt))+(n+1)*sizeof(PetscInt));
288: b->singlemalloc = PETSC_TRUE;
289: b->free_a = PETSC_TRUE;
290: b->free_ij = PETSC_TRUE;
291: fact->preallocated = PETSC_TRUE;
292: fact->assembled = PETSC_TRUE;
293: if (!b->diag) {
294: PetscMalloc1(n+1,&b->diag);
295: PetscLogObjectMemory((PetscObject)fact,(n+1)*sizeof(PetscInt));
296: }
297: bdiag = b->diag;
299: if (n > 0) {
300: PetscMemzero(b->a,bs2*ai[n]*sizeof(MatScalar));
301: }
303: /* set bi and bj with new data structure */
304: bi = b->i;
305: bj = b->j;
307: /* L part */
308: bi[0] = 0;
309: for (i=0; i<n; i++) {
310: nz = adiag[i] - ai[i];
311: bi[i+1] = bi[i] + nz;
312: aj = a->j + ai[i];
313: for (j=0; j<nz; j++) {
314: *bj = aj[j]; bj++;
315: }
316: }
318: /* U part */
319: bi_temp = bi[n];
320: bdiag[n] = bi[n]-1;
321: for (i=n-1; i>=0; i--) {
322: nz = ai[i+1] - adiag[i] - 1;
323: bi_temp = bi_temp + nz + 1;
324: aj = a->j + adiag[i] + 1;
325: for (j=0; j<nz; j++) {
326: *bj = aj[j]; bj++;
327: }
328: /* diag[i] */
329: *bj = i; bj++;
330: bdiag[i] = bi_temp - 1;
331: }
332: return(0);
333: }
335: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
336: {
337: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b;
338: IS isicol;
339: PetscErrorCode ierr;
340: const PetscInt *r,*ic;
341: PetscInt n=a->mbs,*ai=a->i,*aj=a->j,d;
342: PetscInt *bi,*cols,nnz,*cols_lvl;
343: PetscInt *bdiag,prow,fm,nzbd,reallocs=0,dcount=0;
344: PetscInt i,levels,diagonal_fill;
345: PetscBool col_identity,row_identity,both_identity;
346: PetscReal f;
347: PetscInt nlnk,*lnk,*lnk_lvl=NULL;
348: PetscBT lnkbt;
349: PetscInt nzi,*bj,**bj_ptr,**bjlvl_ptr;
350: PetscFreeSpaceList free_space =NULL,current_space=NULL;
351: PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL;
352: PetscBool missing;
353: PetscInt bs=A->rmap->bs,bs2=a->bs2;
356: 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);
357: if (bs>1) { /* check shifttype */
358: 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");
359: }
361: MatMissingDiagonal(A,&missing,&d);
362: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
364: f = info->fill;
365: levels = (PetscInt)info->levels;
366: diagonal_fill = (PetscInt)info->diagonal_fill;
368: ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
370: ISIdentity(isrow,&row_identity);
371: ISIdentity(iscol,&col_identity);
373: both_identity = (PetscBool) (row_identity && col_identity);
375: if (!levels && both_identity) {
376: /* special case: ilu(0) with natural ordering */
377: MatILUFactorSymbolic_SeqBAIJ_ilu0(fact,A,isrow,iscol,info);
378: MatSeqBAIJSetNumericFactorization(fact,both_identity);
380: fact->factortype = MAT_FACTOR_ILU;
381: (fact)->info.factor_mallocs = 0;
382: (fact)->info.fill_ratio_given = info->fill;
383: (fact)->info.fill_ratio_needed = 1.0;
385: b = (Mat_SeqBAIJ*)(fact)->data;
386: b->row = isrow;
387: b->col = iscol;
388: b->icol = isicol;
389: PetscObjectReference((PetscObject)isrow);
390: PetscObjectReference((PetscObject)iscol);
391: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
393: PetscMalloc1((n+1)*bs,&b->solve_work);
394: return(0);
395: }
397: ISGetIndices(isrow,&r);
398: ISGetIndices(isicol,&ic);
400: /* get new row pointers */
401: PetscMalloc1(n+1,&bi);
402: bi[0] = 0;
403: /* bdiag is location of diagonal in factor */
404: PetscMalloc1(n+1,&bdiag);
405: bdiag[0] = 0;
407: PetscMalloc2(n,&bj_ptr,n,&bjlvl_ptr);
409: /* create a linked list for storing column indices of the active row */
410: nlnk = n + 1;
411: PetscIncompleteLLCreate(n,n,nlnk,lnk,lnk_lvl,lnkbt);
413: /* initial FreeSpace size is f*(ai[n]+1) */
414: PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
415: current_space = free_space;
416: PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space_lvl);
417: current_space_lvl = free_space_lvl;
419: for (i=0; i<n; i++) {
420: nzi = 0;
421: /* copy current row into linked list */
422: nnz = ai[r[i]+1] - ai[r[i]];
423: 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);
424: cols = aj + ai[r[i]];
425: lnk[i] = -1; /* marker to indicate if diagonal exists */
426: PetscIncompleteLLInit(nnz,cols,n,ic,nlnk,lnk,lnk_lvl,lnkbt);
427: nzi += nlnk;
429: /* make sure diagonal entry is included */
430: if (diagonal_fill && lnk[i] == -1) {
431: fm = n;
432: while (lnk[fm] < i) fm = lnk[fm];
433: lnk[i] = lnk[fm]; /* insert diagonal into linked list */
434: lnk[fm] = i;
435: lnk_lvl[i] = 0;
436: nzi++; dcount++;
437: }
439: /* add pivot rows into the active row */
440: nzbd = 0;
441: prow = lnk[n];
442: while (prow < i) {
443: nnz = bdiag[prow];
444: cols = bj_ptr[prow] + nnz + 1;
445: cols_lvl = bjlvl_ptr[prow] + nnz + 1;
446: nnz = bi[prow+1] - bi[prow] - nnz - 1;
448: PetscILULLAddSorted(nnz,cols,levels,cols_lvl,prow,nlnk,lnk,lnk_lvl,lnkbt,prow);
449: nzi += nlnk;
450: prow = lnk[prow];
451: nzbd++;
452: }
453: bdiag[i] = nzbd;
454: bi[i+1] = bi[i] + nzi;
456: /* if free space is not available, make more free space */
457: if (current_space->local_remaining<nzi) {
458: nnz = PetscIntMultTruncate(2,PetscIntMultTruncate(nzi,(n - i))); /* estimated and max additional space needed */
459: PetscFreeSpaceGet(nnz,¤t_space);
460: PetscFreeSpaceGet(nnz,¤t_space_lvl);
461: reallocs++;
462: }
464: /* copy data into free_space and free_space_lvl, then initialize lnk */
465: PetscIncompleteLLClean(n,n,nzi,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);
467: bj_ptr[i] = current_space->array;
468: bjlvl_ptr[i] = current_space_lvl->array;
470: /* make sure the active row i has diagonal entry */
471: 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);
473: current_space->array += nzi;
474: current_space->local_used += nzi;
475: current_space->local_remaining -= nzi;
477: current_space_lvl->array += nzi;
478: current_space_lvl->local_used += nzi;
479: current_space_lvl->local_remaining -= nzi;
480: }
482: ISRestoreIndices(isrow,&r);
483: ISRestoreIndices(isicol,&ic);
485: /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
486: PetscMalloc1(bi[n]+1,&bj);
487: PetscFreeSpaceContiguous_LU(&free_space,bj,n,bi,bdiag);
489: PetscIncompleteLLDestroy(lnk,lnkbt);
490: PetscFreeSpaceDestroy(free_space_lvl);
491: PetscFree2(bj_ptr,bjlvl_ptr);
493: #if defined(PETSC_USE_INFO)
494: {
495: PetscReal af = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]);
496: PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);
497: PetscInfo1(A,"Run with -[sub_]pc_factor_fill %g or use \n",(double)af);
498: PetscInfo1(A,"PCFactorSetFill([sub]pc,%g);\n",(double)af);
499: PetscInfo(A,"for best performance.\n");
500: if (diagonal_fill) {
501: PetscInfo1(A,"Detected and replaced %D missing diagonals\n",dcount);
502: }
503: }
504: #endif
506: /* put together the new matrix */
507: MatSeqBAIJSetPreallocation(fact,bs,MAT_SKIP_ALLOCATION,NULL);
508: PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol);
510: b = (Mat_SeqBAIJ*)(fact)->data;
511: b->free_a = PETSC_TRUE;
512: b->free_ij = PETSC_TRUE;
513: b->singlemalloc = PETSC_FALSE;
515: PetscMalloc1(bs2*(bdiag[0]+1),&b->a);
517: b->j = bj;
518: b->i = bi;
519: b->diag = bdiag;
520: b->free_diag = PETSC_TRUE;
521: b->ilen = 0;
522: b->imax = 0;
523: b->row = isrow;
524: b->col = iscol;
525: PetscObjectReference((PetscObject)isrow);
526: PetscObjectReference((PetscObject)iscol);
527: b->icol = isicol;
529: PetscMalloc1(bs*n+bs,&b->solve_work);
530: /* In b structure: Free imax, ilen, old a, old j.
531: Allocate bdiag, solve_work, new a, new j */
532: PetscLogObjectMemory((PetscObject)fact,(bdiag[0]+1) * (sizeof(PetscInt)+bs2*sizeof(PetscScalar)));
533: b->maxnz = b->nz = bdiag[0]+1;
535: fact->info.factor_mallocs = reallocs;
536: fact->info.fill_ratio_given = f;
537: fact->info.fill_ratio_needed = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]);
539: MatSeqBAIJSetNumericFactorization(fact,both_identity);
540: return(0);
541: }
543: /*
544: This code is virtually identical to MatILUFactorSymbolic_SeqAIJ
545: except that the data structure of Mat_SeqAIJ is slightly different.
546: Not a good example of code reuse.
547: */
548: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_inplace(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
549: {
550: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b;
551: IS isicol;
553: const PetscInt *r,*ic,*ai = a->i,*aj = a->j,*xi;
554: PetscInt prow,n = a->mbs,*ainew,*ajnew,jmax,*fill,nz,*im,*ajfill,*flev,*xitmp;
555: PetscInt *dloc,idx,row,m,fm,nzf,nzi,reallocate = 0,dcount = 0;
556: PetscInt incrlev,nnz,i,bs = A->rmap->bs,bs2 = a->bs2,levels,diagonal_fill,dd;
557: PetscBool col_identity,row_identity,both_identity,flg;
558: PetscReal f;
561: MatMissingDiagonal_SeqBAIJ(A,&flg,&dd);
562: if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix A is missing diagonal entry in row %D",dd);
564: f = info->fill;
565: levels = (PetscInt)info->levels;
566: diagonal_fill = (PetscInt)info->diagonal_fill;
568: ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
570: ISIdentity(isrow,&row_identity);
571: ISIdentity(iscol,&col_identity);
572: both_identity = (PetscBool) (row_identity && col_identity);
574: if (!levels && both_identity) { /* special case copy the nonzero structure */
575: MatDuplicateNoCreate_SeqBAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_TRUE);
576: MatSeqBAIJSetNumericFactorization_inplace(fact,both_identity);
578: fact->factortype = MAT_FACTOR_ILU;
579: b = (Mat_SeqBAIJ*)fact->data;
580: b->row = isrow;
581: b->col = iscol;
582: PetscObjectReference((PetscObject)isrow);
583: PetscObjectReference((PetscObject)iscol);
584: b->icol = isicol;
585: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
587: PetscMalloc1((n+1)*bs,&b->solve_work);
588: return(0);
589: }
591: /* general case perform the symbolic factorization */
592: ISGetIndices(isrow,&r);
593: ISGetIndices(isicol,&ic);
595: /* get new row pointers */
596: PetscMalloc1(n+1,&ainew);
597: ainew[0] = 0;
598: /* don't know how many column pointers are needed so estimate */
599: jmax = (PetscInt)(f*ai[n] + 1);
600: PetscMalloc1(jmax,&ajnew);
601: /* ajfill is level of fill for each fill entry */
602: PetscMalloc1(jmax,&ajfill);
603: /* fill is a linked list of nonzeros in active row */
604: PetscMalloc1(n+1,&fill);
605: /* im is level for each filled value */
606: PetscMalloc1(n+1,&im);
607: /* dloc is location of diagonal in factor */
608: PetscMalloc1(n+1,&dloc);
609: dloc[0] = 0;
610: for (prow=0; prow<n; prow++) {
612: /* copy prow into linked list */
613: nzf = nz = ai[r[prow]+1] - ai[r[prow]];
614: 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);
615: xi = aj + ai[r[prow]];
616: fill[n] = n;
617: fill[prow] = -1; /* marker for diagonal entry */
618: while (nz--) {
619: fm = n;
620: idx = ic[*xi++];
621: do {
622: m = fm;
623: fm = fill[m];
624: } while (fm < idx);
625: fill[m] = idx;
626: fill[idx] = fm;
627: im[idx] = 0;
628: }
630: /* make sure diagonal entry is included */
631: if (diagonal_fill && fill[prow] == -1) {
632: fm = n;
633: while (fill[fm] < prow) fm = fill[fm];
634: fill[prow] = fill[fm]; /* insert diagonal into linked list */
635: fill[fm] = prow;
636: im[prow] = 0;
637: nzf++;
638: dcount++;
639: }
641: nzi = 0;
642: row = fill[n];
643: while (row < prow) {
644: incrlev = im[row] + 1;
645: nz = dloc[row];
646: xi = ajnew + ainew[row] + nz + 1;
647: flev = ajfill + ainew[row] + nz + 1;
648: nnz = ainew[row+1] - ainew[row] - nz - 1;
649: fm = row;
650: while (nnz-- > 0) {
651: idx = *xi++;
652: if (*flev + incrlev > levels) {
653: flev++;
654: continue;
655: }
656: do {
657: m = fm;
658: fm = fill[m];
659: } while (fm < idx);
660: if (fm != idx) {
661: im[idx] = *flev + incrlev;
662: fill[m] = idx;
663: fill[idx] = fm;
664: fm = idx;
665: nzf++;
666: } else if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev;
667: flev++;
668: }
669: row = fill[row];
670: nzi++;
671: }
672: /* copy new filled row into permanent storage */
673: ainew[prow+1] = ainew[prow] + nzf;
674: if (ainew[prow+1] > jmax) {
676: /* estimate how much additional space we will need */
677: /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
678: /* just double the memory each time */
679: PetscInt maxadd = jmax;
680: /* maxadd = (int)(((f*ai[n]+1)*(n-prow+5))/n); */
681: if (maxadd < nzf) maxadd = (n-prow)*(nzf+1);
682: jmax += maxadd;
684: /* allocate a longer ajnew and ajfill */
685: PetscMalloc1(jmax,&xitmp);
686: PetscMemcpy(xitmp,ajnew,ainew[prow]*sizeof(PetscInt));
687: PetscFree(ajnew);
688: ajnew = xitmp;
689: PetscMalloc1(jmax,&xitmp);
690: PetscMemcpy(xitmp,ajfill,ainew[prow]*sizeof(PetscInt));
691: PetscFree(ajfill);
692: ajfill = xitmp;
693: reallocate++; /* count how many reallocations are needed */
694: }
695: xitmp = ajnew + ainew[prow];
696: flev = ajfill + ainew[prow];
697: dloc[prow] = nzi;
698: fm = fill[n];
699: while (nzf--) {
700: *xitmp++ = fm;
701: *flev++ = im[fm];
702: fm = fill[fm];
703: }
704: /* make sure row has diagonal entry */
705: 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\
706: try running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill",prow);
707: }
708: PetscFree(ajfill);
709: ISRestoreIndices(isrow,&r);
710: ISRestoreIndices(isicol,&ic);
711: PetscFree(fill);
712: PetscFree(im);
714: #if defined(PETSC_USE_INFO)
715: {
716: PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
717: PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocate,(double)f,(double)af);
718: PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
719: PetscInfo1(A,"PCFactorSetFill(pc,%g);\n",(double)af);
720: PetscInfo(A,"for best performance.\n");
721: if (diagonal_fill) {
722: PetscInfo1(A,"Detected and replaced %D missing diagonals\n",dcount);
723: }
724: }
725: #endif
727: /* put together the new matrix */
728: MatSeqBAIJSetPreallocation(fact,bs,MAT_SKIP_ALLOCATION,NULL);
729: PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol);
730: b = (Mat_SeqBAIJ*)fact->data;
732: b->free_a = PETSC_TRUE;
733: b->free_ij = PETSC_TRUE;
734: b->singlemalloc = PETSC_FALSE;
736: PetscMalloc1(bs2*ainew[n],&b->a);
738: b->j = ajnew;
739: b->i = ainew;
740: for (i=0; i<n; i++) dloc[i] += ainew[i];
741: b->diag = dloc;
742: b->free_diag = PETSC_TRUE;
743: b->ilen = 0;
744: b->imax = 0;
745: b->row = isrow;
746: b->col = iscol;
747: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
749: PetscObjectReference((PetscObject)isrow);
750: PetscObjectReference((PetscObject)iscol);
751: b->icol = isicol;
752: PetscMalloc1(bs*n+bs,&b->solve_work);
753: /* In b structure: Free imax, ilen, old a, old j.
754: Allocate dloc, solve_work, new a, new j */
755: PetscLogObjectMemory((PetscObject)fact,(ainew[n]-n)*(sizeof(PetscInt))+bs2*ainew[n]*sizeof(PetscScalar));
756: b->maxnz = b->nz = ainew[n];
758: fact->info.factor_mallocs = reallocate;
759: fact->info.fill_ratio_given = f;
760: fact->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[prow]);
762: MatSeqBAIJSetNumericFactorization_inplace(fact,both_identity);
763: return(0);
764: }
766: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE(Mat A)
767: {
768: /* Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; */
769: /* int i,*AJ=a->j,nz=a->nz; */
772: /* Undo Column scaling */
773: /* while (nz--) { */
774: /* AJ[i] = AJ[i]/4; */
775: /* } */
776: /* This should really invoke a push/pop logic, but we don't have that yet. */
777: A->ops->setunfactored = NULL;
778: return(0);
779: }
781: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE_usj(Mat A)
782: {
783: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
784: PetscInt *AJ=a->j,nz=a->nz;
785: unsigned short *aj=(unsigned short*)AJ;
788: /* Is this really necessary? */
789: while (nz--) {
790: AJ[nz] = (int)((unsigned int)aj[nz]); /* First extend, then convert to signed. */
791: }
792: A->ops->setunfactored = NULL;
793: return(0);
794: }