Actual source code: baijfact.c

      1: #define PETSCMAT_DLL

      3: /*
      4:     Factorization code for BAIJ format.
      5: */
     6:  #include src/mat/impls/baij/seq/baij.h
     7:  #include src/inline/ilu.h

      9: /* ------------------------------------------------------------*/
     10: /*
     11:       Version for when blocks are 2 by 2
     12: */
     15: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2(Mat A,MatFactorInfo *info,Mat *B)
     16: {
     17:   Mat            C = *B;
     18:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ *)C->data;
     19:   IS             isrow = b->row,isicol = b->icol;
     21:   PetscInt       *r,*ic,i,j,n = a->mbs,*bi = b->i,*bj = b->j;
     22:   PetscInt       *ajtmpold,*ajtmp,nz,row;
     23:   PetscInt       *diag_offset=b->diag,idx,*ai=a->i,*aj=a->j,*pj;
     24:   MatScalar      *pv,*v,*rtmp,m1,m2,m3,m4,*pc,*w,*x,x1,x2,x3,x4;
     25:   MatScalar      p1,p2,p3,p4;
     26:   MatScalar      *ba = b->a,*aa = a->a;

     29:   ISGetIndices(isrow,&r);
     30:   ISGetIndices(isicol,&ic);
     31:   PetscMalloc(4*(n+1)*sizeof(MatScalar),&rtmp);

     33:   for (i=0; i<n; i++) {
     34:     nz    = bi[i+1] - bi[i];
     35:     ajtmp = bj + bi[i];
     36:     for  (j=0; j<nz; j++) {
     37:       x = rtmp+4*ajtmp[j]; x[0] = x[1] = x[2] = x[3] = 0.0;
     38:     }
     39:     /* load in initial (unfactored row) */
     40:     idx      = r[i];
     41:     nz       = ai[idx+1] - ai[idx];
     42:     ajtmpold = aj + ai[idx];
     43:     v        = aa + 4*ai[idx];
     44:     for (j=0; j<nz; j++) {
     45:       x    = rtmp+4*ic[ajtmpold[j]];
     46:       x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3];
     47:       v    += 4;
     48:     }
     49:     row = *ajtmp++;
     50:     while (row < i) {
     51:       pc = rtmp + 4*row;
     52:       p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3];
     53:       if (p1 != 0.0 || p2 != 0.0 || p3 != 0.0 || p4 != 0.0) {
     54:         pv = ba + 4*diag_offset[row];
     55:         pj = bj + diag_offset[row] + 1;
     56:         x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3];
     57:         pc[0] = m1 = p1*x1 + p3*x2;
     58:         pc[1] = m2 = p2*x1 + p4*x2;
     59:         pc[2] = m3 = p1*x3 + p3*x4;
     60:         pc[3] = m4 = p2*x3 + p4*x4;
     61:         nz = bi[row+1] - diag_offset[row] - 1;
     62:         pv += 4;
     63:         for (j=0; j<nz; j++) {
     64:           x1   = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3];
     65:           x    = rtmp + 4*pj[j];
     66:           x[0] -= m1*x1 + m3*x2;
     67:           x[1] -= m2*x1 + m4*x2;
     68:           x[2] -= m1*x3 + m3*x4;
     69:           x[3] -= m2*x3 + m4*x4;
     70:           pv   += 4;
     71:         }
     72:         PetscLogFlops(16*nz+12);
     73:       }
     74:       row = *ajtmp++;
     75:     }
     76:     /* finished row so stick it into b->a */
     77:     pv = ba + 4*bi[i];
     78:     pj = bj + bi[i];
     79:     nz = bi[i+1] - bi[i];
     80:     for (j=0; j<nz; j++) {
     81:       x     = rtmp+4*pj[j];
     82:       pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3];
     83:       pv   += 4;
     84:     }
     85:     /* invert diagonal block */
     86:     w = ba + 4*diag_offset[i];
     87:     Kernel_A_gets_inverse_A_2(w);
     88:   }

     90:   PetscFree(rtmp);
     91:   ISRestoreIndices(isicol,&ic);
     92:   ISRestoreIndices(isrow,&r);
     93:   C->factor = FACTOR_LU;
     94:   C->assembled = PETSC_TRUE;
     95:   PetscLogFlops(1.3333*8*b->mbs); /* from inverting diagonal blocks */
     96:   return(0);
     97: }
     98: /*
     99:       Version for when blocks are 2 by 2 Using natural ordering
    100: */
    103: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering(Mat A,MatFactorInfo *info,Mat *B)
    104: {
    105:   Mat            C = *B;
    106:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ *)C->data;
    108:   PetscInt       i,j,n = a->mbs,*bi = b->i,*bj = b->j;
    109:   PetscInt       *ajtmpold,*ajtmp,nz,row;
    110:   PetscInt       *diag_offset = b->diag,*ai=a->i,*aj=a->j,*pj;
    111:   MatScalar      *pv,*v,*rtmp,*pc,*w,*x;
    112:   MatScalar      p1,p2,p3,p4,m1,m2,m3,m4,x1,x2,x3,x4;
    113:   MatScalar      *ba = b->a,*aa = a->a;

    116:   PetscMalloc(4*(n+1)*sizeof(MatScalar),&rtmp);

    118:   for (i=0; i<n; i++) {
    119:     nz    = bi[i+1] - bi[i];
    120:     ajtmp = bj + bi[i];
    121:     for  (j=0; j<nz; j++) {
    122:       x = rtmp+4*ajtmp[j];
    123:       x[0]  = x[1]  = x[2]  = x[3]  = 0.0;
    124:     }
    125:     /* load in initial (unfactored row) */
    126:     nz       = ai[i+1] - ai[i];
    127:     ajtmpold = aj + ai[i];
    128:     v        = aa + 4*ai[i];
    129:     for (j=0; j<nz; j++) {
    130:       x    = rtmp+4*ajtmpold[j];
    131:       x[0]  = v[0];  x[1]  = v[1];  x[2]  = v[2];  x[3]  = v[3];
    132:       v    += 4;
    133:     }
    134:     row = *ajtmp++;
    135:     while (row < i) {
    136:       pc  = rtmp + 4*row;
    137:       p1  = pc[0];  p2  = pc[1];  p3  = pc[2];  p4  = pc[3];
    138:       if (p1 != 0.0 || p2 != 0.0 || p3 != 0.0 || p4 != 0.0) {
    139:         pv = ba + 4*diag_offset[row];
    140:         pj = bj + diag_offset[row] + 1;
    141:         x1  = pv[0];  x2  = pv[1];  x3  = pv[2];  x4  = pv[3];
    142:         pc[0] = m1 = p1*x1 + p3*x2;
    143:         pc[1] = m2 = p2*x1 + p4*x2;
    144:         pc[2] = m3 = p1*x3 + p3*x4;
    145:         pc[3] = m4 = p2*x3 + p4*x4;
    146:         nz = bi[row+1] - diag_offset[row] - 1;
    147:         pv += 4;
    148:         for (j=0; j<nz; j++) {
    149:           x1   = pv[0];  x2  = pv[1];   x3 = pv[2];  x4  = pv[3];
    150:           x    = rtmp + 4*pj[j];
    151:           x[0] -= m1*x1 + m3*x2;
    152:           x[1] -= m2*x1 + m4*x2;
    153:           x[2] -= m1*x3 + m3*x4;
    154:           x[3] -= m2*x3 + m4*x4;
    155:           pv   += 4;
    156:         }
    157:         PetscLogFlops(16*nz+12);
    158:       }
    159:       row = *ajtmp++;
    160:     }
    161:     /* finished row so stick it into b->a */
    162:     pv = ba + 4*bi[i];
    163:     pj = bj + bi[i];
    164:     nz = bi[i+1] - bi[i];
    165:     for (j=0; j<nz; j++) {
    166:       x      = rtmp+4*pj[j];
    167:       pv[0]  = x[0];  pv[1]  = x[1];  pv[2]  = x[2];  pv[3]  = x[3];
    168:       pv   += 4;
    169:     }
    170:     /* invert diagonal block */
    171:     w = ba + 4*diag_offset[i];
    172:     Kernel_A_gets_inverse_A_2(w);
    173:     /*Kernel_A_gets_inverse_A(bs,w,v_pivots,v_work);*/
    174:   }

    176:   PetscFree(rtmp);
    177:   C->factor    = FACTOR_LU;
    178:   C->assembled = PETSC_TRUE;
    179:   PetscLogFlops(1.3333*8*b->mbs); /* from inverting diagonal blocks */
    180:   return(0);
    181: }

    183: /* ----------------------------------------------------------- */
    184: /*
    185:      Version for when blocks are 1 by 1.
    186: */
    189: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1(Mat A,MatFactorInfo *info,Mat *B)
    190: {
    191:   Mat            C = *B;
    192:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ *)C->data;
    193:   IS             isrow = b->row,isicol = b->icol;
    195:   PetscInt       *r,*ic,i,j,n = a->mbs,*bi = b->i,*bj = b->j;
    196:   PetscInt       *ajtmpold,*ajtmp,nz,row,*ai = a->i,*aj = a->j;
    197:   PetscInt       *diag_offset = b->diag,diag,*pj;
    198:   MatScalar      *pv,*v,*rtmp,multiplier,*pc;
    199:   MatScalar      *ba = b->a,*aa = a->a;

    202:   ISGetIndices(isrow,&r);
    203:   ISGetIndices(isicol,&ic);
    204:   PetscMalloc((n+1)*sizeof(MatScalar),&rtmp);

    206:   for (i=0; i<n; i++) {
    207:     nz    = bi[i+1] - bi[i];
    208:     ajtmp = bj + bi[i];
    209:     for  (j=0; j<nz; j++) rtmp[ajtmp[j]] = 0.0;

    211:     /* load in initial (unfactored row) */
    212:     nz       = ai[r[i]+1] - ai[r[i]];
    213:     ajtmpold = aj + ai[r[i]];
    214:     v        = aa + ai[r[i]];
    215:     for (j=0; j<nz; j++) rtmp[ic[ajtmpold[j]]] =  v[j];

    217:     row = *ajtmp++;
    218:     while (row < i) {
    219:       pc = rtmp + row;
    220:       if (*pc != 0.0) {
    221:         pv         = ba + diag_offset[row];
    222:         pj         = bj + diag_offset[row] + 1;
    223:         multiplier = *pc * *pv++;
    224:         *pc        = multiplier;
    225:         nz         = bi[row+1] - diag_offset[row] - 1;
    226:         for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j];
    227:         PetscLogFlops(1+2*nz);
    228:       }
    229:       row = *ajtmp++;
    230:     }
    231:     /* finished row so stick it into b->a */
    232:     pv = ba + bi[i];
    233:     pj = bj + bi[i];
    234:     nz = bi[i+1] - bi[i];
    235:     for (j=0; j<nz; j++) {pv[j] = rtmp[pj[j]];}
    236:     diag = diag_offset[i] - bi[i];
    237:     /* check pivot entry for current row */
    238:     if (pv[diag] == 0.0) {
    239:       SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot");
    240:     }
    241:     pv[diag] = 1.0/pv[diag];
    242:   }

    244:   PetscFree(rtmp);
    245:   ISRestoreIndices(isicol,&ic);
    246:   ISRestoreIndices(isrow,&r);
    247:   C->factor    = FACTOR_LU;
    248:   C->assembled = PETSC_TRUE;
    249:   PetscLogFlops(C->n);
    250:   return(0);
    251: }

    254: /* ----------------------------------------------------------- */
    257: PetscErrorCode MatLUFactor_SeqBAIJ(Mat A,IS row,IS col,MatFactorInfo *info)
    258: {
    260:   Mat            C;

    263:   MatLUFactorSymbolic(A,row,col,info,&C);
    264:   MatLUFactorNumeric(A,info,&C);
    265:   MatHeaderCopy(A,C);
    266:   PetscLogObjectParent(A,((Mat_SeqBAIJ*)(A->data))->icol);
    267:   return(0);
    268: }

    270:  #include src/mat/impls/sbaij/seq/sbaij.h
    273: PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N(Mat A,MatFactorInfo *info,Mat *B)
    274: {
    276:   Mat            C = *B;
    277:   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data;
    278:   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
    279:   IS             ip=b->row;
    280:   PetscInt       *rip,i,j,mbs=a->mbs,bs=A->bs,*bi=b->i,*bj=b->j,*bcol;
    281:   PetscInt       *ai=a->i,*aj=a->j;
    282:   PetscInt       k,jmin,jmax,*jl,*il,col,nexti,ili,nz;
    283:   MatScalar      *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
    284:   PetscReal      zeropivot,rs,shiftnz;
    285:   PetscTruth     shiftpd;
    286:   ChShift_Ctx    sctx;
    287:   PetscInt       newshift;

    290:   if (bs > 1) {
    291:     if (!a->sbaijMat){
    292:       MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
    293:     }
    294:     (a->sbaijMat)->ops->choleskyfactornumeric(a->sbaijMat,info,B);
    295:     MatDestroy(a->sbaijMat);
    296:     a->sbaijMat = PETSC_NULL;
    297:     return(0);
    298:   }
    299:
    300:   /* initialization */
    301:   shiftnz   = info->shiftnz;
    302:   shiftpd   = info->shiftpd;
    303:   zeropivot = info->zeropivot;

    305:   ISGetIndices(ip,&rip);
    306:   nz   = (2*mbs+1)*sizeof(PetscInt)+mbs*sizeof(MatScalar);
    307:   PetscMalloc(nz,&il);
    308:   jl   = il + mbs;
    309:   rtmp = (MatScalar*)(jl + mbs);

    311:   sctx.shift_amount = 0;
    312:   sctx.nshift       = 0;
    313:   do {
    314:     sctx.chshift = PETSC_FALSE;
    315:     for (i=0; i<mbs; i++) {
    316:       rtmp[i] = 0.0; jl[i] = mbs; il[0] = 0;
    317:     }
    318:
    319:     for (k = 0; k<mbs; k++){
    320:       bval = ba + bi[k];
    321:       /* initialize k-th row by the perm[k]-th row of A */
    322:       jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
    323:       for (j = jmin; j < jmax; j++){
    324:         col = rip[aj[j]];
    325:         if (col >= k){ /* only take upper triangular entry */
    326:           rtmp[col] = aa[j];
    327:           *bval++  = 0.0; /* for in-place factorization */
    328:         }
    329:       }
    330:
    331:       /* shift the diagonal of the matrix */
    332:       if (sctx.nshift) rtmp[k] += sctx.shift_amount;

    334:       /* modify k-th row by adding in those rows i with U(i,k)!=0 */
    335:       dk = rtmp[k];
    336:       i = jl[k]; /* first row to be added to k_th row  */

    338:       while (i < k){
    339:         nexti = jl[i]; /* next row to be added to k_th row */

    341:         /* compute multiplier, update diag(k) and U(i,k) */
    342:         ili = il[i];  /* index of first nonzero element in U(i,k:bms-1) */
    343:         uikdi = - ba[ili]*ba[bi[i]];  /* diagonal(k) */
    344:         dk += uikdi*ba[ili];
    345:         ba[ili] = uikdi; /* -U(i,k) */

    347:         /* add multiple of row i to k-th row */
    348:         jmin = ili + 1; jmax = bi[i+1];
    349:         if (jmin < jmax){
    350:           for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];
    351:           /* update il and jl for row i */
    352:           il[i] = jmin;
    353:           j = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
    354:         }
    355:         i = nexti;
    356:       }

    358:       /* shift the diagonals when zero pivot is detected */
    359:       /* compute rs=sum of abs(off-diagonal) */
    360:       rs   = 0.0;
    361:       jmin = bi[k]+1;
    362:       nz   = bi[k+1] - jmin;
    363:       if (nz){
    364:         bcol = bj + jmin;
    365:         while (nz--){
    366:           rs += PetscAbsScalar(rtmp[*bcol]);
    367:           bcol++;
    368:         }
    369:       }

    371:       sctx.rs = rs;
    372:       sctx.pv = dk;
    373:       MatCholeskyCheckShift_inline(info,sctx,newshift);
    374:       if (newshift == 1){
    375:         break;    /* sctx.shift_amount is updated */
    376:       } else if (newshift == -1){
    377:         SETERRQ4(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %D value %g tolerance %g * rs %g",k,PetscAbsScalar(dk),zeropivot,rs);
    378:       }

    380:       /* copy data into U(k,:) */
    381:       ba[bi[k]] = 1.0/dk; /* U(k,k) */
    382:       jmin = bi[k]+1; jmax = bi[k+1];
    383:       if (jmin < jmax) {
    384:         for (j=jmin; j<jmax; j++){
    385:           col = bj[j]; ba[j] = rtmp[col]; rtmp[col] = 0.0;
    386:         }
    387:         /* add the k-th row into il and jl */
    388:         il[k] = jmin;
    389:         i = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
    390:       }
    391:     }
    392:   } while (sctx.chshift);
    393:   PetscFree(il);

    395:   ISRestoreIndices(ip,&rip);
    396:   C->factor       = FACTOR_CHOLESKY;
    397:   C->assembled    = PETSC_TRUE;
    398:   C->preallocated = PETSC_TRUE;
    399:   PetscLogFlops(C->m);
    400:   if (sctx.nshift){
    401:     if (shiftnz) {
    402:       PetscLogInfo((0,"MatCholeskyFactorNumeric_SeqBAIJ_1: number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,sctx.shift_amount));
    403:     } else if (shiftpd) {
    404:       PetscLogInfo((0,"MatCholeskyFactorNumeric_SeqBAIJ_1: number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,sctx.shift_amount));
    405:     }
    406:   }
    407:   return(0);
    408: }

    412: PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering(Mat A,MatFactorInfo *info,Mat *fact)
    413: {
    414:   Mat            C = *fact;
    415:   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data;
    416:   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
    418:   PetscInt       i,j,am=a->mbs;
    419:   PetscInt       *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
    420:   PetscInt       k,jmin,*jl,*il,nexti,ili,*acol,*bcol,nz;
    421:   MatScalar      *rtmp,*ba=b->a,*aa=a->a,dk,uikdi,*aval,*bval;
    422:   PetscReal      zeropivot,rs,shiftnz;
    423:   PetscTruth     shiftpd;
    424:   ChShift_Ctx    sctx;
    425:   PetscInt       newshift;

    428:   /* initialization */
    429:   shiftnz   = info->shiftnz;
    430:   shiftpd   = info->shiftpd;
    431:   zeropivot = info->zeropivot;

    433:   nz   = (2*am+1)*sizeof(PetscInt)+am*sizeof(MatScalar);
    434:   PetscMalloc(nz,&il);
    435:   jl   = il + am;
    436:   rtmp = (MatScalar*)(jl + am);

    438:   sctx.shift_amount = 0;
    439:   sctx.nshift       = 0;
    440:   do {
    441:     sctx.chshift = PETSC_FALSE;
    442:     for (i=0; i<am; i++) {
    443:       rtmp[i] = 0.0; jl[i] = am; il[0] = 0;
    444:     }

    446:     for (k = 0; k<am; k++){
    447:     /* initialize k-th row with elements nonzero in row perm(k) of A */
    448:       nz   = ai[k+1] - ai[k];
    449:       acol = aj + ai[k];
    450:       aval = aa + ai[k];
    451:       bval = ba + bi[k];
    452:       while (nz -- ){
    453:         if (*acol < k) { /* skip lower triangular entries */
    454:           acol++; aval++;
    455:         } else {
    456:           rtmp[*acol++] = *aval++;
    457:           *bval++       = 0.0; /* for in-place factorization */
    458:         }
    459:       }
    460:
    461:       /* shift the diagonal of the matrix */
    462:       if (sctx.nshift) rtmp[k] += sctx.shift_amount;
    463:
    464:       /* modify k-th row by adding in those rows i with U(i,k)!=0 */
    465:       dk = rtmp[k];
    466:       i  = jl[k]; /* first row to be added to k_th row  */

    468:       while (i < k){
    469:         nexti = jl[i]; /* next row to be added to k_th row */
    470:         /* compute multiplier, update D(k) and U(i,k) */
    471:         ili   = il[i];  /* index of first nonzero element in U(i,k:bms-1) */
    472:         uikdi = - ba[ili]*ba[bi[i]];
    473:         dk   += uikdi*ba[ili];
    474:         ba[ili] = uikdi; /* -U(i,k) */

    476:         /* add multiple of row i to k-th row ... */
    477:         jmin = ili + 1;
    478:         nz   = bi[i+1] - jmin;
    479:         if (nz > 0){
    480:           bcol = bj + jmin;
    481:           bval = ba + jmin;
    482:           while (nz --) rtmp[*bcol++] += uikdi*(*bval++);
    483:           /* update il and jl for i-th row */
    484:           il[i] = jmin;
    485:           j = bj[jmin]; jl[i] = jl[j]; jl[j] = i;
    486:         }
    487:         i = nexti;
    488:       }

    490:       /* shift the diagonals when zero pivot is detected */
    491:       /* compute rs=sum of abs(off-diagonal) */
    492:       rs   = 0.0;
    493:       jmin = bi[k]+1;
    494:       nz   = bi[k+1] - jmin;
    495:       if (nz){
    496:         bcol = bj + jmin;
    497:         while (nz--){
    498:           rs += PetscAbsScalar(rtmp[*bcol]);
    499:           bcol++;
    500:         }
    501:       }

    503:       sctx.rs = rs;
    504:       sctx.pv = dk;
    505:       MatCholeskyCheckShift_inline(info,sctx,newshift);
    506:       if (newshift == 1){
    507:         break;    /* sctx.shift_amount is updated */
    508:       } else if (newshift == -1){
    509:         SETERRQ4(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %D value %g tolerance %g * rs %g",k,PetscAbsScalar(dk),zeropivot,rs);
    510:       }

    512:       /* copy data into U(k,:) */
    513:       ba[bi[k]] = 1.0/dk;
    514:       jmin      = bi[k]+1;
    515:       nz        = bi[k+1] - jmin;
    516:       if (nz){
    517:         bcol = bj + jmin;
    518:         bval = ba + jmin;
    519:         while (nz--){
    520:           *bval++       = rtmp[*bcol];
    521:           rtmp[*bcol++] = 0.0;
    522:         }
    523:         /* add k-th row into il and jl */
    524:         il[k] = jmin;
    525:         i = bj[jmin]; jl[k] = jl[i]; jl[i] = k;
    526:       }
    527:     }
    528:   } while (sctx.chshift);
    529:   PetscFree(il);
    530:
    531:   C->factor       = FACTOR_CHOLESKY;
    532:   C->assembled    = PETSC_TRUE;
    533:   C->preallocated = PETSC_TRUE;
    534:   PetscLogFlops(C->m);
    535:     if (sctx.nshift){
    536:     if (shiftnz) {
    537:       PetscLogInfo((0,"MatCholeskyFactorNumeric_SeqBAIJ_1_NaturalOrdering: number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,sctx.shift_amount));
    538:     } else if (shiftpd) {
    539:       PetscLogInfo((0,"MatCholeskyFactorNumeric_SeqBAIJ_1_NaturalOrdering: number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,sctx.shift_amount));
    540:     }
    541:   }
    542:   return(0);
    543: }

    545:  #include petscbt.h
    546:  #include src/mat/utils/freespace.h
    549: PetscErrorCode MatICCFactorSymbolic_SeqBAIJ(Mat A,IS perm,MatFactorInfo *info,Mat *fact)
    550: {
    551:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
    552:   Mat_SeqSBAIJ   *b;
    553:   Mat            B;
    555:   PetscTruth     perm_identity;
    556:   PetscInt       reallocs=0,*rip,i,*ai=a->i,*aj=a->j,am=a->mbs,bs=A->bs,*ui;
    557:   PetscInt       jmin,jmax,nzk,k,j,*jl,prow,*il,nextprow;
    558:   PetscInt       nlnk,*lnk,*lnk_lvl=PETSC_NULL,ncols,ncols_upper,*cols,*cols_lvl,*uj,**uj_ptr,**uj_lvl_ptr;
    559:   PetscReal      fill=info->fill,levels=info->levels;
    560:   FreeSpaceList  free_space=PETSC_NULL,current_space=PETSC_NULL;
    561:   FreeSpaceList  free_space_lvl=PETSC_NULL,current_space_lvl=PETSC_NULL;
    562:   PetscBT        lnkbt;

    565:   if (bs > 1){
    566:     if (!a->sbaijMat){
    567:       MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
    568:     }
    569:     MatICCFactorSymbolic(a->sbaijMat,perm,info,fact);
    570:     B = *fact;
    571:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
    572:     return(0);
    573:   }

    575:   ISIdentity(perm,&perm_identity);
    576:   ISGetIndices(perm,&rip);

    578:   /* special case that simply copies fill pattern */
    579:   if (!levels && perm_identity) {
    580:     MatMarkDiagonal_SeqBAIJ(A);
    581:     PetscMalloc((am+1)*sizeof(PetscInt),&ui);
    582:     for (i=0; i<am; i++) {
    583:       ui[i] = ai[i+1] - a->diag[i]; /* ui: rowlengths - changes when !perm_identity */
    584:     }
    585:     MatCreate(PETSC_COMM_SELF,fact);
    586:     MatSetSizes(*fact,am,am,am,am);
    587:     B = *fact;
    588:     MatSetType(B,MATSEQSBAIJ);
    589:     MatSeqSBAIJSetPreallocation(B,1,0,ui);

    591:     b  = (Mat_SeqSBAIJ*)B->data;
    592:     uj = b->j;
    593:     for (i=0; i<am; i++) {
    594:       aj = a->j + a->diag[i];
    595:       for (j=0; j<ui[i]; j++){
    596:         *uj++ = *aj++;
    597:       }
    598:       b->ilen[i] = ui[i];
    599:     }
    600:     PetscFree(ui);
    601:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    602:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

    604:     B->ops->solve                 = MatSolve_SeqSBAIJ_1_NaturalOrdering;
    605:     B->ops->solvetranspose        = MatSolve_SeqSBAIJ_1_NaturalOrdering;
    606:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
    607:     return(0);
    608:   }

    610:   /* initialization */
    611:   PetscMalloc((am+1)*sizeof(PetscInt),&ui);
    612:   ui[0] = 0;
    613:   PetscMalloc((2*am+1)*sizeof(PetscInt),&cols_lvl);

    615:   /* jl: linked list for storing indices of the pivot rows
    616:      il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */
    617:   PetscMalloc((2*am+1)*sizeof(PetscInt)+2*am*sizeof(PetscInt*),&jl);
    618:   il         = jl + am;
    619:   uj_ptr     = (PetscInt**)(il + am);
    620:   uj_lvl_ptr = (PetscInt**)(uj_ptr + am);
    621:   for (i=0; i<am; i++){
    622:     jl[i] = am; il[i] = 0;
    623:   }

    625:   /* create and initialize a linked list for storing column indices of the active row k */
    626:   nlnk = am + 1;
    627:   PetscIncompleteLLCreate(am,am,nlnk,lnk,lnk_lvl,lnkbt);

    629:   /* initial FreeSpace size is fill*(ai[am]+1) */
    630:   GetMoreSpace((PetscInt)(fill*(ai[am]+1)),&free_space);
    631:   current_space = free_space;
    632:   GetMoreSpace((PetscInt)(fill*(ai[am]+1)),&free_space_lvl);
    633:   current_space_lvl = free_space_lvl;

    635:   for (k=0; k<am; k++){  /* for each active row k */
    636:     /* initialize lnk by the column indices of row rip[k] of A */
    637:     nzk   = 0;
    638:     ncols = ai[rip[k]+1] - ai[rip[k]];
    639:     ncols_upper = 0;
    640:     cols        = cols_lvl + am;
    641:     for (j=0; j<ncols; j++){
    642:       i = rip[*(aj + ai[rip[k]] + j)];
    643:       if (i >= k){ /* only take upper triangular entry */
    644:         cols[ncols_upper] = i;
    645:         cols_lvl[ncols_upper] = -1;  /* initialize level for nonzero entries */
    646:         ncols_upper++;
    647:       }
    648:     }
    649:     PetscIncompleteLLAdd(ncols_upper,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);
    650:     nzk += nlnk;

    652:     /* update lnk by computing fill-in for each pivot row to be merged in */
    653:     prow = jl[k]; /* 1st pivot row */
    654:
    655:     while (prow < k){
    656:       nextprow = jl[prow];
    657:
    658:       /* merge prow into k-th row */
    659:       jmin = il[prow] + 1;  /* index of the 2nd nzero entry in U(prow,k:am-1) */
    660:       jmax = ui[prow+1];
    661:       ncols = jmax-jmin;
    662:       i     = jmin - ui[prow];
    663:       cols = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */
    664:       for (j=0; j<ncols; j++) cols_lvl[j] = *(uj_lvl_ptr[prow] + i + j);
    665:       PetscIncompleteLLAddSorted(ncols,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);
    666:       nzk += nlnk;

    668:       /* update il and jl for prow */
    669:       if (jmin < jmax){
    670:         il[prow] = jmin;
    671:         j = *cols; jl[prow] = jl[j]; jl[j] = prow;
    672:       }
    673:       prow = nextprow;
    674:     }

    676:     /* if free space is not available, make more free space */
    677:     if (current_space->local_remaining<nzk) {
    678:       i = am - k + 1; /* num of unfactored rows */
    679:       i = PetscMin(i*nzk, i*(i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
    680:       GetMoreSpace(i,&current_space);
    681:       GetMoreSpace(i,&current_space_lvl);
    682:       reallocs++;
    683:     }

    685:     /* copy data into free_space and free_space_lvl, then initialize lnk */
    686:     PetscIncompleteLLClean(am,am,nzk,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);

    688:     /* add the k-th row into il and jl */
    689:     if (nzk-1 > 0){
    690:       i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */
    691:       jl[k] = jl[i]; jl[i] = k;
    692:       il[k] = ui[k] + 1;
    693:     }
    694:     uj_ptr[k]     = current_space->array;
    695:     uj_lvl_ptr[k] = current_space_lvl->array;

    697:     current_space->array           += nzk;
    698:     current_space->local_used      += nzk;
    699:     current_space->local_remaining -= nzk;

    701:     current_space_lvl->array           += nzk;
    702:     current_space_lvl->local_used      += nzk;
    703:     current_space_lvl->local_remaining -= nzk;

    705:     ui[k+1] = ui[k] + nzk;
    706:   }

    708: #if defined(PETSC_USE_DEBUG)
    709:   if (ai[am] != 0) {
    710:     PetscReal af = ((PetscReal)(2*ui[am]-am))/((PetscReal)ai[am]);
    711:     PetscLogInfo((A,"MatICCFactorSymbolic_SeqBAIJ:Reallocs %D Fill ratio:given %g needed %g\n",reallocs,fill,af));
    712:     PetscLogInfo((A,"MatICCFactorSymbolic_SeqBAIJ:Run with -pc_cholesky_fill %g or use \n",af));
    713:     PetscLogInfo((A,"MatICCFactorSymbolic_SeqBAIJ:PCCholeskySetFill(pc,%g) for best performance.\n",af));
    714:   } else {
    715:     PetscLogInfo((A,"MatICCFactorSymbolic_SeqBAIJ:Empty matrix.\n"));
    716:   }
    717: #endif

    719:   ISRestoreIndices(perm,&rip);
    720:   PetscFree(jl);
    721:   PetscFree(cols_lvl);

    723:   /* destroy list of free space and other temporary array(s) */
    724:   PetscMalloc((ui[am]+1)*sizeof(PetscInt),&uj);
    725:   MakeSpaceContiguous(&free_space,uj);
    726:   PetscIncompleteLLDestroy(lnk,lnkbt);
    727:   DestroySpace(free_space_lvl);

    729:   /* put together the new matrix in MATSEQSBAIJ format */
    730:   MatCreate(PETSC_COMM_SELF,fact);
    731:   MatSetSizes(*fact,am,am,am,am);
    732:   B = *fact;
    733:   MatSetType(B,MATSEQSBAIJ);
    734:   MatSeqSBAIJSetPreallocation(B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);

    736:   b = (Mat_SeqSBAIJ*)B->data;
    737:   b->singlemalloc = PETSC_FALSE;
    738:   PetscMalloc((ui[am]+1)*sizeof(MatScalar),&b->a);
    739:   b->j    = uj;
    740:   b->i    = ui;
    741:   b->diag = 0;
    742:   b->ilen = 0;
    743:   b->imax = 0;
    744:   b->row  = perm;
    745:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
    746:   PetscObjectReference((PetscObject)perm);
    747:   b->icol = perm;
    748:   PetscObjectReference((PetscObject)perm);
    749:   PetscMalloc((am+1)*sizeof(PetscScalar),&b->solve_work);
    750:   PetscLogObjectMemory(B,(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));
    751:   b->maxnz = b->nz = ui[am];
    752:
    753:   B->factor                 = FACTOR_CHOLESKY;
    754:   B->info.factor_mallocs    = reallocs;
    755:   B->info.fill_ratio_given  = fill;
    756:   if (ai[am] != 0) {
    757:     B->info.fill_ratio_needed = ((PetscReal)ui[am])/((PetscReal)ai[am]);
    758:   } else {
    759:     B->info.fill_ratio_needed = 0.0;
    760:   }
    761:   if (perm_identity){
    762:     B->ops->solve           = MatSolve_SeqSBAIJ_1_NaturalOrdering;
    763:     B->ops->solvetranspose  = MatSolve_SeqSBAIJ_1_NaturalOrdering;
    764:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
    765:   } else {
    766:     (*fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
    767:   }
    768:   return(0);
    769: }

    773: PetscErrorCode MatCholeskyFactorSymbolic_SeqBAIJ(Mat A,IS perm,MatFactorInfo *info,Mat *fact)
    774: {
    775:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
    776:   Mat_SeqSBAIJ   *b;
    777:   Mat            B;
    779:   PetscTruth     perm_identity;
    780:   PetscReal      fill = info->fill;
    781:   PetscInt       *rip,*riip,i,mbs=a->mbs,bs=A->bs,*ai=a->i,*aj=a->j,reallocs=0,prow;
    782:   PetscInt       *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow;
    783:   PetscInt       nlnk,*lnk,ncols,ncols_upper,*cols,*uj,**ui_ptr,*uj_ptr;
    784:   FreeSpaceList  free_space=PETSC_NULL,current_space=PETSC_NULL;
    785:   PetscBT        lnkbt;
    786:   IS             iperm;

    789:   if (bs > 1) { /* convert to seqsbaij */
    790:     if (!a->sbaijMat){
    791:       MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);
    792:     }
    793:     MatCholeskyFactorSymbolic(a->sbaijMat,perm,info,fact);
    794:     B    = *fact;
    795:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
    796:     return(0);
    797:   }

    799:   /* check whether perm is the identity mapping */
    800:   ISIdentity(perm,&perm_identity);
    801:   ISGetIndices(perm,&rip);

    803:   if (!perm_identity){
    804:     /* check if perm is symmetric! */
    805:     ISInvertPermutation(perm,PETSC_DECIDE,&iperm);
    806:     ISGetIndices(iperm,&riip);
    807:     for (i=0; i<mbs; i++) {
    808:       if (rip[i] != riip[i]) SETERRQ(PETSC_ERR_ARG_INCOMP,"Non-symmetric permutation, must use symmetric permutation");
    809:     }
    810:     ISRestoreIndices(iperm,&riip);
    811:     ISDestroy(iperm);
    812:   }

    814:   /* initialization */
    815:   PetscMalloc((mbs+1)*sizeof(PetscInt),&ui);
    816:   ui[0] = 0;

    818:   /* jl: linked list for storing indices of the pivot rows
    819:      il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */
    820:   PetscMalloc((3*mbs+1)*sizeof(PetscInt)+mbs*sizeof(PetscInt*),&jl);
    821:   il     = jl + mbs;
    822:   cols   = il + mbs;
    823:   ui_ptr = (PetscInt**)(cols + mbs);
    824:   for (i=0; i<mbs; i++){
    825:     jl[i] = mbs; il[i] = 0;
    826:   }

    828:   /* create and initialize a linked list for storing column indices of the active row k */
    829:   nlnk = mbs + 1;
    830:   PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);

    832:   /* initial FreeSpace size is fill*(ai[mbs]+1) */
    833:   GetMoreSpace((PetscInt)(fill*(ai[mbs]+1)),&free_space);
    834:   current_space = free_space;

    836:   for (k=0; k<mbs; k++){  /* for each active row k */
    837:     /* initialize lnk by the column indices of row rip[k] of A */
    838:     nzk   = 0;
    839:     ncols = ai[rip[k]+1] - ai[rip[k]];
    840:     ncols_upper = 0;
    841:     for (j=0; j<ncols; j++){
    842:       i = rip[*(aj + ai[rip[k]] + j)];
    843:       if (i >= k){ /* only take upper triangular entry */
    844:         cols[ncols_upper] = i;
    845:         ncols_upper++;
    846:       }
    847:     }
    848:     PetscLLAdd(ncols_upper,cols,mbs,nlnk,lnk,lnkbt);
    849:     nzk += nlnk;

    851:     /* update lnk by computing fill-in for each pivot row to be merged in */
    852:     prow = jl[k]; /* 1st pivot row */
    853:
    854:     while (prow < k){
    855:       nextprow = jl[prow];
    856:       /* merge prow into k-th row */
    857:       jmin = il[prow] + 1;  /* index of the 2nd nzero entry in U(prow,k:mbs-1) */
    858:       jmax = ui[prow+1];
    859:       ncols = jmax-jmin;
    860:       uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */
    861:       PetscLLAddSorted(ncols,uj_ptr,mbs,nlnk,lnk,lnkbt);
    862:       nzk += nlnk;

    864:       /* update il and jl for prow */
    865:       if (jmin < jmax){
    866:         il[prow] = jmin;
    867:         j = *uj_ptr; jl[prow] = jl[j]; jl[j] = prow;
    868:       }
    869:       prow = nextprow;
    870:     }

    872:     /* if free space is not available, make more free space */
    873:     if (current_space->local_remaining<nzk) {
    874:       i = mbs - k + 1; /* num of unfactored rows */
    875:       i = PetscMin(i*nzk, i*(i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */
    876:       GetMoreSpace(i,&current_space);
    877:       reallocs++;
    878:     }

    880:     /* copy data into free space, then initialize lnk */
    881:     PetscLLClean(mbs,mbs,nzk,lnk,current_space->array,lnkbt);

    883:     /* add the k-th row into il and jl */
    884:     if (nzk-1 > 0){
    885:       i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */
    886:       jl[k] = jl[i]; jl[i] = k;
    887:       il[k] = ui[k] + 1;
    888:     }
    889:     ui_ptr[k] = current_space->array;
    890:     current_space->array           += nzk;
    891:     current_space->local_used      += nzk;
    892:     current_space->local_remaining -= nzk;

    894:     ui[k+1] = ui[k] + nzk;
    895:   }

    897: #if defined(PETSC_USE_DEBUG)
    898:   if (ai[mbs] != 0) {
    899:     PetscReal af = ((PetscReal)ui[mbs])/((PetscReal)ai[mbs]);
    900:     PetscLogInfo((A,"MatCholeskyFactorSymbolic_SeqSBAIJ:Reallocs %D Fill ratio:given %g needed %g\n",reallocs,fill,af));
    901:     PetscLogInfo((A,"MatCholeskyFactorSymbolic_SeqSBAIJ:Run with -pc_cholesky_fill %g or use \n",af));
    902:     PetscLogInfo((A,"MatCholeskyFactorSymbolic_SeqSBAIJ:PCCholeskySetFill(pc,%g) for best performance.\n",af));
    903:   } else {
    904:     PetscLogInfo((A,"MatCholeskyFactorSymbolic_SeqSBAIJ:Empty matrix.\n"));
    905:   }
    906: #endif

    908:   ISRestoreIndices(perm,&rip);
    909:   PetscFree(jl);

    911:   /* destroy list of free space and other temporary array(s) */
    912:   PetscMalloc((ui[mbs]+1)*sizeof(PetscInt),&uj);
    913:   MakeSpaceContiguous(&free_space,uj);
    914:   PetscLLDestroy(lnk,lnkbt);

    916:   /* put together the new matrix in MATSEQSBAIJ format */
    917:   MatCreate(PETSC_COMM_SELF,fact);
    918:   MatSetSizes(*fact,mbs,mbs,mbs,mbs);
    919:   B    = *fact;
    920:   MatSetType(B,MATSEQSBAIJ);
    921:   MatSeqSBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,PETSC_NULL);

    923:   b = (Mat_SeqSBAIJ*)B->data;
    924:   b->singlemalloc = PETSC_FALSE;
    925:   PetscMalloc((ui[mbs]+1)*sizeof(MatScalar),&b->a);
    926:   b->j    = uj;
    927:   b->i    = ui;
    928:   b->diag = 0;
    929:   b->ilen = 0;
    930:   b->imax = 0;
    931:   b->row  = perm;
    932:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
    933:   PetscObjectReference((PetscObject)perm);
    934:   b->icol = perm;
    935:   PetscObjectReference((PetscObject)perm);
    936:   PetscMalloc((mbs+1)*sizeof(PetscScalar),&b->solve_work);
    937:   PetscLogObjectMemory(B,(ui[mbs]-mbs)*(sizeof(PetscInt)+sizeof(MatScalar)));
    938:   b->maxnz = b->nz = ui[mbs];
    939:
    940:   B->factor                 = FACTOR_CHOLESKY;
    941:   B->info.factor_mallocs    = reallocs;
    942:   B->info.fill_ratio_given  = fill;
    943:   if (ai[mbs] != 0) {
    944:     B->info.fill_ratio_needed = ((PetscReal)ui[mbs])/((PetscReal)ai[mbs]);
    945:   } else {
    946:     B->info.fill_ratio_needed = 0.0;
    947:   }
    948:   if (perm_identity){
    949:     B->ops->solve           = MatSolve_SeqSBAIJ_1_NaturalOrdering;
    950:     B->ops->solvetranspose  = MatSolve_SeqSBAIJ_1_NaturalOrdering;
    951:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering;
    952:   } else {
    953:     B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N;
    954:   }
    955:   return(0);
    956: }