Actual source code: baijfact2.c

      1: /*
      2:     Factorization code for BAIJ format.
      3: */

      5: #include <../src/mat/impls/baij/seq/baij.h>
      6: #include <petsc/private/kernels/blockinvert.h>
      7: #include <petscbt.h>
      8: #include <../src/mat/utils/freespace.h>

     10: extern PetscErrorCode MatDuplicateNoCreate_SeqBAIJ(Mat, Mat, MatDuplicateOption, PetscBool);

     12: /*
     13:    This is not much faster than MatLUFactorNumeric_SeqBAIJ_N() but the solve is faster at least sometimes
     14: */
     15: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_15_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info)
     16: {
     17:   Mat              C = B;
     18:   Mat_SeqBAIJ     *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
     19:   PetscInt         i, j, k, ipvt[15];
     20:   const PetscInt   n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *ajtmp, *bjtmp, *bdiag = b->diag, *pj;
     21:   PetscInt         nz, nzL, row;
     22:   MatScalar       *rtmp, *pc, *mwork, *pv, *vv, work[225];
     23:   const MatScalar *v, *aa = a->a;
     24:   PetscInt         bs2 = a->bs2, bs = A->rmap->bs, flg;
     25:   PetscInt         sol_ver;
     26:   PetscBool        allowzeropivot, zeropivotdetected;

     28:   PetscFunctionBegin;
     29:   allowzeropivot = PetscNot(A->erroriffailure);
     30:   PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)A)->prefix, "-sol_ver", &sol_ver, NULL));

     32:   /* generate work space needed by the factorization */
     33:   PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
     34:   PetscCall(PetscArrayzero(rtmp, bs2 * n));

     36:   for (i = 0; i < n; i++) {
     37:     /* zero rtmp */
     38:     /* L part */
     39:     nz    = bi[i + 1] - bi[i];
     40:     bjtmp = bj + bi[i];
     41:     for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));

     43:     /* U part */
     44:     nz    = bdiag[i] - bdiag[i + 1];
     45:     bjtmp = bj + bdiag[i + 1] + 1;
     46:     for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));

     48:     /* load in initial (unfactored row) */
     49:     nz    = ai[i + 1] - ai[i];
     50:     ajtmp = aj + ai[i];
     51:     v     = aa + bs2 * ai[i];
     52:     for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ajtmp[j], v + bs2 * j, bs2));

     54:     /* elimination */
     55:     bjtmp = bj + bi[i];
     56:     nzL   = bi[i + 1] - bi[i];
     57:     for (k = 0; k < nzL; k++) {
     58:       row = bjtmp[k];
     59:       pc  = rtmp + bs2 * row;
     60:       for (flg = 0, j = 0; j < bs2; j++) {
     61:         if (pc[j] != 0.0) {
     62:           flg = 1;
     63:           break;
     64:         }
     65:       }
     66:       if (flg) {
     67:         pv = b->a + bs2 * bdiag[row];
     68:         PetscKernel_A_gets_A_times_B(bs, pc, pv, mwork);
     69:         /* PetscCall(PetscKernel_A_gets_A_times_B_15(pc,pv,mwork)); */
     70:         pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
     71:         pv = b->a + bs2 * (bdiag[row + 1] + 1);
     72:         nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
     73:         for (j = 0; j < nz; j++) {
     74:           vv = rtmp + bs2 * pj[j];
     75:           PetscKernel_A_gets_A_minus_B_times_C(bs, vv, pc, pv);
     76:           /* PetscCall(PetscKernel_A_gets_A_minus_B_times_C_15(vv,pc,pv)); */
     77:           pv += bs2;
     78:         }
     79:         PetscCall(PetscLogFlops(2.0 * bs2 * bs * (nz + 1) - bs2)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
     80:       }
     81:     }

     83:     /* finished row so stick it into b->a */
     84:     /* L part */
     85:     pv = b->a + bs2 * bi[i];
     86:     pj = b->j + bi[i];
     87:     nz = bi[i + 1] - bi[i];
     88:     for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));

     90:     /* Mark diagonal and invert diagonal for simpler triangular solves */
     91:     pv = b->a + bs2 * bdiag[i];
     92:     pj = b->j + bdiag[i];
     93:     PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
     94:     PetscCall(PetscKernel_A_gets_inverse_A_15(pv, ipvt, work, info->shiftamount, allowzeropivot, &zeropivotdetected));
     95:     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

     97:     /* U part */
     98:     pv = b->a + bs2 * (bdiag[i + 1] + 1);
     99:     pj = b->j + bdiag[i + 1] + 1;
    100:     nz = bdiag[i] - bdiag[i + 1] - 1;
    101:     for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
    102:   }

    104:   PetscCall(PetscFree2(rtmp, mwork));

    106:   C->ops->solve          = MatSolve_SeqBAIJ_15_NaturalOrdering_ver1;
    107:   C->ops->solvetranspose = MatSolve_SeqBAIJ_N_NaturalOrdering;
    108:   C->assembled           = PETSC_TRUE;

    110:   PetscCall(PetscLogFlops(1.333333333333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */
    111:   PetscFunctionReturn(PETSC_SUCCESS);
    112: }

    114: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_N(Mat B, Mat A, const MatFactorInfo *info)
    115: {
    116:   Mat             C = B;
    117:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
    118:   IS              isrow = b->row, isicol = b->icol;
    119:   const PetscInt *r, *ic;
    120:   PetscInt        i, j, k, n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
    121:   PetscInt       *ajtmp, *bjtmp, nz, nzL, row, *bdiag = b->diag, *pj;
    122:   MatScalar      *rtmp, *pc, *mwork, *v, *pv, *aa     = a->a;
    123:   PetscInt        bs = A->rmap->bs, bs2 = a->bs2, *v_pivots, flg;
    124:   MatScalar      *v_work;
    125:   PetscBool       col_identity, row_identity, both_identity;
    126:   PetscBool       allowzeropivot, zeropivotdetected;

    128:   PetscFunctionBegin;
    129:   PetscCall(ISGetIndices(isrow, &r));
    130:   PetscCall(ISGetIndices(isicol, &ic));
    131:   allowzeropivot = PetscNot(A->erroriffailure);

    133:   PetscCall(PetscCalloc1(bs2 * n, &rtmp));

    135:   /* generate work space needed by dense LU factorization */
    136:   PetscCall(PetscMalloc3(bs, &v_work, bs2, &mwork, bs, &v_pivots));

    138:   for (i = 0; i < n; i++) {
    139:     /* zero rtmp */
    140:     /* L part */
    141:     nz    = bi[i + 1] - bi[i];
    142:     bjtmp = bj + bi[i];
    143:     for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));

    145:     /* U part */
    146:     nz    = bdiag[i] - bdiag[i + 1];
    147:     bjtmp = bj + bdiag[i + 1] + 1;
    148:     for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));

    150:     /* load in initial (unfactored row) */
    151:     nz    = ai[r[i] + 1] - ai[r[i]];
    152:     ajtmp = aj + ai[r[i]];
    153:     v     = aa + bs2 * ai[r[i]];
    154:     for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ic[ajtmp[j]], v + bs2 * j, bs2));

    156:     /* elimination */
    157:     bjtmp = bj + bi[i];
    158:     nzL   = bi[i + 1] - bi[i];
    159:     for (k = 0; k < nzL; k++) {
    160:       row = bjtmp[k];
    161:       pc  = rtmp + bs2 * row;
    162:       for (flg = 0, j = 0; j < bs2; j++) {
    163:         if (pc[j] != 0.0) {
    164:           flg = 1;
    165:           break;
    166:         }
    167:       }
    168:       if (flg) {
    169:         pv = b->a + bs2 * bdiag[row];
    170:         PetscKernel_A_gets_A_times_B(bs, pc, pv, mwork); /* *pc = *pc * (*pv); */
    171:         pj = b->j + bdiag[row + 1] + 1;                  /* beginning of U(row,:) */
    172:         pv = b->a + bs2 * (bdiag[row + 1] + 1);
    173:         nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
    174:         for (j = 0; j < nz; j++) PetscKernel_A_gets_A_minus_B_times_C(bs, rtmp + bs2 * pj[j], pc, pv + bs2 * j);
    175:         PetscCall(PetscLogFlops(2.0 * bs2 * bs * (nz + 1) - bs2)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
    176:       }
    177:     }

    179:     /* finished row so stick it into b->a */
    180:     /* L part */
    181:     pv = b->a + bs2 * bi[i];
    182:     pj = b->j + bi[i];
    183:     nz = bi[i + 1] - bi[i];
    184:     for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));

    186:     /* Mark diagonal and invert diagonal for simpler triangular solves */
    187:     pv = b->a + bs2 * bdiag[i];
    188:     pj = b->j + bdiag[i];
    189:     PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));

    191:     PetscCall(PetscKernel_A_gets_inverse_A(bs, pv, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
    192:     if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

    194:     /* U part */
    195:     pv = b->a + bs2 * (bdiag[i + 1] + 1);
    196:     pj = b->j + bdiag[i + 1] + 1;
    197:     nz = bdiag[i] - bdiag[i + 1] - 1;
    198:     for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
    199:   }

    201:   PetscCall(PetscFree(rtmp));
    202:   PetscCall(PetscFree3(v_work, mwork, v_pivots));
    203:   PetscCall(ISRestoreIndices(isicol, &ic));
    204:   PetscCall(ISRestoreIndices(isrow, &r));

    206:   PetscCall(ISIdentity(isrow, &row_identity));
    207:   PetscCall(ISIdentity(isicol, &col_identity));

    209:   both_identity = (PetscBool)(row_identity && col_identity);
    210:   if (both_identity) {
    211:     switch (bs) {
    212:     case 9:
    213: #if defined(PETSC_HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && defined(PETSC_USE_REAL_DOUBLE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_64BIT_INDICES)
    214:       C->ops->solve = MatSolve_SeqBAIJ_9_NaturalOrdering;
    215: #else
    216:       C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
    217: #endif
    218:       break;
    219:     case 11:
    220:       C->ops->solve = MatSolve_SeqBAIJ_11_NaturalOrdering;
    221:       break;
    222:     case 12:
    223:       C->ops->solve = MatSolve_SeqBAIJ_12_NaturalOrdering;
    224:       break;
    225:     case 13:
    226:       C->ops->solve = MatSolve_SeqBAIJ_13_NaturalOrdering;
    227:       break;
    228:     case 14:
    229:       C->ops->solve = MatSolve_SeqBAIJ_14_NaturalOrdering;
    230:       break;
    231:     default:
    232:       C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering;
    233:       break;
    234:     }
    235:   } else {
    236:     C->ops->solve = MatSolve_SeqBAIJ_N;
    237:   }
    238:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N;

    240:   C->assembled = PETSC_TRUE;

    242:   PetscCall(PetscLogFlops(1.333333333333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */
    243:   PetscFunctionReturn(PETSC_SUCCESS);
    244: }

    246: /*
    247:    ilu(0) with natural ordering under new data structure.
    248:    See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed description
    249:    because this code is almost identical to MatILUFactorSymbolic_SeqAIJ_ilu0_inplace().
    250: */

    252: static PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_ilu0(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
    253: {
    254:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b;
    255:   PetscInt     n = a->mbs, *ai = a->i, *aj, *adiag = a->diag, bs2 = a->bs2;
    256:   PetscInt     i, j, nz, *bi, *bj, *bdiag, bi_temp;

    258:   PetscFunctionBegin;
    259:   PetscCall(MatDuplicateNoCreate_SeqBAIJ(fact, A, MAT_DO_NOT_COPY_VALUES, PETSC_FALSE));
    260:   b = (Mat_SeqBAIJ *)fact->data;

    262:   /* allocate matrix arrays for new data structure */
    263:   PetscCall(PetscShmgetAllocateArray(bs2 * ai[n] + 1, sizeof(PetscScalar), (void **)&b->a));
    264:   PetscCall(PetscShmgetAllocateArray(ai[n] + 1, sizeof(PetscInt), (void **)&b->j));
    265:   PetscCall(PetscShmgetAllocateArray(n + 1, sizeof(PetscInt), (void **)&b->i));
    266:   b->free_a          = PETSC_TRUE;
    267:   b->free_ij         = PETSC_TRUE;
    268:   fact->preallocated = PETSC_TRUE;
    269:   fact->assembled    = PETSC_TRUE;
    270:   if (!b->diag) PetscCall(PetscMalloc1(n + 1, &b->diag));
    271:   bdiag = b->diag;

    273:   if (n > 0) PetscCall(PetscArrayzero(b->a, bs2 * ai[n]));

    275:   /* set bi and bj with new data structure */
    276:   bi = b->i;
    277:   bj = b->j;

    279:   /* L part */
    280:   bi[0] = 0;
    281:   for (i = 0; i < n; i++) {
    282:     nz        = adiag[i] - ai[i];
    283:     bi[i + 1] = bi[i] + nz;
    284:     aj        = a->j + ai[i];
    285:     for (j = 0; j < nz; j++) {
    286:       *bj = aj[j];
    287:       bj++;
    288:     }
    289:   }

    291:   /* U part */
    292:   bi_temp  = bi[n];
    293:   bdiag[n] = bi[n] - 1;
    294:   for (i = n - 1; i >= 0; i--) {
    295:     nz      = ai[i + 1] - adiag[i] - 1;
    296:     bi_temp = bi_temp + nz + 1;
    297:     aj      = a->j + adiag[i] + 1;
    298:     for (j = 0; j < nz; j++) {
    299:       *bj = aj[j];
    300:       bj++;
    301:     }
    302:     /* diag[i] */
    303:     *bj = i;
    304:     bj++;
    305:     bdiag[i] = bi_temp - 1;
    306:   }
    307:   PetscFunctionReturn(PETSC_SUCCESS);
    308: }

    310: PetscErrorCode MatILUFactorSymbolic_SeqBAIJ(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
    311: {
    312:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ *)A->data, *b;
    313:   IS                 isicol;
    314:   const PetscInt    *r, *ic;
    315:   PetscInt           n = a->mbs, *ai = a->i, *aj = a->j, d;
    316:   PetscInt          *bi, *cols, nnz, *cols_lvl;
    317:   PetscInt          *bdiag, prow, fm, nzbd, reallocs = 0, dcount = 0;
    318:   PetscInt           i, levels, diagonal_fill;
    319:   PetscBool          col_identity, row_identity, both_identity;
    320:   PetscReal          f;
    321:   PetscInt           nlnk, *lnk, *lnk_lvl = NULL;
    322:   PetscBT            lnkbt;
    323:   PetscInt           nzi, *bj, **bj_ptr, **bjlvl_ptr;
    324:   PetscFreeSpaceList free_space = NULL, current_space = NULL;
    325:   PetscFreeSpaceList free_space_lvl = NULL, current_space_lvl = NULL;
    326:   PetscBool          missing;
    327:   PetscInt           bs = A->rmap->bs, bs2 = a->bs2;

    329:   PetscFunctionBegin;
    330:   PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must be square matrix, rows %" PetscInt_FMT " columns %" PetscInt_FMT, A->rmap->n, A->cmap->n);
    331:   if (bs > 1) { /* check shifttype */
    332:     PetscCheck(info->shifttype != (PetscReal)MAT_SHIFT_NONZERO && info->shifttype != (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only MAT_SHIFT_NONE and MAT_SHIFT_INBLOCKS are supported for BAIJ matrix");
    333:   }

    335:   PetscCall(MatMissingDiagonal(A, &missing, &d));
    336:   PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry %" PetscInt_FMT, d);

    338:   f             = info->fill;
    339:   levels        = (PetscInt)info->levels;
    340:   diagonal_fill = (PetscInt)info->diagonal_fill;

    342:   PetscCall(ISInvertPermutation(iscol, PETSC_DECIDE, &isicol));

    344:   PetscCall(ISIdentity(isrow, &row_identity));
    345:   PetscCall(ISIdentity(iscol, &col_identity));

    347:   both_identity = (PetscBool)(row_identity && col_identity);

    349:   if (!levels && both_identity) {
    350:     /* special case: ilu(0) with natural ordering */
    351:     PetscCall(MatILUFactorSymbolic_SeqBAIJ_ilu0(fact, A, isrow, iscol, info));
    352:     PetscCall(MatSeqBAIJSetNumericFactorization(fact, both_identity));

    354:     fact->factortype             = MAT_FACTOR_ILU;
    355:     fact->info.factor_mallocs    = 0;
    356:     fact->info.fill_ratio_given  = info->fill;
    357:     fact->info.fill_ratio_needed = 1.0;

    359:     b       = (Mat_SeqBAIJ *)fact->data;
    360:     b->row  = isrow;
    361:     b->col  = iscol;
    362:     b->icol = isicol;
    363:     PetscCall(PetscObjectReference((PetscObject)isrow));
    364:     PetscCall(PetscObjectReference((PetscObject)iscol));
    365:     b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;

    367:     PetscCall(PetscMalloc1((n + 1) * bs, &b->solve_work));
    368:     PetscFunctionReturn(PETSC_SUCCESS);
    369:   }

    371:   PetscCall(ISGetIndices(isrow, &r));
    372:   PetscCall(ISGetIndices(isicol, &ic));

    374:   /* get new row pointers */
    375:   PetscCall(PetscMalloc1(n + 1, &bi));
    376:   bi[0] = 0;
    377:   /* bdiag is location of diagonal in factor */
    378:   PetscCall(PetscMalloc1(n + 1, &bdiag));
    379:   bdiag[0] = 0;

    381:   PetscCall(PetscMalloc2(n, &bj_ptr, n, &bjlvl_ptr));

    383:   /* create a linked list for storing column indices of the active row */
    384:   nlnk = n + 1;
    385:   PetscCall(PetscIncompleteLLCreate(n, n, nlnk, lnk, lnk_lvl, lnkbt));

    387:   /* initial FreeSpace size is f*(ai[n]+1) */
    388:   PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space));
    389:   current_space = free_space;
    390:   PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f, ai[n] + 1), &free_space_lvl));
    391:   current_space_lvl = free_space_lvl;

    393:   for (i = 0; i < n; i++) {
    394:     nzi = 0;
    395:     /* copy current row into linked list */
    396:     nnz = ai[r[i] + 1] - ai[r[i]];
    397:     PetscCheck(nnz, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Empty row in matrix: row in original ordering %" PetscInt_FMT " in permuted ordering %" PetscInt_FMT, r[i], i);
    398:     cols   = aj + ai[r[i]];
    399:     lnk[i] = -1; /* marker to indicate if diagonal exists */
    400:     PetscCall(PetscIncompleteLLInit(nnz, cols, n, ic, &nlnk, lnk, lnk_lvl, lnkbt));
    401:     nzi += nlnk;

    403:     /* make sure diagonal entry is included */
    404:     if (diagonal_fill && lnk[i] == -1) {
    405:       fm = n;
    406:       while (lnk[fm] < i) fm = lnk[fm];
    407:       lnk[i]     = lnk[fm]; /* insert diagonal into linked list */
    408:       lnk[fm]    = i;
    409:       lnk_lvl[i] = 0;
    410:       nzi++;
    411:       dcount++;
    412:     }

    414:     /* add pivot rows into the active row */
    415:     nzbd = 0;
    416:     prow = lnk[n];
    417:     while (prow < i) {
    418:       nnz      = bdiag[prow];
    419:       cols     = bj_ptr[prow] + nnz + 1;
    420:       cols_lvl = bjlvl_ptr[prow] + nnz + 1;
    421:       nnz      = bi[prow + 1] - bi[prow] - nnz - 1;

    423:       PetscCall(PetscILULLAddSorted(nnz, cols, levels, cols_lvl, prow, &nlnk, lnk, lnk_lvl, lnkbt, prow));
    424:       nzi += nlnk;
    425:       prow = lnk[prow];
    426:       nzbd++;
    427:     }
    428:     bdiag[i]  = nzbd;
    429:     bi[i + 1] = bi[i] + nzi;

    431:     /* if free space is not available, make more free space */
    432:     if (current_space->local_remaining < nzi) {
    433:       nnz = PetscIntMultTruncate(2, PetscIntMultTruncate(nzi, n - i)); /* estimated and max additional space needed */
    434:       PetscCall(PetscFreeSpaceGet(nnz, &current_space));
    435:       PetscCall(PetscFreeSpaceGet(nnz, &current_space_lvl));
    436:       reallocs++;
    437:     }

    439:     /* copy data into free_space and free_space_lvl, then initialize lnk */
    440:     PetscCall(PetscIncompleteLLClean(n, n, nzi, lnk, lnk_lvl, current_space->array, current_space_lvl->array, lnkbt));

    442:     bj_ptr[i]    = current_space->array;
    443:     bjlvl_ptr[i] = current_space_lvl->array;

    445:     /* make sure the active row i has diagonal entry */
    446:     PetscCheck(*(bj_ptr[i] + bdiag[i]) == i, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Row %" PetscInt_FMT " has missing diagonal in factored matrix, try running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill", i);

    448:     current_space->array += nzi;
    449:     current_space->local_used += nzi;
    450:     current_space->local_remaining -= nzi;

    452:     current_space_lvl->array += nzi;
    453:     current_space_lvl->local_used += nzi;
    454:     current_space_lvl->local_remaining -= nzi;
    455:   }

    457:   PetscCall(ISRestoreIndices(isrow, &r));
    458:   PetscCall(ISRestoreIndices(isicol, &ic));

    460:   /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */
    461:   PetscCall(PetscMalloc1(bi[n] + 1, &bj));
    462:   PetscCall(PetscFreeSpaceContiguous_LU(&free_space, bj, n, bi, bdiag));

    464:   PetscCall(PetscIncompleteLLDestroy(lnk, lnkbt));
    465:   PetscCall(PetscFreeSpaceDestroy(free_space_lvl));
    466:   PetscCall(PetscFree2(bj_ptr, bjlvl_ptr));

    468: #if defined(PETSC_USE_INFO)
    469:   {
    470:     PetscReal af = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);
    471:     PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocs, (double)f, (double)af));
    472:     PetscCall(PetscInfo(A, "Run with -[sub_]pc_factor_fill %g or use \n", (double)af));
    473:     PetscCall(PetscInfo(A, "PCFactorSetFill([sub]pc,%g);\n", (double)af));
    474:     PetscCall(PetscInfo(A, "for best performance.\n"));
    475:     if (diagonal_fill) PetscCall(PetscInfo(A, "Detected and replaced %" PetscInt_FMT " missing diagonals\n", dcount));
    476:   }
    477: #endif

    479:   /* put together the new matrix */
    480:   PetscCall(MatSeqBAIJSetPreallocation(fact, bs, MAT_SKIP_ALLOCATION, NULL));

    482:   b          = (Mat_SeqBAIJ *)fact->data;
    483:   b->free_ij = PETSC_TRUE;
    484:   PetscCall(PetscShmgetAllocateArray(bs2 * (bdiag[0] + 1), sizeof(PetscScalar), (void **)&b->a));
    485:   b->free_a = PETSC_TRUE;

    487:   b->j         = bj;
    488:   b->i         = bi;
    489:   b->diag      = bdiag;
    490:   b->free_diag = PETSC_TRUE;
    491:   b->ilen      = NULL;
    492:   b->imax      = NULL;
    493:   b->row       = isrow;
    494:   b->col       = iscol;
    495:   PetscCall(PetscObjectReference((PetscObject)isrow));
    496:   PetscCall(PetscObjectReference((PetscObject)iscol));
    497:   b->icol = isicol;

    499:   PetscCall(PetscMalloc1(bs * n + bs, &b->solve_work));
    500:   /* In b structure:  Free imax, ilen, old a, old j.
    501:      Allocate bdiag, solve_work, new a, new j */
    502:   b->maxnz = b->nz = bdiag[0] + 1;

    504:   fact->info.factor_mallocs    = reallocs;
    505:   fact->info.fill_ratio_given  = f;
    506:   fact->info.fill_ratio_needed = ((PetscReal)(bdiag[0] + 1)) / ((PetscReal)ai[n]);

    508:   PetscCall(MatSeqBAIJSetNumericFactorization(fact, both_identity));
    509:   PetscFunctionReturn(PETSC_SUCCESS);
    510: }

    512: #if 0
    513: // unused
    514: /*
    515:      This code is virtually identical to MatILUFactorSymbolic_SeqAIJ
    516:    except that the data structure of Mat_SeqAIJ is slightly different.
    517:    Not a good example of code reuse.
    518: */
    519: static PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_inplace(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
    520: {
    521:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ *)A->data, *b;
    522:   IS              isicol;
    523:   const PetscInt *r, *ic, *ai = a->i, *aj = a->j, *xi;
    524:   PetscInt        prow, n = a->mbs, *ainew, *ajnew, jmax, *fill, nz, *im, *ajfill, *flev, *xitmp;
    525:   PetscInt       *dloc, idx, row, m, fm, nzf, nzi, reallocate = 0, dcount = 0;
    526:   PetscInt        incrlev, nnz, i, bs = A->rmap->bs, bs2 = a->bs2, levels, diagonal_fill, dd;
    527:   PetscBool       col_identity, row_identity, both_identity, flg;
    528:   PetscReal       f;

    530:   PetscFunctionBegin;
    531:   PetscCall(MatMissingDiagonal_SeqBAIJ(A, &flg, &dd));
    532:   PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix A is missing diagonal entry in row %" PetscInt_FMT, dd);

    534:   f             = info->fill;
    535:   levels        = (PetscInt)info->levels;
    536:   diagonal_fill = (PetscInt)info->diagonal_fill;

    538:   PetscCall(ISInvertPermutation(iscol, PETSC_DECIDE, &isicol));

    540:   PetscCall(ISIdentity(isrow, &row_identity));
    541:   PetscCall(ISIdentity(iscol, &col_identity));
    542:   both_identity = (PetscBool)(row_identity && col_identity);

    544:   if (!levels && both_identity) { /* special case copy the nonzero structure */
    545:     PetscCall(MatDuplicateNoCreate_SeqBAIJ(fact, A, MAT_DO_NOT_COPY_VALUES, PETSC_TRUE));
    546:     PetscCall(MatSeqBAIJSetNumericFactorization_inplace(fact, both_identity));

    548:     fact->factortype = MAT_FACTOR_ILU;
    549:     b                = (Mat_SeqBAIJ *)fact->data;
    550:     b->row           = isrow;
    551:     b->col           = iscol;
    552:     PetscCall(PetscObjectReference((PetscObject)isrow));
    553:     PetscCall(PetscObjectReference((PetscObject)iscol));
    554:     b->icol          = isicol;
    555:     b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;

    557:     PetscCall(PetscMalloc1((n + 1) * bs, &b->solve_work));
    558:     PetscFunctionReturn(PETSC_SUCCESS);
    559:   }

    561:   /* general case perform the symbolic factorization */
    562:   PetscCall(ISGetIndices(isrow, &r));
    563:   PetscCall(ISGetIndices(isicol, &ic));

    565:   /* get new row pointers */
    566:   PetscCall(PetscMalloc1(n + 1, &ainew));
    567:   ainew[0] = 0;
    568:   /* don't know how many column pointers are needed so estimate */
    569:   jmax = (PetscInt)(f * ai[n] + 1);
    570:   PetscCall(PetscMalloc1(jmax, &ajnew));
    571:   /* ajfill is level of fill for each fill entry */
    572:   PetscCall(PetscMalloc1(jmax, &ajfill));
    573:   /* fill is a linked list of nonzeros in active row */
    574:   PetscCall(PetscMalloc1(n + 1, &fill));
    575:   /* im is level for each filled value */
    576:   PetscCall(PetscMalloc1(n + 1, &im));
    577:   /* dloc is location of diagonal in factor */
    578:   PetscCall(PetscMalloc1(n + 1, &dloc));
    579:   dloc[0] = 0;
    580:   for (prow = 0; prow < n; prow++) {
    581:     /* copy prow into linked list */
    582:     nzf = nz = ai[r[prow] + 1] - ai[r[prow]];
    583:     PetscCheck(nz, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Empty row in matrix: row in original ordering %" PetscInt_FMT " in permuted ordering %" PetscInt_FMT, r[prow], prow);
    584:     xi         = aj + ai[r[prow]];
    585:     fill[n]    = n;
    586:     fill[prow] = -1; /* marker for diagonal entry */
    587:     while (nz--) {
    588:       fm  = n;
    589:       idx = ic[*xi++];
    590:       do {
    591:         m  = fm;
    592:         fm = fill[m];
    593:       } while (fm < idx);
    594:       fill[m]   = idx;
    595:       fill[idx] = fm;
    596:       im[idx]   = 0;
    597:     }

    599:     /* make sure diagonal entry is included */
    600:     if (diagonal_fill && fill[prow] == -1) {
    601:       fm = n;
    602:       while (fill[fm] < prow) fm = fill[fm];
    603:       fill[prow] = fill[fm]; /* insert diagonal into linked list */
    604:       fill[fm]   = prow;
    605:       im[prow]   = 0;
    606:       nzf++;
    607:       dcount++;
    608:     }

    610:     nzi = 0;
    611:     row = fill[n];
    612:     while (row < prow) {
    613:       incrlev = im[row] + 1;
    614:       nz      = dloc[row];
    615:       xi      = ajnew + ainew[row] + nz + 1;
    616:       flev    = ajfill + ainew[row] + nz + 1;
    617:       nnz     = ainew[row + 1] - ainew[row] - nz - 1;
    618:       fm      = row;
    619:       while (nnz-- > 0) {
    620:         idx = *xi++;
    621:         if (*flev + incrlev > levels) {
    622:           flev++;
    623:           continue;
    624:         }
    625:         do {
    626:           m  = fm;
    627:           fm = fill[m];
    628:         } while (fm < idx);
    629:         if (fm != idx) {
    630:           im[idx]   = *flev + incrlev;
    631:           fill[m]   = idx;
    632:           fill[idx] = fm;
    633:           fm        = idx;
    634:           nzf++;
    635:         } else if (im[idx] > *flev + incrlev) im[idx] = *flev + incrlev;
    636:         flev++;
    637:       }
    638:       row = fill[row];
    639:       nzi++;
    640:     }
    641:     /* copy new filled row into permanent storage */
    642:     ainew[prow + 1] = ainew[prow] + nzf;
    643:     if (ainew[prow + 1] > jmax) {
    644:       /* estimate how much additional space we will need */
    645:       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
    646:       /* just double the memory each time */
    647:       PetscInt maxadd = jmax;
    648:       /* maxadd = (int)(((f*ai[n]+1)*(n-prow+5))/n); */
    649:       if (maxadd < nzf) maxadd = (n - prow) * (nzf + 1);
    650:       jmax += maxadd;

    652:       /* allocate a longer ajnew and ajfill */
    653:       PetscCall(PetscMalloc1(jmax, &xitmp));
    654:       PetscCall(PetscArraycpy(xitmp, ajnew, ainew[prow]));
    655:       PetscCall(PetscFree(ajnew));
    656:       ajnew = xitmp;
    657:       PetscCall(PetscMalloc1(jmax, &xitmp));
    658:       PetscCall(PetscArraycpy(xitmp, ajfill, ainew[prow]));
    659:       PetscCall(PetscFree(ajfill));
    660:       ajfill = xitmp;
    661:       reallocate++; /* count how many reallocations are needed */
    662:     }
    663:     xitmp      = ajnew + ainew[prow];
    664:     flev       = ajfill + ainew[prow];
    665:     dloc[prow] = nzi;
    666:     fm         = fill[n];
    667:     while (nzf--) {
    668:       *xitmp++ = fm;
    669:       *flev++  = im[fm];
    670:       fm       = fill[fm];
    671:     }
    672:     /* make sure row has diagonal entry */
    673:     PetscCheck(ajnew[ainew[prow] + dloc[prow]] == prow, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Row %" PetscInt_FMT " has missing diagonal in factored matrix, try running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill", prow);
    674:   }
    675:   PetscCall(PetscFree(ajfill));
    676:   PetscCall(ISRestoreIndices(isrow, &r));
    677:   PetscCall(ISRestoreIndices(isicol, &ic));
    678:   PetscCall(PetscFree(fill));
    679:   PetscCall(PetscFree(im));

    681:   #if defined(PETSC_USE_INFO)
    682:   {
    683:     PetscReal af = ((PetscReal)ainew[n]) / ((PetscReal)ai[n]);
    684:     PetscCall(PetscInfo(A, "Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n", reallocate, (double)f, (double)af));
    685:     PetscCall(PetscInfo(A, "Run with -pc_factor_fill %g or use \n", (double)af));
    686:     PetscCall(PetscInfo(A, "PCFactorSetFill(pc,%g);\n", (double)af));
    687:     PetscCall(PetscInfo(A, "for best performance.\n"));
    688:     if (diagonal_fill) PetscCall(PetscInfo(A, "Detected and replaced %" PetscInt_FMT " missing diagonals\n", dcount));
    689:   }
    690:   #endif

    692:   /* put together the new matrix */
    693:   PetscCall(MatSeqBAIJSetPreallocation(fact, bs, MAT_SKIP_ALLOCATION, NULL));
    694:   b = (Mat_SeqBAIJ *)fact->data;
    695:   b->free_ij      = PETSC_TRUE;
    696:   PetscCall(PetscShmgetAllocateArray(bs2 * ainew[n],sizeof(PetscScalar),(void **)&b->a));
    697:   b->free_a       = PETSC_TRUE;
    698:   b->j = ajnew;
    699:   b->i = ainew;
    700:   for (i = 0; i < n; i++) dloc[i] += ainew[i];
    701:   b->diag          = dloc;
    702:   b->free_diag     = PETSC_TRUE;
    703:   b->ilen          = NULL;
    704:   b->imax          = NULL;
    705:   b->row           = isrow;
    706:   b->col           = iscol;
    707:   b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;

    709:   PetscCall(PetscObjectReference((PetscObject)isrow));
    710:   PetscCall(PetscObjectReference((PetscObject)iscol));
    711:   b->icol = isicol;
    712:   PetscCall(PetscMalloc1(bs * n + bs, &b->solve_work));
    713:   /* In b structure:  Free imax, ilen, old a, old j.
    714:      Allocate dloc, solve_work, new a, new j */
    715:   b->maxnz = b->nz = ainew[n];

    717:   fact->info.factor_mallocs    = reallocate;
    718:   fact->info.fill_ratio_given  = f;
    719:   fact->info.fill_ratio_needed = ((PetscReal)ainew[n]) / ((PetscReal)ai[prow]);

    721:   PetscCall(MatSeqBAIJSetNumericFactorization_inplace(fact, both_identity));
    722:   PetscFunctionReturn(PETSC_SUCCESS);
    723: }
    724: #endif

    726: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE(Mat A)
    727: {
    728:   /* Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; */
    729:   /* int i,*AJ=a->j,nz=a->nz; */

    731:   PetscFunctionBegin;
    732:   /* Undo Column scaling */
    733:   /*    while (nz--) { */
    734:   /*      AJ[i] = AJ[i]/4; */
    735:   /*    } */
    736:   /* This should really invoke a push/pop logic, but we don't have that yet. */
    737:   A->ops->setunfactored = NULL;
    738:   PetscFunctionReturn(PETSC_SUCCESS);
    739: }

    741: PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE_usj(Mat A)
    742: {
    743:   Mat_SeqBAIJ    *a  = (Mat_SeqBAIJ *)A->data;
    744:   PetscInt       *AJ = a->j, nz = a->nz;
    745:   unsigned short *aj = (unsigned short *)AJ;

    747:   PetscFunctionBegin;
    748:   /* Is this really necessary? */
    749:   while (nz--) { AJ[nz] = (int)((unsigned int)aj[nz]); /* First extend, then convert to signed. */ }
    750:   A->ops->setunfactored = NULL;
    751:   PetscFunctionReturn(PETSC_SUCCESS);
    752: }