Actual source code: baijfact2.c

petsc-master 2019-05-24
Report Typos and Errors

  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,&current_space);
467:       PetscFreeSpaceGet(nnz,&current_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: }