Actual source code: baijfact2.c
petsc-3.9.4 2018-09-11
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 9:
240: #if defined(PETSC_HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && defined(PETSC_USE_REAL_DOUBLE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_64BIT_INDICES)
241: C->ops->solve = MatSolve_SeqBAIJ_9_NaturalOrdering;
242: #else
243: C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
244: #endif
245: break;
246: case 11:
247: C->ops->solve = MatSolve_SeqBAIJ_11_NaturalOrdering;
248: break;
249: case 12:
250: C->ops->solve = MatSolve_SeqBAIJ_12_NaturalOrdering;
251: break;
252: case 13:
253: C->ops->solve = MatSolve_SeqBAIJ_13_NaturalOrdering;
254: break;
255: case 14:
256: C->ops->solve = MatSolve_SeqBAIJ_14_NaturalOrdering;
257: break;
258: default:
259: C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
260: break;
261: }
262: } else {
263: C->ops->solve = MatSolve_SeqBAIJ_N;
264: }
265: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N;
267: C->assembled = PETSC_TRUE;
269: PetscLogFlops(1.333333333333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
270: return(0);
271: }
273: /*
274: ilu(0) with natural ordering under new data structure.
275: See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed description
276: because this code is almost identical to MatILUFactorSymbolic_SeqAIJ_ilu0_inplace().
277: */
279: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_ilu0(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
280: {
282: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b;
284: PetscInt n=a->mbs,*ai=a->i,*aj,*adiag=a->diag,bs2 = a->bs2;
285: PetscInt i,j,nz,*bi,*bj,*bdiag,bi_temp;
288: MatDuplicateNoCreate_SeqBAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_FALSE);
289: b = (Mat_SeqBAIJ*)(fact)->data;
291: /* allocate matrix arrays for new data structure */
292: PetscMalloc3(bs2*ai[n]+1,&b->a,ai[n]+1,&b->j,n+1,&b->i);
293: PetscLogObjectMemory((PetscObject)fact,ai[n]*(bs2*sizeof(PetscScalar)+sizeof(PetscInt))+(n+1)*sizeof(PetscInt));
295: b->singlemalloc = PETSC_TRUE;
296: b->free_a = PETSC_TRUE;
297: b->free_ij = PETSC_TRUE;
298: fact->preallocated = PETSC_TRUE;
299: fact->assembled = PETSC_TRUE;
300: if (!b->diag) {
301: PetscMalloc1(n+1,&b->diag);
302: PetscLogObjectMemory((PetscObject)fact,(n+1)*sizeof(PetscInt));
303: }
304: bdiag = b->diag;
306: if (n > 0) {
307: PetscMemzero(b->a,bs2*ai[n]*sizeof(MatScalar));
308: }
310: /* set bi and bj with new data structure */
311: bi = b->i;
312: bj = b->j;
314: /* L part */
315: bi[0] = 0;
316: for (i=0; i<n; i++) {
317: nz = adiag[i] - ai[i];
318: bi[i+1] = bi[i] + nz;
319: aj = a->j + ai[i];
320: for (j=0; j<nz; j++) {
321: *bj = aj[j]; bj++;
322: }
323: }
325: /* U part */
326: bi_temp = bi[n];
327: bdiag[n] = bi[n]-1;
328: for (i=n-1; i>=0; i--) {
329: nz = ai[i+1] - adiag[i] - 1;
330: bi_temp = bi_temp + nz + 1;
331: aj = a->j + adiag[i] + 1;
332: for (j=0; j<nz; j++) {
333: *bj = aj[j]; bj++;
334: }
335: /* diag[i] */
336: *bj = i; bj++;
337: bdiag[i] = bi_temp - 1;
338: }
339: return(0);
340: }
342: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
343: {
344: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b;
345: IS isicol;
346: PetscErrorCode ierr;
347: const PetscInt *r,*ic;
348: PetscInt n=a->mbs,*ai=a->i,*aj=a->j,d;
349: PetscInt *bi,*cols,nnz,*cols_lvl;
350: PetscInt *bdiag,prow,fm,nzbd,reallocs=0,dcount=0;
351: PetscInt i,levels,diagonal_fill;
352: PetscBool col_identity,row_identity,both_identity;
353: PetscReal f;
354: PetscInt nlnk,*lnk,*lnk_lvl=NULL;
355: PetscBT lnkbt;
356: PetscInt nzi,*bj,**bj_ptr,**bjlvl_ptr;
357: PetscFreeSpaceList free_space =NULL,current_space=NULL;
358: PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL;
359: PetscBool missing;
360: PetscInt bs=A->rmap->bs,bs2=a->bs2;
363: 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);
364: if (bs>1) { /* check shifttype */
365: 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");
366: }
368: MatMissingDiagonal(A,&missing,&d);
369: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d);
371: f = info->fill;
372: levels = (PetscInt)info->levels;
373: diagonal_fill = (PetscInt)info->diagonal_fill;
375: ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
377: ISIdentity(isrow,&row_identity);
378: ISIdentity(iscol,&col_identity);
380: both_identity = (PetscBool) (row_identity && col_identity);
382: if (!levels && both_identity) {
383: /* special case: ilu(0) with natural ordering */
384: MatILUFactorSymbolic_SeqBAIJ_ilu0(fact,A,isrow,iscol,info);
385: MatSeqBAIJSetNumericFactorization(fact,both_identity);
387: fact->factortype = MAT_FACTOR_ILU;
388: (fact)->info.factor_mallocs = 0;
389: (fact)->info.fill_ratio_given = info->fill;
390: (fact)->info.fill_ratio_needed = 1.0;
392: b = (Mat_SeqBAIJ*)(fact)->data;
393: b->row = isrow;
394: b->col = iscol;
395: b->icol = isicol;
396: PetscObjectReference((PetscObject)isrow);
397: PetscObjectReference((PetscObject)iscol);
398: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
400: PetscMalloc1((n+1)*bs,&b->solve_work);
401: return(0);
402: }
404: ISGetIndices(isrow,&r);
405: ISGetIndices(isicol,&ic);
407: /* get new row pointers */
408: PetscMalloc1(n+1,&bi);
409: bi[0] = 0;
410: /* bdiag is location of diagonal in factor */
411: PetscMalloc1(n+1,&bdiag);
412: bdiag[0] = 0;
414: PetscMalloc2(n,&bj_ptr,n,&bjlvl_ptr);
416: /* create a linked list for storing column indices of the active row */
417: nlnk = n + 1;
418: PetscIncompleteLLCreate(n,n,nlnk,lnk,lnk_lvl,lnkbt);
420: /* initial FreeSpace size is f*(ai[n]+1) */
421: PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);
422: current_space = free_space;
423: PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space_lvl);
424: current_space_lvl = free_space_lvl;
426: for (i=0; i<n; i++) {
427: nzi = 0;
428: /* copy current row into linked list */
429: nnz = ai[r[i]+1] - ai[r[i]];
430: 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);
431: cols = aj + ai[r[i]];
432: lnk[i] = -1; /* marker to indicate if diagonal exists */
433: PetscIncompleteLLInit(nnz,cols,n,ic,nlnk,lnk,lnk_lvl,lnkbt);
434: nzi += nlnk;
436: /* make sure diagonal entry is included */
437: if (diagonal_fill && lnk[i] == -1) {
438: fm = n;
439: while (lnk[fm] < i) fm = lnk[fm];
440: lnk[i] = lnk[fm]; /* insert diagonal into linked list */
441: lnk[fm] = i;
442: lnk_lvl[i] = 0;
443: nzi++; dcount++;
444: }
446: /* add pivot rows into the active row */
447: nzbd = 0;
448: prow = lnk[n];
449: while (prow < i) {
450: nnz = bdiag[prow];
451: cols = bj_ptr[prow] + nnz + 1;
452: cols_lvl = bjlvl_ptr[prow] + nnz + 1;
453: nnz = bi[prow+1] - bi[prow] - nnz - 1;
455: PetscILULLAddSorted(nnz,cols,levels,cols_lvl,prow,nlnk,lnk,lnk_lvl,lnkbt,prow);
456: nzi += nlnk;
457: prow = lnk[prow];
458: nzbd++;
459: }
460: bdiag[i] = nzbd;
461: bi[i+1] = bi[i] + nzi;
463: /* if free space is not available, make more free space */
464: if (current_space->local_remaining<nzi) {
465: nnz = PetscIntMultTruncate(2,PetscIntMultTruncate(nzi,(n - i))); /* estimated and max additional space needed */
466: PetscFreeSpaceGet(nnz,¤t_space);
467: PetscFreeSpaceGet(nnz,¤t_space_lvl);
468: reallocs++;
469: }
471: /* copy data into free_space and free_space_lvl, then initialize lnk */
472: PetscIncompleteLLClean(n,n,nzi,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);
474: bj_ptr[i] = current_space->array;
475: bjlvl_ptr[i] = current_space_lvl->array;
477: /* make sure the active row i has diagonal entry */
478: 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);
480: current_space->array += nzi;
481: current_space->local_used += nzi;
482: current_space->local_remaining -= nzi;
484: current_space_lvl->array += nzi;
485: current_space_lvl->local_used += nzi;
486: current_space_lvl->local_remaining -= nzi;
487: }
489: ISRestoreIndices(isrow,&r);
490: ISRestoreIndices(isicol,&ic);
492: /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
493: PetscMalloc1(bi[n]+1,&bj);
494: PetscFreeSpaceContiguous_LU(&free_space,bj,n,bi,bdiag);
496: PetscIncompleteLLDestroy(lnk,lnkbt);
497: PetscFreeSpaceDestroy(free_space_lvl);
498: PetscFree2(bj_ptr,bjlvl_ptr);
500: #if defined(PETSC_USE_INFO)
501: {
502: PetscReal af = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]);
503: PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);
504: PetscInfo1(A,"Run with -[sub_]pc_factor_fill %g or use \n",(double)af);
505: PetscInfo1(A,"PCFactorSetFill([sub]pc,%g);\n",(double)af);
506: PetscInfo(A,"for best performance.\n");
507: if (diagonal_fill) {
508: PetscInfo1(A,"Detected and replaced %D missing diagonals\n",dcount);
509: }
510: }
511: #endif
513: /* put together the new matrix */
514: MatSeqBAIJSetPreallocation(fact,bs,MAT_SKIP_ALLOCATION,NULL);
515: PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol);
517: b = (Mat_SeqBAIJ*)(fact)->data;
518: b->free_a = PETSC_TRUE;
519: b->free_ij = PETSC_TRUE;
520: b->singlemalloc = PETSC_FALSE;
522: PetscMalloc1(bs2*(bdiag[0]+1),&b->a);
524: b->j = bj;
525: b->i = bi;
526: b->diag = bdiag;
527: b->free_diag = PETSC_TRUE;
528: b->ilen = 0;
529: b->imax = 0;
530: b->row = isrow;
531: b->col = iscol;
532: PetscObjectReference((PetscObject)isrow);
533: PetscObjectReference((PetscObject)iscol);
534: b->icol = isicol;
536: PetscMalloc1(bs*n+bs,&b->solve_work);
537: /* In b structure: Free imax, ilen, old a, old j.
538: Allocate bdiag, solve_work, new a, new j */
539: PetscLogObjectMemory((PetscObject)fact,(bdiag[0]+1) * (sizeof(PetscInt)+bs2*sizeof(PetscScalar)));
540: b->maxnz = b->nz = bdiag[0]+1;
542: fact->info.factor_mallocs = reallocs;
543: fact->info.fill_ratio_given = f;
544: fact->info.fill_ratio_needed = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]);
546: MatSeqBAIJSetNumericFactorization(fact,both_identity);
547: return(0);
548: }
550: /*
551: This code is virtually identical to MatILUFactorSymbolic_SeqAIJ
552: except that the data structure of Mat_SeqAIJ is slightly different.
553: Not a good example of code reuse.
554: */
555: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_inplace(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
556: {
557: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b;
558: IS isicol;
560: const PetscInt *r,*ic,*ai = a->i,*aj = a->j,*xi;
561: PetscInt prow,n = a->mbs,*ainew,*ajnew,jmax,*fill,nz,*im,*ajfill,*flev,*xitmp;
562: PetscInt *dloc,idx,row,m,fm,nzf,nzi,reallocate = 0,dcount = 0;
563: PetscInt incrlev,nnz,i,bs = A->rmap->bs,bs2 = a->bs2,levels,diagonal_fill,dd;
564: PetscBool col_identity,row_identity,both_identity,flg;
565: PetscReal f;
568: MatMissingDiagonal_SeqBAIJ(A,&flg,&dd);
569: if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix A is missing diagonal entry in row %D",dd);
571: f = info->fill;
572: levels = (PetscInt)info->levels;
573: diagonal_fill = (PetscInt)info->diagonal_fill;
575: ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
577: ISIdentity(isrow,&row_identity);
578: ISIdentity(iscol,&col_identity);
579: both_identity = (PetscBool) (row_identity && col_identity);
581: if (!levels && both_identity) { /* special case copy the nonzero structure */
582: MatDuplicateNoCreate_SeqBAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_TRUE);
583: MatSeqBAIJSetNumericFactorization_inplace(fact,both_identity);
585: fact->factortype = MAT_FACTOR_ILU;
586: b = (Mat_SeqBAIJ*)fact->data;
587: b->row = isrow;
588: b->col = iscol;
589: PetscObjectReference((PetscObject)isrow);
590: PetscObjectReference((PetscObject)iscol);
591: b->icol = isicol;
592: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
594: PetscMalloc1((n+1)*bs,&b->solve_work);
595: return(0);
596: }
598: /* general case perform the symbolic factorization */
599: ISGetIndices(isrow,&r);
600: ISGetIndices(isicol,&ic);
602: /* get new row pointers */
603: PetscMalloc1(n+1,&ainew);
604: ainew[0] = 0;
605: /* don't know how many column pointers are needed so estimate */
606: jmax = (PetscInt)(f*ai[n] + 1);
607: PetscMalloc1(jmax,&ajnew);
608: /* ajfill is level of fill for each fill entry */
609: PetscMalloc1(jmax,&ajfill);
610: /* fill is a linked list of nonzeros in active row */
611: PetscMalloc1(n+1,&fill);
612: /* im is level for each filled value */
613: PetscMalloc1(n+1,&im);
614: /* dloc is location of diagonal in factor */
615: PetscMalloc1(n+1,&dloc);
616: dloc[0] = 0;
617: for (prow=0; prow<n; prow++) {
619: /* copy prow into linked list */
620: nzf = nz = ai[r[prow]+1] - ai[r[prow]];
621: 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);
622: xi = aj + ai[r[prow]];
623: fill[n] = n;
624: fill[prow] = -1; /* marker for diagonal entry */
625: while (nz--) {
626: fm = n;
627: idx = ic[*xi++];
628: do {
629: m = fm;
630: fm = fill[m];
631: } while (fm < idx);
632: fill[m] = idx;
633: fill[idx] = fm;
634: im[idx] = 0;
635: }
637: /* make sure diagonal entry is included */
638: if (diagonal_fill && fill[prow] == -1) {
639: fm = n;
640: while (fill[fm] < prow) fm = fill[fm];
641: fill[prow] = fill[fm]; /* insert diagonal into linked list */
642: fill[fm] = prow;
643: im[prow] = 0;
644: nzf++;
645: dcount++;
646: }
648: nzi = 0;
649: row = fill[n];
650: while (row < prow) {
651: incrlev = im[row] + 1;
652: nz = dloc[row];
653: xi = ajnew + ainew[row] + nz + 1;
654: flev = ajfill + ainew[row] + nz + 1;
655: nnz = ainew[row+1] - ainew[row] - nz - 1;
656: fm = row;
657: while (nnz-- > 0) {
658: idx = *xi++;
659: if (*flev + incrlev > levels) {
660: flev++;
661: continue;
662: }
663: do {
664: m = fm;
665: fm = fill[m];
666: } while (fm < idx);
667: if (fm != idx) {
668: im[idx] = *flev + incrlev;
669: fill[m] = idx;
670: fill[idx] = fm;
671: fm = idx;
672: nzf++;
673: } else if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev;
674: flev++;
675: }
676: row = fill[row];
677: nzi++;
678: }
679: /* copy new filled row into permanent storage */
680: ainew[prow+1] = ainew[prow] + nzf;
681: if (ainew[prow+1] > jmax) {
683: /* estimate how much additional space we will need */
684: /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
685: /* just double the memory each time */
686: PetscInt maxadd = jmax;
687: /* maxadd = (int)(((f*ai[n]+1)*(n-prow+5))/n); */
688: if (maxadd < nzf) maxadd = (n-prow)*(nzf+1);
689: jmax += maxadd;
691: /* allocate a longer ajnew and ajfill */
692: PetscMalloc1(jmax,&xitmp);
693: PetscMemcpy(xitmp,ajnew,ainew[prow]*sizeof(PetscInt));
694: PetscFree(ajnew);
695: ajnew = xitmp;
696: PetscMalloc1(jmax,&xitmp);
697: PetscMemcpy(xitmp,ajfill,ainew[prow]*sizeof(PetscInt));
698: PetscFree(ajfill);
699: ajfill = xitmp;
700: reallocate++; /* count how many reallocations are needed */
701: }
702: xitmp = ajnew + ainew[prow];
703: flev = ajfill + ainew[prow];
704: dloc[prow] = nzi;
705: fm = fill[n];
706: while (nzf--) {
707: *xitmp++ = fm;
708: *flev++ = im[fm];
709: fm = fill[fm];
710: }
711: /* make sure row has diagonal entry */
712: 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\
713: try running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill",prow);
714: }
715: PetscFree(ajfill);
716: ISRestoreIndices(isrow,&r);
717: ISRestoreIndices(isicol,&ic);
718: PetscFree(fill);
719: PetscFree(im);
721: #if defined(PETSC_USE_INFO)
722: {
723: PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
724: PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocate,(double)f,(double)af);
725: PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);
726: PetscInfo1(A,"PCFactorSetFill(pc,%g);\n",(double)af);
727: PetscInfo(A,"for best performance.\n");
728: if (diagonal_fill) {
729: PetscInfo1(A,"Detected and replaced %D missing diagonals\n",dcount);
730: }
731: }
732: #endif
734: /* put together the new matrix */
735: MatSeqBAIJSetPreallocation(fact,bs,MAT_SKIP_ALLOCATION,NULL);
736: PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol);
737: b = (Mat_SeqBAIJ*)fact->data;
739: b->free_a = PETSC_TRUE;
740: b->free_ij = PETSC_TRUE;
741: b->singlemalloc = PETSC_FALSE;
743: PetscMalloc1(bs2*ainew[n],&b->a);
745: b->j = ajnew;
746: b->i = ainew;
747: for (i=0; i<n; i++) dloc[i] += ainew[i];
748: b->diag = dloc;
749: b->free_diag = PETSC_TRUE;
750: b->ilen = 0;
751: b->imax = 0;
752: b->row = isrow;
753: b->col = iscol;
754: b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
756: PetscObjectReference((PetscObject)isrow);
757: PetscObjectReference((PetscObject)iscol);
758: b->icol = isicol;
759: PetscMalloc1(bs*n+bs,&b->solve_work);
760: /* In b structure: Free imax, ilen, old a, old j.
761: Allocate dloc, solve_work, new a, new j */
762: PetscLogObjectMemory((PetscObject)fact,(ainew[n]-n)*(sizeof(PetscInt))+bs2*ainew[n]*sizeof(PetscScalar));
763: b->maxnz = b->nz = ainew[n];
765: fact->info.factor_mallocs = reallocate;
766: fact->info.fill_ratio_given = f;
767: fact->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[prow]);
769: MatSeqBAIJSetNumericFactorization_inplace(fact,both_identity);
770: return(0);
771: }
773: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE(Mat A)
774: {
775: /* Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; */
776: /* int i,*AJ=a->j,nz=a->nz; */
779: /* Undo Column scaling */
780: /* while (nz--) { */
781: /* AJ[i] = AJ[i]/4; */
782: /* } */
783: /* This should really invoke a push/pop logic, but we don't have that yet. */
784: A->ops->setunfactored = NULL;
785: return(0);
786: }
788: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE_usj(Mat A)
789: {
790: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
791: PetscInt *AJ=a->j,nz=a->nz;
792: unsigned short *aj=(unsigned short*)AJ;
795: /* Is this really necessary? */
796: while (nz--) {
797: AJ[nz] = (int)((unsigned int)aj[nz]); /* First extend, then convert to signed. */
798: }
799: A->ops->setunfactored = NULL;
800: return(0);
801: }