Actual source code: matmatmult.c

      1: /*
      2:   Defines matrix-matrix product routines for pairs of SeqAIJ matrices
      3:           C = A * B
      4: */

      6: #include <../src/mat/impls/aij/seq/aij.h>
      7: #include <../src/mat/utils/freespace.h>
      8: #include <petscbt.h>
      9: #include <petsc/private/isimpl.h>
     10: #include <../src/mat/impls/dense/seq/dense.h>

     12: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C)
     13: {
     14:   PetscFunctionBegin;
     15:   if (C->ops->matmultnumeric) PetscCall((*C->ops->matmultnumeric)(A, B, C));
     16:   else PetscCall(MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(A, B, C));
     17:   PetscFunctionReturn(PETSC_SUCCESS);
     18: }

     20: /* Modified from MatCreateSeqAIJWithArrays() */
     21: PETSC_INTERN PetscErrorCode MatSetSeqAIJWithArrays_private(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], MatType mtype, Mat mat)
     22: {
     23:   PetscInt    ii;
     24:   Mat_SeqAIJ *aij;
     25:   PetscBool   isseqaij, ofree_a, ofree_ij;

     27:   PetscFunctionBegin;
     28:   PetscCheck(m <= 0 || !i[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0");
     29:   PetscCall(MatSetSizes(mat, m, n, m, n));

     31:   if (!mtype) {
     32:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQAIJ, &isseqaij));
     33:     if (!isseqaij) PetscCall(MatSetType(mat, MATSEQAIJ));
     34:   } else {
     35:     PetscCall(MatSetType(mat, mtype));
     36:   }

     38:   aij      = (Mat_SeqAIJ *)mat->data;
     39:   ofree_a  = aij->free_a;
     40:   ofree_ij = aij->free_ij;
     41:   /* changes the free flags */
     42:   PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat, MAT_SKIP_ALLOCATION, NULL));

     44:   PetscCall(PetscFree(aij->ilen));
     45:   PetscCall(PetscFree(aij->imax));
     46:   PetscCall(PetscMalloc1(m, &aij->imax));
     47:   PetscCall(PetscMalloc1(m, &aij->ilen));
     48:   for (ii = 0, aij->nonzerorowcnt = 0, aij->rmax = 0; ii < m; ii++) {
     49:     const PetscInt rnz = i[ii + 1] - i[ii];
     50:     aij->nonzerorowcnt += !!rnz;
     51:     aij->rmax     = PetscMax(aij->rmax, rnz);
     52:     aij->ilen[ii] = aij->imax[ii] = i[ii + 1] - i[ii];
     53:   }
     54:   aij->maxnz = i[m];
     55:   aij->nz    = i[m];

     57:   if (ofree_a) PetscCall(PetscShmgetDeallocateArray((void **)&aij->a));
     58:   if (ofree_ij) PetscCall(PetscShmgetDeallocateArray((void **)&aij->j));
     59:   if (ofree_ij) PetscCall(PetscShmgetDeallocateArray((void **)&aij->i));

     61:   aij->i       = i;
     62:   aij->j       = j;
     63:   aij->a       = a;
     64:   aij->nonew   = -1; /* this indicates that inserting a new value in the matrix that generates a new nonzero is an error */
     65:   aij->free_a  = PETSC_FALSE;
     66:   aij->free_ij = PETSC_FALSE;
     67:   PetscCall(MatCheckCompressedRow(mat, aij->nonzerorowcnt, &aij->compressedrow, aij->i, m, 0.6));
     68:   // Always build the diag info when i, j are set
     69:   PetscCall(MatMarkDiagonal_SeqAIJ(mat));
     70:   PetscFunctionReturn(PETSC_SUCCESS);
     71: }

     73: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C)
     74: {
     75:   Mat_Product        *product = C->product;
     76:   MatProductAlgorithm alg;
     77:   PetscBool           flg;

     79:   PetscFunctionBegin;
     80:   if (product) {
     81:     alg = product->alg;
     82:   } else {
     83:     alg = "sorted";
     84:   }
     85:   /* sorted */
     86:   PetscCall(PetscStrcmp(alg, "sorted", &flg));
     87:   if (flg) {
     88:     PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(A, B, fill, C));
     89:     PetscFunctionReturn(PETSC_SUCCESS);
     90:   }

     92:   /* scalable */
     93:   PetscCall(PetscStrcmp(alg, "scalable", &flg));
     94:   if (flg) {
     95:     PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A, B, fill, C));
     96:     PetscFunctionReturn(PETSC_SUCCESS);
     97:   }

     99:   /* scalable_fast */
    100:   PetscCall(PetscStrcmp(alg, "scalable_fast", &flg));
    101:   if (flg) {
    102:     PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(A, B, fill, C));
    103:     PetscFunctionReturn(PETSC_SUCCESS);
    104:   }

    106:   /* heap */
    107:   PetscCall(PetscStrcmp(alg, "heap", &flg));
    108:   if (flg) {
    109:     PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(A, B, fill, C));
    110:     PetscFunctionReturn(PETSC_SUCCESS);
    111:   }

    113:   /* btheap */
    114:   PetscCall(PetscStrcmp(alg, "btheap", &flg));
    115:   if (flg) {
    116:     PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(A, B, fill, C));
    117:     PetscFunctionReturn(PETSC_SUCCESS);
    118:   }

    120:   /* llcondensed */
    121:   PetscCall(PetscStrcmp(alg, "llcondensed", &flg));
    122:   if (flg) {
    123:     PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(A, B, fill, C));
    124:     PetscFunctionReturn(PETSC_SUCCESS);
    125:   }

    127:   /* rowmerge */
    128:   PetscCall(PetscStrcmp(alg, "rowmerge", &flg));
    129:   if (flg) {
    130:     PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(A, B, fill, C));
    131:     PetscFunctionReturn(PETSC_SUCCESS);
    132:   }

    134: #if defined(PETSC_HAVE_HYPRE)
    135:   PetscCall(PetscStrcmp(alg, "hypre", &flg));
    136:   if (flg) {
    137:     PetscCall(MatMatMultSymbolic_AIJ_AIJ_wHYPRE(A, B, fill, C));
    138:     PetscFunctionReturn(PETSC_SUCCESS);
    139:   }
    140: #endif

    142:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Mat Product Algorithm is not supported");
    143: }

    145: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(Mat A, Mat B, PetscReal fill, Mat C)
    146: {
    147:   Mat_SeqAIJ        *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
    148:   PetscInt          *ai = a->i, *bi = b->i, *ci, *cj;
    149:   PetscInt           am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
    150:   PetscReal          afill;
    151:   PetscInt           i, j, anzi, brow, bnzj, cnzi, *bj, *aj, *lnk, ndouble = 0, Crmax;
    152:   PetscHMapI         ta;
    153:   PetscBT            lnkbt;
    154:   PetscFreeSpaceList free_space = NULL, current_space = NULL;

    156:   PetscFunctionBegin;
    157:   /* Get ci and cj */
    158:   /* Allocate ci array, arrays for fill computation and */
    159:   /* free space for accumulating nonzero column info */
    160:   PetscCall(PetscMalloc1(am + 2, &ci));
    161:   ci[0] = 0;

    163:   /* create and initialize a linked list */
    164:   PetscCall(PetscHMapICreateWithSize(bn, &ta));
    165:   MatRowMergeMax_SeqAIJ(b, bm, ta);
    166:   PetscCall(PetscHMapIGetSize(ta, &Crmax));
    167:   PetscCall(PetscHMapIDestroy(&ta));

    169:   PetscCall(PetscLLCondensedCreate(Crmax, bn, &lnk, &lnkbt));

    171:   /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
    172:   PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space));

    174:   current_space = free_space;

    176:   /* Determine ci and cj */
    177:   for (i = 0; i < am; i++) {
    178:     anzi = ai[i + 1] - ai[i];
    179:     aj   = a->j + ai[i];
    180:     for (j = 0; j < anzi; j++) {
    181:       brow = aj[j];
    182:       bnzj = bi[brow + 1] - bi[brow];
    183:       bj   = b->j + bi[brow];
    184:       /* add non-zero cols of B into the sorted linked list lnk */
    185:       PetscCall(PetscLLCondensedAddSorted(bnzj, bj, lnk, lnkbt));
    186:     }
    187:     /* add possible missing diagonal entry */
    188:     if (C->force_diagonals) PetscCall(PetscLLCondensedAddSorted(1, &i, lnk, lnkbt));
    189:     cnzi = lnk[0];

    191:     /* If free space is not available, make more free space */
    192:     /* Double the amount of total space in the list */
    193:     if (current_space->local_remaining < cnzi) {
    194:       PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(cnzi, current_space->total_array_size), &current_space));
    195:       ndouble++;
    196:     }

    198:     /* Copy data into free space, then initialize lnk */
    199:     PetscCall(PetscLLCondensedClean(bn, cnzi, current_space->array, lnk, lnkbt));

    201:     current_space->array += cnzi;
    202:     current_space->local_used += cnzi;
    203:     current_space->local_remaining -= cnzi;

    205:     ci[i + 1] = ci[i] + cnzi;
    206:   }

    208:   /* Column indices are in the list of free space */
    209:   /* Allocate space for cj, initialize cj, and */
    210:   /* destroy list of free space and other temporary array(s) */
    211:   PetscCall(PetscMalloc1(ci[am] + 1, &cj));
    212:   PetscCall(PetscFreeSpaceContiguous(&free_space, cj));
    213:   PetscCall(PetscLLCondensedDestroy(lnk, lnkbt));

    215:   /* put together the new symbolic matrix */
    216:   PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C));
    217:   PetscCall(MatSetBlockSizesFromMats(C, A, B));

    219:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    220:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    221:   c          = (Mat_SeqAIJ *)C->data;
    222:   c->free_a  = PETSC_FALSE;
    223:   c->free_ij = PETSC_TRUE;
    224:   c->nonew   = 0;

    226:   /* fast, needs non-scalable O(bn) array 'abdense' */
    227:   C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

    229:   /* set MatInfo */
    230:   afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
    231:   if (afill < 1.0) afill = 1.0;
    232:   C->info.mallocs           = ndouble;
    233:   C->info.fill_ratio_given  = fill;
    234:   C->info.fill_ratio_needed = afill;

    236: #if defined(PETSC_USE_INFO)
    237:   if (ci[am]) {
    238:     PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
    239:     PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
    240:   } else {
    241:     PetscCall(PetscInfo(C, "Empty matrix product\n"));
    242:   }
    243: #endif
    244:   PetscFunctionReturn(PETSC_SUCCESS);
    245: }

    247: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(Mat A, Mat B, Mat C)
    248: {
    249:   PetscLogDouble     flops = 0.0;
    250:   Mat_SeqAIJ        *a     = (Mat_SeqAIJ *)A->data;
    251:   Mat_SeqAIJ        *b     = (Mat_SeqAIJ *)B->data;
    252:   Mat_SeqAIJ        *c     = (Mat_SeqAIJ *)C->data;
    253:   PetscInt          *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *bjj, *ci = c->i, *cj = c->j;
    254:   PetscInt           am = A->rmap->n, cm = C->rmap->n;
    255:   PetscInt           i, j, k, anzi, bnzi, cnzi, brow;
    256:   PetscScalar       *ca, valtmp;
    257:   PetscScalar       *ab_dense;
    258:   PetscContainer     cab_dense;
    259:   const PetscScalar *aa, *ba, *baj;

    261:   PetscFunctionBegin;
    262:   PetscCall(MatSeqAIJGetArrayRead(A, &aa));
    263:   PetscCall(MatSeqAIJGetArrayRead(B, &ba));
    264:   if (!c->a) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
    265:     PetscCall(PetscMalloc1(ci[cm] + 1, &ca));
    266:     c->a      = ca;
    267:     c->free_a = PETSC_TRUE;
    268:   } else ca = c->a;

    270:   /* TODO this should be done in the symbolic phase */
    271:   /* However, this function is so heavily used (sometimes in an hidden way through multnumeric function pointers
    272:      that is hard to eradicate) */
    273:   PetscCall(PetscObjectQuery((PetscObject)C, "__PETSc__ab_dense", (PetscObject *)&cab_dense));
    274:   if (!cab_dense) {
    275:     PetscCall(PetscMalloc1(B->cmap->N, &ab_dense));
    276:     PetscCall(PetscObjectContainerCompose((PetscObject)C, "__PETSc__ab_dense", ab_dense, PetscCtxDestroyDefault));
    277:   } else PetscCall(PetscContainerGetPointer(cab_dense, (void **)&ab_dense));
    278:   PetscCall(PetscArrayzero(ab_dense, B->cmap->N));

    280:   /* clean old values in C */
    281:   PetscCall(PetscArrayzero(ca, ci[cm]));
    282:   /* Traverse A row-wise. */
    283:   /* Build the ith row in C by summing over nonzero columns in A, */
    284:   /* the rows of B corresponding to nonzeros of A. */
    285:   for (i = 0; i < am; i++) {
    286:     anzi = ai[i + 1] - ai[i];
    287:     for (j = 0; j < anzi; j++) {
    288:       brow = aj[j];
    289:       bnzi = bi[brow + 1] - bi[brow];
    290:       bjj  = PetscSafePointerPlusOffset(bj, bi[brow]);
    291:       baj  = PetscSafePointerPlusOffset(ba, bi[brow]);
    292:       /* perform dense axpy */
    293:       valtmp = aa[j];
    294:       for (k = 0; k < bnzi; k++) ab_dense[bjj[k]] += valtmp * baj[k];
    295:       flops += 2 * bnzi;
    296:     }
    297:     aj = PetscSafePointerPlusOffset(aj, anzi);
    298:     aa = PetscSafePointerPlusOffset(aa, anzi);

    300:     cnzi = ci[i + 1] - ci[i];
    301:     for (k = 0; k < cnzi; k++) {
    302:       ca[k] += ab_dense[cj[k]];
    303:       ab_dense[cj[k]] = 0.0; /* zero ab_dense */
    304:     }
    305:     flops += cnzi;
    306:     cj = PetscSafePointerPlusOffset(cj, cnzi);
    307:     ca += cnzi;
    308:   }
    309: #if defined(PETSC_HAVE_DEVICE)
    310:   if (C->offloadmask != PETSC_OFFLOAD_UNALLOCATED) C->offloadmask = PETSC_OFFLOAD_CPU;
    311: #endif
    312:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
    313:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
    314:   PetscCall(PetscLogFlops(flops));
    315:   PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
    316:   PetscCall(MatSeqAIJRestoreArrayRead(B, &ba));
    317:   PetscFunctionReturn(PETSC_SUCCESS);
    318: }

    320: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(Mat A, Mat B, Mat C)
    321: {
    322:   PetscLogDouble     flops = 0.0;
    323:   Mat_SeqAIJ        *a     = (Mat_SeqAIJ *)A->data;
    324:   Mat_SeqAIJ        *b     = (Mat_SeqAIJ *)B->data;
    325:   Mat_SeqAIJ        *c     = (Mat_SeqAIJ *)C->data;
    326:   PetscInt          *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *bjj, *ci = c->i, *cj = c->j;
    327:   PetscInt           am = A->rmap->N, cm = C->rmap->N;
    328:   PetscInt           i, j, k, anzi, bnzi, cnzi, brow;
    329:   PetscScalar       *ca = c->a, valtmp;
    330:   const PetscScalar *aa, *ba, *baj;
    331:   PetscInt           nextb;

    333:   PetscFunctionBegin;
    334:   PetscCall(MatSeqAIJGetArrayRead(A, &aa));
    335:   PetscCall(MatSeqAIJGetArrayRead(B, &ba));
    336:   if (!ca) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
    337:     PetscCall(PetscMalloc1(ci[cm] + 1, &ca));
    338:     c->a      = ca;
    339:     c->free_a = PETSC_TRUE;
    340:   }

    342:   /* clean old values in C */
    343:   PetscCall(PetscArrayzero(ca, ci[cm]));
    344:   /* Traverse A row-wise. */
    345:   /* Build the ith row in C by summing over nonzero columns in A, */
    346:   /* the rows of B corresponding to nonzeros of A. */
    347:   for (i = 0; i < am; i++) {
    348:     anzi = ai[i + 1] - ai[i];
    349:     cnzi = ci[i + 1] - ci[i];
    350:     for (j = 0; j < anzi; j++) {
    351:       brow = aj[j];
    352:       bnzi = bi[brow + 1] - bi[brow];
    353:       bjj  = bj + bi[brow];
    354:       baj  = ba + bi[brow];
    355:       /* perform sparse axpy */
    356:       valtmp = aa[j];
    357:       nextb  = 0;
    358:       for (k = 0; nextb < bnzi; k++) {
    359:         if (cj[k] == bjj[nextb]) { /* ccol == bcol */
    360:           ca[k] += valtmp * baj[nextb++];
    361:         }
    362:       }
    363:       flops += 2 * bnzi;
    364:     }
    365:     aj += anzi;
    366:     aa += anzi;
    367:     cj += cnzi;
    368:     ca += cnzi;
    369:   }
    370: #if defined(PETSC_HAVE_DEVICE)
    371:   if (C->offloadmask != PETSC_OFFLOAD_UNALLOCATED) C->offloadmask = PETSC_OFFLOAD_CPU;
    372: #endif
    373:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
    374:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
    375:   PetscCall(PetscLogFlops(flops));
    376:   PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
    377:   PetscCall(MatSeqAIJRestoreArrayRead(B, &ba));
    378:   PetscFunctionReturn(PETSC_SUCCESS);
    379: }

    381: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(Mat A, Mat B, PetscReal fill, Mat C)
    382: {
    383:   Mat_SeqAIJ        *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
    384:   PetscInt          *ai = a->i, *bi = b->i, *ci, *cj;
    385:   PetscInt           am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
    386:   MatScalar         *ca;
    387:   PetscReal          afill;
    388:   PetscInt           i, j, anzi, brow, bnzj, cnzi, *bj, *aj, *lnk, ndouble = 0, Crmax;
    389:   PetscHMapI         ta;
    390:   PetscFreeSpaceList free_space = NULL, current_space = NULL;

    392:   PetscFunctionBegin;
    393:   /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_fast() */
    394:   /* Allocate arrays for fill computation and free space for accumulating nonzero column */
    395:   PetscCall(PetscMalloc1(am + 2, &ci));
    396:   ci[0] = 0;

    398:   /* create and initialize a linked list */
    399:   PetscCall(PetscHMapICreateWithSize(bn, &ta));
    400:   MatRowMergeMax_SeqAIJ(b, bm, ta);
    401:   PetscCall(PetscHMapIGetSize(ta, &Crmax));
    402:   PetscCall(PetscHMapIDestroy(&ta));

    404:   PetscCall(PetscLLCondensedCreate_fast(Crmax, &lnk));

    406:   /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
    407:   PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space));
    408:   current_space = free_space;

    410:   /* Determine ci and cj */
    411:   for (i = 0; i < am; i++) {
    412:     anzi = ai[i + 1] - ai[i];
    413:     aj   = a->j + ai[i];
    414:     for (j = 0; j < anzi; j++) {
    415:       brow = aj[j];
    416:       bnzj = bi[brow + 1] - bi[brow];
    417:       bj   = b->j + bi[brow];
    418:       /* add non-zero cols of B into the sorted linked list lnk */
    419:       PetscCall(PetscLLCondensedAddSorted_fast(bnzj, bj, lnk));
    420:     }
    421:     /* add possible missing diagonal entry */
    422:     if (C->force_diagonals) PetscCall(PetscLLCondensedAddSorted_fast(1, &i, lnk));
    423:     cnzi = lnk[1];

    425:     /* If free space is not available, make more free space */
    426:     /* Double the amount of total space in the list */
    427:     if (current_space->local_remaining < cnzi) {
    428:       PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(cnzi, current_space->total_array_size), &current_space));
    429:       ndouble++;
    430:     }

    432:     /* Copy data into free space, then initialize lnk */
    433:     PetscCall(PetscLLCondensedClean_fast(cnzi, current_space->array, lnk));

    435:     current_space->array += cnzi;
    436:     current_space->local_used += cnzi;
    437:     current_space->local_remaining -= cnzi;

    439:     ci[i + 1] = ci[i] + cnzi;
    440:   }

    442:   /* Column indices are in the list of free space */
    443:   /* Allocate space for cj, initialize cj, and */
    444:   /* destroy list of free space and other temporary array(s) */
    445:   PetscCall(PetscMalloc1(ci[am] + 1, &cj));
    446:   PetscCall(PetscFreeSpaceContiguous(&free_space, cj));
    447:   PetscCall(PetscLLCondensedDestroy_fast(lnk));

    449:   /* Allocate space for ca */
    450:   PetscCall(PetscCalloc1(ci[am] + 1, &ca));

    452:   /* put together the new symbolic matrix */
    453:   PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, ca, ((PetscObject)A)->type_name, C));
    454:   PetscCall(MatSetBlockSizesFromMats(C, A, B));

    456:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    457:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    458:   c          = (Mat_SeqAIJ *)C->data;
    459:   c->free_a  = PETSC_TRUE;
    460:   c->free_ij = PETSC_TRUE;
    461:   c->nonew   = 0;

    463:   /* slower, less memory */
    464:   C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable;

    466:   /* set MatInfo */
    467:   afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
    468:   if (afill < 1.0) afill = 1.0;
    469:   C->info.mallocs           = ndouble;
    470:   C->info.fill_ratio_given  = fill;
    471:   C->info.fill_ratio_needed = afill;

    473: #if defined(PETSC_USE_INFO)
    474:   if (ci[am]) {
    475:     PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
    476:     PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
    477:   } else {
    478:     PetscCall(PetscInfo(C, "Empty matrix product\n"));
    479:   }
    480: #endif
    481:   PetscFunctionReturn(PETSC_SUCCESS);
    482: }

    484: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(Mat A, Mat B, PetscReal fill, Mat C)
    485: {
    486:   Mat_SeqAIJ        *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
    487:   PetscInt          *ai = a->i, *bi = b->i, *ci, *cj;
    488:   PetscInt           am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
    489:   MatScalar         *ca;
    490:   PetscReal          afill;
    491:   PetscInt           i, j, anzi, brow, bnzj, cnzi, *bj, *aj, *lnk, ndouble = 0, Crmax;
    492:   PetscHMapI         ta;
    493:   PetscFreeSpaceList free_space = NULL, current_space = NULL;

    495:   PetscFunctionBegin;
    496:   /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_Scalalbe() */
    497:   /* Allocate arrays for fill computation and free space for accumulating nonzero column */
    498:   PetscCall(PetscMalloc1(am + 2, &ci));
    499:   ci[0] = 0;

    501:   /* create and initialize a linked list */
    502:   PetscCall(PetscHMapICreateWithSize(bn, &ta));
    503:   MatRowMergeMax_SeqAIJ(b, bm, ta);
    504:   PetscCall(PetscHMapIGetSize(ta, &Crmax));
    505:   PetscCall(PetscHMapIDestroy(&ta));
    506:   PetscCall(PetscLLCondensedCreate_Scalable(Crmax, &lnk));

    508:   /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
    509:   PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space));
    510:   current_space = free_space;

    512:   /* Determine ci and cj */
    513:   for (i = 0; i < am; i++) {
    514:     anzi = ai[i + 1] - ai[i];
    515:     aj   = a->j + ai[i];
    516:     for (j = 0; j < anzi; j++) {
    517:       brow = aj[j];
    518:       bnzj = bi[brow + 1] - bi[brow];
    519:       bj   = b->j + bi[brow];
    520:       /* add non-zero cols of B into the sorted linked list lnk */
    521:       PetscCall(PetscLLCondensedAddSorted_Scalable(bnzj, bj, lnk));
    522:     }
    523:     /* add possible missing diagonal entry */
    524:     if (C->force_diagonals) PetscCall(PetscLLCondensedAddSorted_Scalable(1, &i, lnk));

    526:     cnzi = lnk[0];

    528:     /* If free space is not available, make more free space */
    529:     /* Double the amount of total space in the list */
    530:     if (current_space->local_remaining < cnzi) {
    531:       PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(cnzi, current_space->total_array_size), &current_space));
    532:       ndouble++;
    533:     }

    535:     /* Copy data into free space, then initialize lnk */
    536:     PetscCall(PetscLLCondensedClean_Scalable(cnzi, current_space->array, lnk));

    538:     current_space->array += cnzi;
    539:     current_space->local_used += cnzi;
    540:     current_space->local_remaining -= cnzi;

    542:     ci[i + 1] = ci[i] + cnzi;
    543:   }

    545:   /* Column indices are in the list of free space */
    546:   /* Allocate space for cj, initialize cj, and */
    547:   /* destroy list of free space and other temporary array(s) */
    548:   PetscCall(PetscMalloc1(ci[am] + 1, &cj));
    549:   PetscCall(PetscFreeSpaceContiguous(&free_space, cj));
    550:   PetscCall(PetscLLCondensedDestroy_Scalable(lnk));

    552:   /* Allocate space for ca */
    553:   PetscCall(PetscCalloc1(ci[am] + 1, &ca));

    555:   /* put together the new symbolic matrix */
    556:   PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, ca, ((PetscObject)A)->type_name, C));
    557:   PetscCall(MatSetBlockSizesFromMats(C, A, B));

    559:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    560:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    561:   c          = (Mat_SeqAIJ *)C->data;
    562:   c->free_a  = PETSC_TRUE;
    563:   c->free_ij = PETSC_TRUE;
    564:   c->nonew   = 0;

    566:   /* slower, less memory */
    567:   C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable;

    569:   /* set MatInfo */
    570:   afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
    571:   if (afill < 1.0) afill = 1.0;
    572:   C->info.mallocs           = ndouble;
    573:   C->info.fill_ratio_given  = fill;
    574:   C->info.fill_ratio_needed = afill;

    576: #if defined(PETSC_USE_INFO)
    577:   if (ci[am]) {
    578:     PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
    579:     PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
    580:   } else {
    581:     PetscCall(PetscInfo(C, "Empty matrix product\n"));
    582:   }
    583: #endif
    584:   PetscFunctionReturn(PETSC_SUCCESS);
    585: }

    587: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(Mat A, Mat B, PetscReal fill, Mat C)
    588: {
    589:   Mat_SeqAIJ        *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
    590:   const PetscInt    *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j;
    591:   PetscInt          *ci, *cj, *bb;
    592:   PetscInt           am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
    593:   PetscReal          afill;
    594:   PetscInt           i, j, col, ndouble = 0;
    595:   PetscFreeSpaceList free_space = NULL, current_space = NULL;
    596:   PetscHeap          h;

    598:   PetscFunctionBegin;
    599:   /* Get ci and cj - by merging sorted rows using a heap */
    600:   /* Allocate arrays for fill computation and free space for accumulating nonzero column */
    601:   PetscCall(PetscMalloc1(am + 2, &ci));
    602:   ci[0] = 0;

    604:   /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
    605:   PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space));
    606:   current_space = free_space;

    608:   PetscCall(PetscHeapCreate(a->rmax, &h));
    609:   PetscCall(PetscMalloc1(a->rmax, &bb));

    611:   /* Determine ci and cj */
    612:   for (i = 0; i < am; i++) {
    613:     const PetscInt  anzi = ai[i + 1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
    614:     const PetscInt *acol = aj + ai[i];        /* column indices of nonzero entries in this row */
    615:     ci[i + 1]            = ci[i];
    616:     /* Populate the min heap */
    617:     for (j = 0; j < anzi; j++) {
    618:       bb[j] = bi[acol[j]];           /* bb points at the start of the row */
    619:       if (bb[j] < bi[acol[j] + 1]) { /* Add if row is nonempty */
    620:         PetscCall(PetscHeapAdd(h, j, bj[bb[j]++]));
    621:       }
    622:     }
    623:     /* Pick off the min element, adding it to free space */
    624:     PetscCall(PetscHeapPop(h, &j, &col));
    625:     while (j >= 0) {
    626:       if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
    627:         PetscCall(PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2, current_space->total_array_size), 16 << 20), &current_space));
    628:         ndouble++;
    629:       }
    630:       *(current_space->array++) = col;
    631:       current_space->local_used++;
    632:       current_space->local_remaining--;
    633:       ci[i + 1]++;

    635:       /* stash if anything else remains in this row of B */
    636:       if (bb[j] < bi[acol[j] + 1]) PetscCall(PetscHeapStash(h, j, bj[bb[j]++]));
    637:       while (1) { /* pop and stash any other rows of B that also had an entry in this column */
    638:         PetscInt j2, col2;
    639:         PetscCall(PetscHeapPeek(h, &j2, &col2));
    640:         if (col2 != col) break;
    641:         PetscCall(PetscHeapPop(h, &j2, &col2));
    642:         if (bb[j2] < bi[acol[j2] + 1]) PetscCall(PetscHeapStash(h, j2, bj[bb[j2]++]));
    643:       }
    644:       /* Put any stashed elements back into the min heap */
    645:       PetscCall(PetscHeapUnstash(h));
    646:       PetscCall(PetscHeapPop(h, &j, &col));
    647:     }
    648:   }
    649:   PetscCall(PetscFree(bb));
    650:   PetscCall(PetscHeapDestroy(&h));

    652:   /* Column indices are in the list of free space */
    653:   /* Allocate space for cj, initialize cj, and */
    654:   /* destroy list of free space and other temporary array(s) */
    655:   PetscCall(PetscMalloc1(ci[am], &cj));
    656:   PetscCall(PetscFreeSpaceContiguous(&free_space, cj));

    658:   /* put together the new symbolic matrix */
    659:   PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C));
    660:   PetscCall(MatSetBlockSizesFromMats(C, A, B));

    662:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    663:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    664:   c          = (Mat_SeqAIJ *)C->data;
    665:   c->free_a  = PETSC_TRUE;
    666:   c->free_ij = PETSC_TRUE;
    667:   c->nonew   = 0;

    669:   C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

    671:   /* set MatInfo */
    672:   afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
    673:   if (afill < 1.0) afill = 1.0;
    674:   C->info.mallocs           = ndouble;
    675:   C->info.fill_ratio_given  = fill;
    676:   C->info.fill_ratio_needed = afill;

    678: #if defined(PETSC_USE_INFO)
    679:   if (ci[am]) {
    680:     PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
    681:     PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
    682:   } else {
    683:     PetscCall(PetscInfo(C, "Empty matrix product\n"));
    684:   }
    685: #endif
    686:   PetscFunctionReturn(PETSC_SUCCESS);
    687: }

    689: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(Mat A, Mat B, PetscReal fill, Mat C)
    690: {
    691:   Mat_SeqAIJ        *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
    692:   const PetscInt    *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j;
    693:   PetscInt          *ci, *cj, *bb;
    694:   PetscInt           am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
    695:   PetscReal          afill;
    696:   PetscInt           i, j, col, ndouble = 0;
    697:   PetscFreeSpaceList free_space = NULL, current_space = NULL;
    698:   PetscHeap          h;
    699:   PetscBT            bt;

    701:   PetscFunctionBegin;
    702:   /* Get ci and cj - using a heap for the sorted rows, but use BT so that each index is only added once */
    703:   /* Allocate arrays for fill computation and free space for accumulating nonzero column */
    704:   PetscCall(PetscMalloc1(am + 2, &ci));
    705:   ci[0] = 0;

    707:   /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
    708:   PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space));

    710:   current_space = free_space;

    712:   PetscCall(PetscHeapCreate(a->rmax, &h));
    713:   PetscCall(PetscMalloc1(a->rmax, &bb));
    714:   PetscCall(PetscBTCreate(bn, &bt));

    716:   /* Determine ci and cj */
    717:   for (i = 0; i < am; i++) {
    718:     const PetscInt  anzi = ai[i + 1] - ai[i];    /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
    719:     const PetscInt *acol = aj + ai[i];           /* column indices of nonzero entries in this row */
    720:     const PetscInt *fptr = current_space->array; /* Save beginning of the row so we can clear the BT later */
    721:     ci[i + 1]            = ci[i];
    722:     /* Populate the min heap */
    723:     for (j = 0; j < anzi; j++) {
    724:       PetscInt brow = acol[j];
    725:       for (bb[j] = bi[brow]; bb[j] < bi[brow + 1]; bb[j]++) {
    726:         PetscInt bcol = bj[bb[j]];
    727:         if (!PetscBTLookupSet(bt, bcol)) { /* new entry */
    728:           PetscCall(PetscHeapAdd(h, j, bcol));
    729:           bb[j]++;
    730:           break;
    731:         }
    732:       }
    733:     }
    734:     /* Pick off the min element, adding it to free space */
    735:     PetscCall(PetscHeapPop(h, &j, &col));
    736:     while (j >= 0) {
    737:       if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
    738:         fptr = NULL;                            /* need PetscBTMemzero */
    739:         PetscCall(PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2, current_space->total_array_size), 16 << 20), &current_space));
    740:         ndouble++;
    741:       }
    742:       *(current_space->array++) = col;
    743:       current_space->local_used++;
    744:       current_space->local_remaining--;
    745:       ci[i + 1]++;

    747:       /* stash if anything else remains in this row of B */
    748:       for (; bb[j] < bi[acol[j] + 1]; bb[j]++) {
    749:         PetscInt bcol = bj[bb[j]];
    750:         if (!PetscBTLookupSet(bt, bcol)) { /* new entry */
    751:           PetscCall(PetscHeapAdd(h, j, bcol));
    752:           bb[j]++;
    753:           break;
    754:         }
    755:       }
    756:       PetscCall(PetscHeapPop(h, &j, &col));
    757:     }
    758:     if (fptr) { /* Clear the bits for this row */
    759:       for (; fptr < current_space->array; fptr++) PetscCall(PetscBTClear(bt, *fptr));
    760:     } else { /* We reallocated so we don't remember (easily) how to clear only the bits we changed */
    761:       PetscCall(PetscBTMemzero(bn, bt));
    762:     }
    763:   }
    764:   PetscCall(PetscFree(bb));
    765:   PetscCall(PetscHeapDestroy(&h));
    766:   PetscCall(PetscBTDestroy(&bt));

    768:   /* Column indices are in the list of free space */
    769:   /* Allocate space for cj, initialize cj, and */
    770:   /* destroy list of free space and other temporary array(s) */
    771:   PetscCall(PetscMalloc1(ci[am], &cj));
    772:   PetscCall(PetscFreeSpaceContiguous(&free_space, cj));

    774:   /* put together the new symbolic matrix */
    775:   PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C));
    776:   PetscCall(MatSetBlockSizesFromMats(C, A, B));

    778:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    779:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    780:   c          = (Mat_SeqAIJ *)C->data;
    781:   c->free_a  = PETSC_TRUE;
    782:   c->free_ij = PETSC_TRUE;
    783:   c->nonew   = 0;

    785:   C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

    787:   /* set MatInfo */
    788:   afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
    789:   if (afill < 1.0) afill = 1.0;
    790:   C->info.mallocs           = ndouble;
    791:   C->info.fill_ratio_given  = fill;
    792:   C->info.fill_ratio_needed = afill;

    794: #if defined(PETSC_USE_INFO)
    795:   if (ci[am]) {
    796:     PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
    797:     PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
    798:   } else {
    799:     PetscCall(PetscInfo(C, "Empty matrix product\n"));
    800:   }
    801: #endif
    802:   PetscFunctionReturn(PETSC_SUCCESS);
    803: }

    805: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(Mat A, Mat B, PetscReal fill, Mat C)
    806: {
    807:   Mat_SeqAIJ     *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
    808:   const PetscInt *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j, *inputi, *inputj, *inputcol, *inputcol_L1;
    809:   PetscInt       *ci, *cj, *outputj, worki_L1[9], worki_L2[9];
    810:   PetscInt        c_maxmem, a_maxrownnz = 0, a_rownnz;
    811:   const PetscInt  workcol[8] = {0, 1, 2, 3, 4, 5, 6, 7};
    812:   const PetscInt  am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
    813:   const PetscInt *brow_ptr[8], *brow_end[8];
    814:   PetscInt        window[8];
    815:   PetscInt        window_min, old_window_min, ci_nnz, outputi_nnz = 0, L1_nrows, L2_nrows;
    816:   PetscInt        i, k, ndouble = 0, L1_rowsleft, rowsleft;
    817:   PetscReal       afill;
    818:   PetscInt       *workj_L1, *workj_L2, *workj_L3;
    819:   PetscInt        L1_nnz, L2_nnz;

    821:   /* Step 1: Get upper bound on memory required for allocation.
    822:              Because of the way virtual memory works,
    823:              only the memory pages that are actually needed will be physically allocated. */
    824:   PetscFunctionBegin;
    825:   PetscCall(PetscMalloc1(am + 1, &ci));
    826:   for (i = 0; i < am; i++) {
    827:     const PetscInt  anzi = ai[i + 1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
    828:     const PetscInt *acol = aj + ai[i];        /* column indices of nonzero entries in this row */
    829:     a_rownnz             = 0;
    830:     for (k = 0; k < anzi; ++k) {
    831:       a_rownnz += bi[acol[k] + 1] - bi[acol[k]];
    832:       if (a_rownnz > bn) {
    833:         a_rownnz = bn;
    834:         break;
    835:       }
    836:     }
    837:     a_maxrownnz = PetscMax(a_maxrownnz, a_rownnz);
    838:   }
    839:   /* temporary work areas for merging rows */
    840:   PetscCall(PetscMalloc1(a_maxrownnz * 8, &workj_L1));
    841:   PetscCall(PetscMalloc1(a_maxrownnz * 8, &workj_L2));
    842:   PetscCall(PetscMalloc1(a_maxrownnz, &workj_L3));

    844:   /* This should be enough for almost all matrices. If not, memory is reallocated later. */
    845:   c_maxmem = 8 * (ai[am] + bi[bm]);
    846:   /* Step 2: Populate pattern for C */
    847:   PetscCall(PetscMalloc1(c_maxmem, &cj));

    849:   ci_nnz      = 0;
    850:   ci[0]       = 0;
    851:   worki_L1[0] = 0;
    852:   worki_L2[0] = 0;
    853:   for (i = 0; i < am; i++) {
    854:     const PetscInt  anzi = ai[i + 1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
    855:     const PetscInt *acol = aj + ai[i];        /* column indices of nonzero entries in this row */
    856:     rowsleft             = anzi;
    857:     inputcol_L1          = acol;
    858:     L2_nnz               = 0;
    859:     L2_nrows             = 1; /* Number of rows to be merged on Level 3. output of L3 already exists -> initial value 1   */
    860:     worki_L2[1]          = 0;
    861:     outputi_nnz          = 0;

    863:     /* If the number of indices in C so far + the max number of columns in the next row > c_maxmem  -> allocate more memory */
    864:     while (ci_nnz + a_maxrownnz > c_maxmem) {
    865:       c_maxmem *= 2;
    866:       ndouble++;
    867:       PetscCall(PetscRealloc(sizeof(PetscInt) * c_maxmem, &cj));
    868:     }

    870:     while (rowsleft) {
    871:       L1_rowsleft = PetscMin(64, rowsleft); /* In the inner loop max 64 rows of B can be merged */
    872:       L1_nrows    = 0;
    873:       L1_nnz      = 0;
    874:       inputcol    = inputcol_L1;
    875:       inputi      = bi;
    876:       inputj      = bj;

    878:       /* The following macro is used to specialize for small rows in A.
    879:          This helps with compiler unrolling, improving performance substantially.
    880:           Input:  inputj   inputi  inputcol  bn
    881:           Output: outputj  outputi_nnz                       */
    882: #define MatMatMultSymbolic_RowMergeMacro(ANNZ) \
    883:   do { \
    884:     window_min  = bn; \
    885:     outputi_nnz = 0; \
    886:     for (k = 0; k < ANNZ; ++k) { \
    887:       brow_ptr[k] = inputj + inputi[inputcol[k]]; \
    888:       brow_end[k] = inputj + inputi[inputcol[k] + 1]; \
    889:       window[k]   = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \
    890:       window_min  = PetscMin(window[k], window_min); \
    891:     } \
    892:     while (window_min < bn) { \
    893:       outputj[outputi_nnz++] = window_min; \
    894:       /* advance front and compute new minimum */ \
    895:       old_window_min = window_min; \
    896:       window_min     = bn; \
    897:       for (k = 0; k < ANNZ; ++k) { \
    898:         if (window[k] == old_window_min) { \
    899:           brow_ptr[k]++; \
    900:           window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \
    901:         } \
    902:         window_min = PetscMin(window[k], window_min); \
    903:       } \
    904:     } \
    905:   } while (0)

    907:       /************** L E V E L  1 ***************/
    908:       /* Merge up to 8 rows of B to L1 work array*/
    909:       while (L1_rowsleft) {
    910:         outputi_nnz = 0;
    911:         if (anzi > 8) outputj = workj_L1 + L1_nnz; /* Level 1 rowmerge*/
    912:         else outputj = cj + ci_nnz;                /* Merge directly to C */

    914:         switch (L1_rowsleft) {
    915:         case 1:
    916:           brow_ptr[0] = inputj + inputi[inputcol[0]];
    917:           brow_end[0] = inputj + inputi[inputcol[0] + 1];
    918:           for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
    919:           inputcol += L1_rowsleft;
    920:           rowsleft -= L1_rowsleft;
    921:           L1_rowsleft = 0;
    922:           break;
    923:         case 2:
    924:           MatMatMultSymbolic_RowMergeMacro(2);
    925:           inputcol += L1_rowsleft;
    926:           rowsleft -= L1_rowsleft;
    927:           L1_rowsleft = 0;
    928:           break;
    929:         case 3:
    930:           MatMatMultSymbolic_RowMergeMacro(3);
    931:           inputcol += L1_rowsleft;
    932:           rowsleft -= L1_rowsleft;
    933:           L1_rowsleft = 0;
    934:           break;
    935:         case 4:
    936:           MatMatMultSymbolic_RowMergeMacro(4);
    937:           inputcol += L1_rowsleft;
    938:           rowsleft -= L1_rowsleft;
    939:           L1_rowsleft = 0;
    940:           break;
    941:         case 5:
    942:           MatMatMultSymbolic_RowMergeMacro(5);
    943:           inputcol += L1_rowsleft;
    944:           rowsleft -= L1_rowsleft;
    945:           L1_rowsleft = 0;
    946:           break;
    947:         case 6:
    948:           MatMatMultSymbolic_RowMergeMacro(6);
    949:           inputcol += L1_rowsleft;
    950:           rowsleft -= L1_rowsleft;
    951:           L1_rowsleft = 0;
    952:           break;
    953:         case 7:
    954:           MatMatMultSymbolic_RowMergeMacro(7);
    955:           inputcol += L1_rowsleft;
    956:           rowsleft -= L1_rowsleft;
    957:           L1_rowsleft = 0;
    958:           break;
    959:         default:
    960:           MatMatMultSymbolic_RowMergeMacro(8);
    961:           inputcol += 8;
    962:           rowsleft -= 8;
    963:           L1_rowsleft -= 8;
    964:           break;
    965:         }
    966:         inputcol_L1 = inputcol;
    967:         L1_nnz += outputi_nnz;
    968:         worki_L1[++L1_nrows] = L1_nnz;
    969:       }

    971:       /********************** L E V E L  2 ************************/
    972:       /* Merge from L1 work array to either C or to L2 work array */
    973:       if (anzi > 8) {
    974:         inputi      = worki_L1;
    975:         inputj      = workj_L1;
    976:         inputcol    = workcol;
    977:         outputi_nnz = 0;

    979:         if (anzi <= 64) outputj = cj + ci_nnz; /* Merge from L1 work array to C */
    980:         else outputj = workj_L2 + L2_nnz;      /* Merge from L1 work array to L2 work array */

    982:         switch (L1_nrows) {
    983:         case 1:
    984:           brow_ptr[0] = inputj + inputi[inputcol[0]];
    985:           brow_end[0] = inputj + inputi[inputcol[0] + 1];
    986:           for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
    987:           break;
    988:         case 2:
    989:           MatMatMultSymbolic_RowMergeMacro(2);
    990:           break;
    991:         case 3:
    992:           MatMatMultSymbolic_RowMergeMacro(3);
    993:           break;
    994:         case 4:
    995:           MatMatMultSymbolic_RowMergeMacro(4);
    996:           break;
    997:         case 5:
    998:           MatMatMultSymbolic_RowMergeMacro(5);
    999:           break;
    1000:         case 6:
    1001:           MatMatMultSymbolic_RowMergeMacro(6);
    1002:           break;
    1003:         case 7:
    1004:           MatMatMultSymbolic_RowMergeMacro(7);
    1005:           break;
    1006:         case 8:
    1007:           MatMatMultSymbolic_RowMergeMacro(8);
    1008:           break;
    1009:         default:
    1010:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatMatMult logic error: Not merging 1-8 rows from L1 work array!");
    1011:         }
    1012:         L2_nnz += outputi_nnz;
    1013:         worki_L2[++L2_nrows] = L2_nnz;

    1015:         /************************ L E V E L  3 **********************/
    1016:         /* Merge from L2 work array to either C or to L2 work array */
    1017:         if (anzi > 64 && (L2_nrows == 8 || rowsleft == 0)) {
    1018:           inputi      = worki_L2;
    1019:           inputj      = workj_L2;
    1020:           inputcol    = workcol;
    1021:           outputi_nnz = 0;
    1022:           if (rowsleft) outputj = workj_L3;
    1023:           else outputj = cj + ci_nnz;
    1024:           switch (L2_nrows) {
    1025:           case 1:
    1026:             brow_ptr[0] = inputj + inputi[inputcol[0]];
    1027:             brow_end[0] = inputj + inputi[inputcol[0] + 1];
    1028:             for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
    1029:             break;
    1030:           case 2:
    1031:             MatMatMultSymbolic_RowMergeMacro(2);
    1032:             break;
    1033:           case 3:
    1034:             MatMatMultSymbolic_RowMergeMacro(3);
    1035:             break;
    1036:           case 4:
    1037:             MatMatMultSymbolic_RowMergeMacro(4);
    1038:             break;
    1039:           case 5:
    1040:             MatMatMultSymbolic_RowMergeMacro(5);
    1041:             break;
    1042:           case 6:
    1043:             MatMatMultSymbolic_RowMergeMacro(6);
    1044:             break;
    1045:           case 7:
    1046:             MatMatMultSymbolic_RowMergeMacro(7);
    1047:             break;
    1048:           case 8:
    1049:             MatMatMultSymbolic_RowMergeMacro(8);
    1050:             break;
    1051:           default:
    1052:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatMatMult logic error: Not merging 1-8 rows from L2 work array!");
    1053:           }
    1054:           L2_nrows    = 1;
    1055:           L2_nnz      = outputi_nnz;
    1056:           worki_L2[1] = outputi_nnz;
    1057:           /* Copy to workj_L2 */
    1058:           if (rowsleft) {
    1059:             for (k = 0; k < outputi_nnz; ++k) workj_L2[k] = outputj[k];
    1060:           }
    1061:         }
    1062:       }
    1063:     } /* while (rowsleft) */
    1064: #undef MatMatMultSymbolic_RowMergeMacro

    1066:     /* terminate current row */
    1067:     ci_nnz += outputi_nnz;
    1068:     ci[i + 1] = ci_nnz;
    1069:   }

    1071:   /* Step 3: Create the new symbolic matrix */
    1072:   PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C));
    1073:   PetscCall(MatSetBlockSizesFromMats(C, A, B));

    1075:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    1076:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    1077:   c          = (Mat_SeqAIJ *)C->data;
    1078:   c->free_a  = PETSC_TRUE;
    1079:   c->free_ij = PETSC_TRUE;
    1080:   c->nonew   = 0;

    1082:   C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

    1084:   /* set MatInfo */
    1085:   afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
    1086:   if (afill < 1.0) afill = 1.0;
    1087:   C->info.mallocs           = ndouble;
    1088:   C->info.fill_ratio_given  = fill;
    1089:   C->info.fill_ratio_needed = afill;

    1091: #if defined(PETSC_USE_INFO)
    1092:   if (ci[am]) {
    1093:     PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
    1094:     PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
    1095:   } else {
    1096:     PetscCall(PetscInfo(C, "Empty matrix product\n"));
    1097:   }
    1098: #endif

    1100:   /* Step 4: Free temporary work areas */
    1101:   PetscCall(PetscFree(workj_L1));
    1102:   PetscCall(PetscFree(workj_L2));
    1103:   PetscCall(PetscFree(workj_L3));
    1104:   PetscFunctionReturn(PETSC_SUCCESS);
    1105: }

    1107: /* concatenate unique entries and then sort */
    1108: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(Mat A, Mat B, PetscReal fill, Mat C)
    1109: {
    1110:   Mat_SeqAIJ     *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
    1111:   const PetscInt *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j;
    1112:   PetscInt       *ci, *cj, bcol;
    1113:   PetscInt        am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
    1114:   PetscReal       afill;
    1115:   PetscInt        i, j, ndouble = 0;
    1116:   PetscSegBuffer  seg, segrow;
    1117:   char           *seen;

    1119:   PetscFunctionBegin;
    1120:   PetscCall(PetscMalloc1(am + 1, &ci));
    1121:   ci[0] = 0;

    1123:   /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
    1124:   PetscCall(PetscSegBufferCreate(sizeof(PetscInt), (PetscInt)(fill * (ai[am] + bi[bm])), &seg));
    1125:   PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 100, &segrow));
    1126:   PetscCall(PetscCalloc1(bn, &seen));

    1128:   /* Determine ci and cj */
    1129:   for (i = 0; i < am; i++) {
    1130:     const PetscInt  anzi = ai[i + 1] - ai[i];                     /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
    1131:     const PetscInt *acol = PetscSafePointerPlusOffset(aj, ai[i]); /* column indices of nonzero entries in this row */
    1132:     PetscInt packlen     = 0, *PETSC_RESTRICT crow;

    1134:     /* Pack segrow */
    1135:     for (j = 0; j < anzi; j++) {
    1136:       PetscInt brow = acol[j], bjstart = bi[brow], bjend = bi[brow + 1], k;
    1137:       for (k = bjstart; k < bjend; k++) {
    1138:         bcol = bj[k];
    1139:         if (!seen[bcol]) { /* new entry */
    1140:           PetscInt *PETSC_RESTRICT slot;
    1141:           PetscCall(PetscSegBufferGetInts(segrow, 1, &slot));
    1142:           *slot      = bcol;
    1143:           seen[bcol] = 1;
    1144:           packlen++;
    1145:         }
    1146:       }
    1147:     }

    1149:     /* Check i-th diagonal entry */
    1150:     if (C->force_diagonals && !seen[i]) {
    1151:       PetscInt *PETSC_RESTRICT slot;
    1152:       PetscCall(PetscSegBufferGetInts(segrow, 1, &slot));
    1153:       *slot   = i;
    1154:       seen[i] = 1;
    1155:       packlen++;
    1156:     }

    1158:     PetscCall(PetscSegBufferGetInts(seg, packlen, &crow));
    1159:     PetscCall(PetscSegBufferExtractTo(segrow, crow));
    1160:     PetscCall(PetscSortInt(packlen, crow));
    1161:     ci[i + 1] = ci[i] + packlen;
    1162:     for (j = 0; j < packlen; j++) seen[crow[j]] = 0;
    1163:   }
    1164:   PetscCall(PetscSegBufferDestroy(&segrow));
    1165:   PetscCall(PetscFree(seen));

    1167:   /* Column indices are in the segmented buffer */
    1168:   PetscCall(PetscSegBufferExtractAlloc(seg, &cj));
    1169:   PetscCall(PetscSegBufferDestroy(&seg));

    1171:   /* put together the new symbolic matrix */
    1172:   PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C));
    1173:   PetscCall(MatSetBlockSizesFromMats(C, A, B));

    1175:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    1176:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    1177:   c          = (Mat_SeqAIJ *)C->data;
    1178:   c->free_a  = PETSC_TRUE;
    1179:   c->free_ij = PETSC_TRUE;
    1180:   c->nonew   = 0;

    1182:   C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

    1184:   /* set MatInfo */
    1185:   afill = (PetscReal)ci[am] / PetscMax(ai[am] + bi[bm], 1) + 1.e-5;
    1186:   if (afill < 1.0) afill = 1.0;
    1187:   C->info.mallocs           = ndouble;
    1188:   C->info.fill_ratio_given  = fill;
    1189:   C->info.fill_ratio_needed = afill;

    1191: #if defined(PETSC_USE_INFO)
    1192:   if (ci[am]) {
    1193:     PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
    1194:     PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
    1195:   } else {
    1196:     PetscCall(PetscInfo(C, "Empty matrix product\n"));
    1197:   }
    1198: #endif
    1199:   PetscFunctionReturn(PETSC_SUCCESS);
    1200: }

    1202: static PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(void *data)
    1203: {
    1204:   Mat_MatMatTransMult *abt = (Mat_MatMatTransMult *)data;

    1206:   PetscFunctionBegin;
    1207:   PetscCall(MatTransposeColoringDestroy(&abt->matcoloring));
    1208:   PetscCall(MatDestroy(&abt->Bt_den));
    1209:   PetscCall(MatDestroy(&abt->ABt_den));
    1210:   PetscCall(PetscFree(abt));
    1211:   PetscFunctionReturn(PETSC_SUCCESS);
    1212: }

    1214: PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C)
    1215: {
    1216:   Mat                  Bt;
    1217:   Mat_MatMatTransMult *abt;
    1218:   Mat_Product         *product = C->product;
    1219:   char                *alg;

    1221:   PetscFunctionBegin;
    1222:   PetscCheck(product, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct");
    1223:   PetscCheck(!product->data, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Extra product struct not empty");

    1225:   /* create symbolic Bt */
    1226:   PetscCall(MatTransposeSymbolic(B, &Bt));
    1227:   PetscCall(MatSetBlockSizes(Bt, A->cmap->bs, B->cmap->bs));
    1228:   PetscCall(MatSetType(Bt, ((PetscObject)A)->type_name));

    1230:   /* get symbolic C=A*Bt */
    1231:   PetscCall(PetscStrallocpy(product->alg, &alg));
    1232:   PetscCall(MatProductSetAlgorithm(C, "sorted")); /* set algorithm for C = A*Bt */
    1233:   PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(A, Bt, fill, C));
    1234:   PetscCall(MatProductSetAlgorithm(C, alg)); /* resume original algorithm for ABt product */
    1235:   PetscCall(PetscFree(alg));

    1237:   /* create a supporting struct for reuse intermediate dense matrices with matcoloring */
    1238:   PetscCall(PetscNew(&abt));

    1240:   product->data    = abt;
    1241:   product->destroy = MatDestroy_SeqAIJ_MatMatMultTrans;

    1243:   C->ops->mattransposemultnumeric = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ;

    1245:   abt->usecoloring = PETSC_FALSE;
    1246:   PetscCall(PetscStrcmp(product->alg, "color", &abt->usecoloring));
    1247:   if (abt->usecoloring) {
    1248:     /* Create MatTransposeColoring from symbolic C=A*B^T */
    1249:     MatTransposeColoring matcoloring;
    1250:     MatColoring          coloring;
    1251:     ISColoring           iscoloring;
    1252:     Mat                  Bt_dense, C_dense;

    1254:     /* inode causes memory problem */
    1255:     PetscCall(MatSetOption(C, MAT_USE_INODES, PETSC_FALSE));

    1257:     PetscCall(MatColoringCreate(C, &coloring));
    1258:     PetscCall(MatColoringSetDistance(coloring, 2));
    1259:     PetscCall(MatColoringSetType(coloring, MATCOLORINGSL));
    1260:     PetscCall(MatColoringSetFromOptions(coloring));
    1261:     PetscCall(MatColoringApply(coloring, &iscoloring));
    1262:     PetscCall(MatColoringDestroy(&coloring));
    1263:     PetscCall(MatTransposeColoringCreate(C, iscoloring, &matcoloring));

    1265:     abt->matcoloring = matcoloring;

    1267:     PetscCall(ISColoringDestroy(&iscoloring));

    1269:     /* Create Bt_dense and C_dense = A*Bt_dense */
    1270:     PetscCall(MatCreate(PETSC_COMM_SELF, &Bt_dense));
    1271:     PetscCall(MatSetSizes(Bt_dense, A->cmap->n, matcoloring->ncolors, A->cmap->n, matcoloring->ncolors));
    1272:     PetscCall(MatSetType(Bt_dense, MATSEQDENSE));
    1273:     PetscCall(MatSeqDenseSetPreallocation(Bt_dense, NULL));

    1275:     Bt_dense->assembled = PETSC_TRUE;
    1276:     abt->Bt_den         = Bt_dense;

    1278:     PetscCall(MatCreate(PETSC_COMM_SELF, &C_dense));
    1279:     PetscCall(MatSetSizes(C_dense, A->rmap->n, matcoloring->ncolors, A->rmap->n, matcoloring->ncolors));
    1280:     PetscCall(MatSetType(C_dense, MATSEQDENSE));
    1281:     PetscCall(MatSeqDenseSetPreallocation(C_dense, NULL));

    1283:     Bt_dense->assembled = PETSC_TRUE;
    1284:     abt->ABt_den        = C_dense;

    1286: #if defined(PETSC_USE_INFO)
    1287:     {
    1288:       Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data;
    1289:       PetscCall(PetscInfo(C, "Use coloring of C=A*B^T; B^T: %" PetscInt_FMT " %" PetscInt_FMT ", Bt_dense: %" PetscInt_FMT ",%" PetscInt_FMT "; Cnz %" PetscInt_FMT " / (cm*ncolors %" PetscInt_FMT ") = %g\n", B->cmap->n, B->rmap->n, Bt_dense->rmap->n,
    1290:                           Bt_dense->cmap->n, c->nz, A->rmap->n * matcoloring->ncolors, (double)(((PetscReal)c->nz) / ((PetscReal)(A->rmap->n * matcoloring->ncolors)))));
    1291:     }
    1292: #endif
    1293:   }
    1294:   /* clean up */
    1295:   PetscCall(MatDestroy(&Bt));
    1296:   PetscFunctionReturn(PETSC_SUCCESS);
    1297: }

    1299: PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C)
    1300: {
    1301:   Mat_SeqAIJ          *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c = (Mat_SeqAIJ *)C->data;
    1302:   PetscInt            *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, anzi, bnzj, nexta, nextb, *acol, *bcol, brow;
    1303:   PetscInt             cm = C->rmap->n, *ci = c->i, *cj = c->j, i, j, cnzi, *ccol;
    1304:   PetscLogDouble       flops = 0.0;
    1305:   MatScalar           *aa = a->a, *aval, *ba = b->a, *bval, *ca, *cval;
    1306:   Mat_MatMatTransMult *abt;
    1307:   Mat_Product         *product = C->product;

    1309:   PetscFunctionBegin;
    1310:   PetscCheck(product, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct");
    1311:   abt = (Mat_MatMatTransMult *)product->data;
    1312:   PetscCheck(abt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct");
    1313:   /* clear old values in C */
    1314:   if (!c->a) {
    1315:     PetscCall(PetscCalloc1(ci[cm] + 1, &ca));
    1316:     c->a      = ca;
    1317:     c->free_a = PETSC_TRUE;
    1318:   } else {
    1319:     ca = c->a;
    1320:     PetscCall(PetscArrayzero(ca, ci[cm] + 1));
    1321:   }

    1323:   if (abt->usecoloring) {
    1324:     MatTransposeColoring matcoloring = abt->matcoloring;
    1325:     Mat                  Bt_dense, C_dense = abt->ABt_den;

    1327:     /* Get Bt_dense by Apply MatTransposeColoring to B */
    1328:     Bt_dense = abt->Bt_den;
    1329:     PetscCall(MatTransColoringApplySpToDen(matcoloring, B, Bt_dense));

    1331:     /* C_dense = A*Bt_dense */
    1332:     PetscCall(MatMatMultNumeric_SeqAIJ_SeqDense(A, Bt_dense, C_dense));

    1334:     /* Recover C from C_dense */
    1335:     PetscCall(MatTransColoringApplyDenToSp(matcoloring, C_dense, C));
    1336:     PetscFunctionReturn(PETSC_SUCCESS);
    1337:   }

    1339:   for (i = 0; i < cm; i++) {
    1340:     anzi = ai[i + 1] - ai[i];
    1341:     acol = PetscSafePointerPlusOffset(aj, ai[i]);
    1342:     aval = PetscSafePointerPlusOffset(aa, ai[i]);
    1343:     cnzi = ci[i + 1] - ci[i];
    1344:     ccol = PetscSafePointerPlusOffset(cj, ci[i]);
    1345:     cval = ca + ci[i];
    1346:     for (j = 0; j < cnzi; j++) {
    1347:       brow = ccol[j];
    1348:       bnzj = bi[brow + 1] - bi[brow];
    1349:       bcol = bj + bi[brow];
    1350:       bval = ba + bi[brow];

    1352:       /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */
    1353:       nexta = 0;
    1354:       nextb = 0;
    1355:       while (nexta < anzi && nextb < bnzj) {
    1356:         while (nexta < anzi && acol[nexta] < bcol[nextb]) nexta++;
    1357:         if (nexta == anzi) break;
    1358:         while (nextb < bnzj && acol[nexta] > bcol[nextb]) nextb++;
    1359:         if (nextb == bnzj) break;
    1360:         if (acol[nexta] == bcol[nextb]) {
    1361:           cval[j] += aval[nexta] * bval[nextb];
    1362:           nexta++;
    1363:           nextb++;
    1364:           flops += 2;
    1365:         }
    1366:       }
    1367:     }
    1368:   }
    1369:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
    1370:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
    1371:   PetscCall(PetscLogFlops(flops));
    1372:   PetscFunctionReturn(PETSC_SUCCESS);
    1373: }

    1375: PetscErrorCode MatDestroy_SeqAIJ_MatTransMatMult(void *data)
    1376: {
    1377:   Mat_MatTransMatMult *atb = (Mat_MatTransMatMult *)data;

    1379:   PetscFunctionBegin;
    1380:   PetscCall(MatDestroy(&atb->At));
    1381:   if (atb->destroy) PetscCall((*atb->destroy)(atb->data));
    1382:   PetscCall(PetscFree(atb));
    1383:   PetscFunctionReturn(PETSC_SUCCESS);
    1384: }

    1386: PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C)
    1387: {
    1388:   Mat          At      = NULL;
    1389:   Mat_Product *product = C->product;
    1390:   PetscBool    flg, def, square;

    1392:   PetscFunctionBegin;
    1393:   MatCheckProduct(C, 4);
    1394:   square = (PetscBool)(A == B && A->symmetric == PETSC_BOOL3_TRUE);
    1395:   /* outerproduct */
    1396:   PetscCall(PetscStrcmp(product->alg, "outerproduct", &flg));
    1397:   if (flg) {
    1398:     /* create symbolic At */
    1399:     if (!square) {
    1400:       PetscCall(MatTransposeSymbolic(A, &At));
    1401:       PetscCall(MatSetBlockSizes(At, A->cmap->bs, B->cmap->bs));
    1402:       PetscCall(MatSetType(At, ((PetscObject)A)->type_name));
    1403:     }
    1404:     /* get symbolic C=At*B */
    1405:     PetscCall(MatProductSetAlgorithm(C, "sorted"));
    1406:     PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(square ? A : At, B, fill, C));

    1408:     /* clean up */
    1409:     if (!square) PetscCall(MatDestroy(&At));

    1411:     C->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ; /* outerproduct */
    1412:     PetscCall(MatProductSetAlgorithm(C, "outerproduct"));
    1413:     PetscFunctionReturn(PETSC_SUCCESS);
    1414:   }

    1416:   /* matmatmult */
    1417:   PetscCall(PetscStrcmp(product->alg, "default", &def));
    1418:   PetscCall(PetscStrcmp(product->alg, "at*b", &flg));
    1419:   if (flg || def) {
    1420:     Mat_MatTransMatMult *atb;

    1422:     PetscCheck(!product->data, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Extra product struct not empty");
    1423:     PetscCall(PetscNew(&atb));
    1424:     if (!square) PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &At));
    1425:     PetscCall(MatProductSetAlgorithm(C, "sorted"));
    1426:     PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(square ? A : At, B, fill, C));
    1427:     PetscCall(MatProductSetAlgorithm(C, "at*b"));
    1428:     product->data    = atb;
    1429:     product->destroy = MatDestroy_SeqAIJ_MatTransMatMult;
    1430:     atb->At          = At;

    1432:     C->ops->mattransposemultnumeric = NULL; /* see MatProductNumeric_AtB_SeqAIJ_SeqAIJ */
    1433:     PetscFunctionReturn(PETSC_SUCCESS);
    1434:   }

    1436:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Mat Product Algorithm is not supported");
    1437: }

    1439: PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C)
    1440: {
    1441:   Mat_SeqAIJ    *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c = (Mat_SeqAIJ *)C->data;
    1442:   PetscInt       am = A->rmap->n, anzi, *ai = a->i, *aj = a->j, *bi = b->i, *bj, bnzi, nextb;
    1443:   PetscInt       cm = C->rmap->n, *ci = c->i, *cj = c->j, crow, *cjj, i, j, k;
    1444:   PetscLogDouble flops = 0.0;
    1445:   MatScalar     *aa    = a->a, *ba, *ca, *caj;

    1447:   PetscFunctionBegin;
    1448:   if (!c->a) {
    1449:     PetscCall(PetscCalloc1(ci[cm] + 1, &ca));

    1451:     c->a      = ca;
    1452:     c->free_a = PETSC_TRUE;
    1453:   } else {
    1454:     ca = c->a;
    1455:     PetscCall(PetscArrayzero(ca, ci[cm]));
    1456:   }

    1458:   /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */
    1459:   for (i = 0; i < am; i++) {
    1460:     bj   = b->j + bi[i];
    1461:     ba   = b->a + bi[i];
    1462:     bnzi = bi[i + 1] - bi[i];
    1463:     anzi = ai[i + 1] - ai[i];
    1464:     for (j = 0; j < anzi; j++) {
    1465:       nextb = 0;
    1466:       crow  = *aj++;
    1467:       cjj   = cj + ci[crow];
    1468:       caj   = ca + ci[crow];
    1469:       /* perform sparse axpy operation.  Note cjj includes bj. */
    1470:       for (k = 0; nextb < bnzi; k++) {
    1471:         if (cjj[k] == *(bj + nextb)) { /* ccol == bcol */
    1472:           caj[k] += (*aa) * (*(ba + nextb));
    1473:           nextb++;
    1474:         }
    1475:       }
    1476:       flops += 2 * bnzi;
    1477:       aa++;
    1478:     }
    1479:   }

    1481:   /* Assemble the final matrix and clean up */
    1482:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
    1483:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
    1484:   PetscCall(PetscLogFlops(flops));
    1485:   PetscFunctionReturn(PETSC_SUCCESS);
    1486: }

    1488: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
    1489: {
    1490:   PetscFunctionBegin;
    1491:   PetscCall(MatMatMultSymbolic_SeqDense_SeqDense(A, B, 0.0, C));
    1492:   C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
    1493:   PetscFunctionReturn(PETSC_SUCCESS);
    1494: }

    1496: PETSC_INTERN PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A, Mat B, Mat C, const PetscBool add)
    1497: {
    1498:   Mat_SeqAIJ        *a = (Mat_SeqAIJ *)A->data;
    1499:   PetscScalar       *c, r1, r2, r3, r4, *c1, *c2, *c3, *c4;
    1500:   const PetscScalar *aa, *b, *b1, *b2, *b3, *b4, *av;
    1501:   const PetscInt    *aj;
    1502:   PetscInt           cm = C->rmap->n, cn = B->cmap->n, bm, am = A->rmap->n;
    1503:   PetscInt           clda;
    1504:   PetscInt           am4, bm4, col, i, j, n;

    1506:   PetscFunctionBegin;
    1507:   if (!cm || !cn) PetscFunctionReturn(PETSC_SUCCESS);
    1508:   PetscCall(MatSeqAIJGetArrayRead(A, &av));
    1509:   if (add) {
    1510:     PetscCall(MatDenseGetArray(C, &c));
    1511:   } else {
    1512:     PetscCall(MatDenseGetArrayWrite(C, &c));
    1513:   }
    1514:   PetscCall(MatDenseGetArrayRead(B, &b));
    1515:   PetscCall(MatDenseGetLDA(B, &bm));
    1516:   PetscCall(MatDenseGetLDA(C, &clda));
    1517:   am4 = 4 * clda;
    1518:   bm4 = 4 * bm;
    1519:   if (b) {
    1520:     b1 = b;
    1521:     b2 = b1 + bm;
    1522:     b3 = b2 + bm;
    1523:     b4 = b3 + bm;
    1524:   } else b1 = b2 = b3 = b4 = NULL;
    1525:   c1 = c;
    1526:   c2 = c1 + clda;
    1527:   c3 = c2 + clda;
    1528:   c4 = c3 + clda;
    1529:   for (col = 0; col < (cn / 4) * 4; col += 4) { /* over columns of C */
    1530:     for (i = 0; i < am; i++) {                  /* over rows of A in those columns */
    1531:       r1 = r2 = r3 = r4 = 0.0;
    1532:       n                 = a->i[i + 1] - a->i[i];
    1533:       aj                = PetscSafePointerPlusOffset(a->j, a->i[i]);
    1534:       aa                = PetscSafePointerPlusOffset(av, a->i[i]);
    1535:       for (j = 0; j < n; j++) {
    1536:         const PetscScalar aatmp = aa[j];
    1537:         const PetscInt    ajtmp = aj[j];
    1538:         r1 += aatmp * b1[ajtmp];
    1539:         r2 += aatmp * b2[ajtmp];
    1540:         r3 += aatmp * b3[ajtmp];
    1541:         r4 += aatmp * b4[ajtmp];
    1542:       }
    1543:       if (add) {
    1544:         c1[i] += r1;
    1545:         c2[i] += r2;
    1546:         c3[i] += r3;
    1547:         c4[i] += r4;
    1548:       } else {
    1549:         c1[i] = r1;
    1550:         c2[i] = r2;
    1551:         c3[i] = r3;
    1552:         c4[i] = r4;
    1553:       }
    1554:     }
    1555:     if (b) {
    1556:       b1 += bm4;
    1557:       b2 += bm4;
    1558:       b3 += bm4;
    1559:       b4 += bm4;
    1560:     }
    1561:     c1 += am4;
    1562:     c2 += am4;
    1563:     c3 += am4;
    1564:     c4 += am4;
    1565:   }
    1566:   /* process remaining columns */
    1567:   if (col != cn) {
    1568:     PetscInt rc = cn - col;

    1570:     if (rc == 1) {
    1571:       for (i = 0; i < am; i++) {
    1572:         r1 = 0.0;
    1573:         n  = a->i[i + 1] - a->i[i];
    1574:         aj = PetscSafePointerPlusOffset(a->j, a->i[i]);
    1575:         aa = PetscSafePointerPlusOffset(av, a->i[i]);
    1576:         for (j = 0; j < n; j++) r1 += aa[j] * b1[aj[j]];
    1577:         if (add) c1[i] += r1;
    1578:         else c1[i] = r1;
    1579:       }
    1580:     } else if (rc == 2) {
    1581:       for (i = 0; i < am; i++) {
    1582:         r1 = r2 = 0.0;
    1583:         n       = a->i[i + 1] - a->i[i];
    1584:         aj      = PetscSafePointerPlusOffset(a->j, a->i[i]);
    1585:         aa      = PetscSafePointerPlusOffset(av, a->i[i]);
    1586:         for (j = 0; j < n; j++) {
    1587:           const PetscScalar aatmp = aa[j];
    1588:           const PetscInt    ajtmp = aj[j];
    1589:           r1 += aatmp * b1[ajtmp];
    1590:           r2 += aatmp * b2[ajtmp];
    1591:         }
    1592:         if (add) {
    1593:           c1[i] += r1;
    1594:           c2[i] += r2;
    1595:         } else {
    1596:           c1[i] = r1;
    1597:           c2[i] = r2;
    1598:         }
    1599:       }
    1600:     } else {
    1601:       for (i = 0; i < am; i++) {
    1602:         r1 = r2 = r3 = 0.0;
    1603:         n            = a->i[i + 1] - a->i[i];
    1604:         aj           = PetscSafePointerPlusOffset(a->j, a->i[i]);
    1605:         aa           = PetscSafePointerPlusOffset(av, a->i[i]);
    1606:         for (j = 0; j < n; j++) {
    1607:           const PetscScalar aatmp = aa[j];
    1608:           const PetscInt    ajtmp = aj[j];
    1609:           r1 += aatmp * b1[ajtmp];
    1610:           r2 += aatmp * b2[ajtmp];
    1611:           r3 += aatmp * b3[ajtmp];
    1612:         }
    1613:         if (add) {
    1614:           c1[i] += r1;
    1615:           c2[i] += r2;
    1616:           c3[i] += r3;
    1617:         } else {
    1618:           c1[i] = r1;
    1619:           c2[i] = r2;
    1620:           c3[i] = r3;
    1621:         }
    1622:       }
    1623:     }
    1624:   }
    1625:   PetscCall(PetscLogFlops(cn * (2.0 * a->nz)));
    1626:   if (add) {
    1627:     PetscCall(MatDenseRestoreArray(C, &c));
    1628:   } else {
    1629:     PetscCall(MatDenseRestoreArrayWrite(C, &c));
    1630:   }
    1631:   PetscCall(MatDenseRestoreArrayRead(B, &b));
    1632:   PetscCall(MatSeqAIJRestoreArrayRead(A, &av));
    1633:   PetscFunctionReturn(PETSC_SUCCESS);
    1634: }

    1636: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A, Mat B, Mat C)
    1637: {
    1638:   PetscFunctionBegin;
    1639:   PetscCheck(B->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number columns in A %" PetscInt_FMT " not equal rows in B %" PetscInt_FMT, A->cmap->n, B->rmap->n);
    1640:   PetscCheck(A->rmap->n == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number rows in C %" PetscInt_FMT " not equal rows in A %" PetscInt_FMT, C->rmap->n, A->rmap->n);
    1641:   PetscCheck(B->cmap->n == C->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number columns in B %" PetscInt_FMT " not equal columns in C %" PetscInt_FMT, B->cmap->n, C->cmap->n);

    1643:   PetscCall(MatMatMultNumericAdd_SeqAIJ_SeqDense(A, B, C, PETSC_FALSE));
    1644:   PetscFunctionReturn(PETSC_SUCCESS);
    1645: }

    1647: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_AB(Mat C)
    1648: {
    1649:   PetscFunctionBegin;
    1650:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqAIJ_SeqDense;
    1651:   C->ops->productsymbolic = MatProductSymbolic_AB;
    1652:   PetscFunctionReturn(PETSC_SUCCESS);
    1653: }

    1655: PETSC_INTERN PetscErrorCode MatTMatTMultSymbolic_SeqAIJ_SeqDense(Mat, Mat, PetscReal, Mat);

    1657: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_AtB(Mat C)
    1658: {
    1659:   PetscFunctionBegin;
    1660:   C->ops->transposematmultsymbolic = MatTMatTMultSymbolic_SeqAIJ_SeqDense;
    1661:   C->ops->productsymbolic          = MatProductSymbolic_AtB;
    1662:   PetscFunctionReturn(PETSC_SUCCESS);
    1663: }

    1665: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_ABt(Mat C)
    1666: {
    1667:   PetscFunctionBegin;
    1668:   C->ops->mattransposemultsymbolic = MatTMatTMultSymbolic_SeqAIJ_SeqDense;
    1669:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
    1670:   PetscFunctionReturn(PETSC_SUCCESS);
    1671: }

    1673: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense(Mat C)
    1674: {
    1675:   Mat_Product *product = C->product;

    1677:   PetscFunctionBegin;
    1678:   switch (product->type) {
    1679:   case MATPRODUCT_AB:
    1680:     PetscCall(MatProductSetFromOptions_SeqAIJ_SeqDense_AB(C));
    1681:     break;
    1682:   case MATPRODUCT_AtB:
    1683:     PetscCall(MatProductSetFromOptions_SeqAIJ_SeqDense_AtB(C));
    1684:     break;
    1685:   case MATPRODUCT_ABt:
    1686:     PetscCall(MatProductSetFromOptions_SeqAIJ_SeqDense_ABt(C));
    1687:     break;
    1688:   default:
    1689:     break;
    1690:   }
    1691:   PetscFunctionReturn(PETSC_SUCCESS);
    1692: }

    1694: static PetscErrorCode MatProductSetFromOptions_SeqXBAIJ_SeqDense_AB(Mat C)
    1695: {
    1696:   Mat_Product *product = C->product;
    1697:   Mat          A       = product->A;
    1698:   PetscBool    baij;

    1700:   PetscFunctionBegin;
    1701:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQBAIJ, &baij));
    1702:   if (!baij) { /* A is seqsbaij */
    1703:     PetscBool sbaij;
    1704:     PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &sbaij));
    1705:     PetscCheck(sbaij, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Mat must be either seqbaij or seqsbaij format");

    1707:     C->ops->matmultsymbolic = MatMatMultSymbolic_SeqSBAIJ_SeqDense;
    1708:   } else { /* A is seqbaij */
    1709:     C->ops->matmultsymbolic = MatMatMultSymbolic_SeqBAIJ_SeqDense;
    1710:   }

    1712:   C->ops->productsymbolic = MatProductSymbolic_AB;
    1713:   PetscFunctionReturn(PETSC_SUCCESS);
    1714: }

    1716: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqXBAIJ_SeqDense(Mat C)
    1717: {
    1718:   Mat_Product *product = C->product;

    1720:   PetscFunctionBegin;
    1721:   MatCheckProduct(C, 1);
    1722:   PetscCheck(product->A, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing A");
    1723:   if (product->type == MATPRODUCT_AB || (product->type == MATPRODUCT_AtB && product->A->symmetric == PETSC_BOOL3_TRUE)) PetscCall(MatProductSetFromOptions_SeqXBAIJ_SeqDense_AB(C));
    1724:   PetscFunctionReturn(PETSC_SUCCESS);
    1725: }

    1727: static PetscErrorCode MatProductSetFromOptions_SeqDense_SeqAIJ_AB(Mat C)
    1728: {
    1729:   PetscFunctionBegin;
    1730:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqAIJ;
    1731:   C->ops->productsymbolic = MatProductSymbolic_AB;
    1732:   PetscFunctionReturn(PETSC_SUCCESS);
    1733: }

    1735: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense_SeqAIJ(Mat C)
    1736: {
    1737:   Mat_Product *product = C->product;

    1739:   PetscFunctionBegin;
    1740:   if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_SeqDense_SeqAIJ_AB(C));
    1741:   PetscFunctionReturn(PETSC_SUCCESS);
    1742: }

    1744: PetscErrorCode MatTransColoringApplySpToDen_SeqAIJ(MatTransposeColoring coloring, Mat B, Mat Btdense)
    1745: {
    1746:   Mat_SeqAIJ   *b       = (Mat_SeqAIJ *)B->data;
    1747:   Mat_SeqDense *btdense = (Mat_SeqDense *)Btdense->data;
    1748:   PetscInt     *bi = b->i, *bj = b->j;
    1749:   PetscInt      m = Btdense->rmap->n, n = Btdense->cmap->n, j, k, l, col, anz, *btcol, brow, ncolumns;
    1750:   MatScalar    *btval, *btval_den, *ba = b->a;
    1751:   PetscInt     *columns = coloring->columns, *colorforcol = coloring->colorforcol, ncolors = coloring->ncolors;

    1753:   PetscFunctionBegin;
    1754:   btval_den = btdense->v;
    1755:   PetscCall(PetscArrayzero(btval_den, m * n));
    1756:   for (k = 0; k < ncolors; k++) {
    1757:     ncolumns = coloring->ncolumns[k];
    1758:     for (l = 0; l < ncolumns; l++) { /* insert a row of B to a column of Btdense */
    1759:       col   = *(columns + colorforcol[k] + l);
    1760:       btcol = bj + bi[col];
    1761:       btval = ba + bi[col];
    1762:       anz   = bi[col + 1] - bi[col];
    1763:       for (j = 0; j < anz; j++) {
    1764:         brow            = btcol[j];
    1765:         btval_den[brow] = btval[j];
    1766:       }
    1767:     }
    1768:     btval_den += m;
    1769:   }
    1770:   PetscFunctionReturn(PETSC_SUCCESS);
    1771: }

    1773: PetscErrorCode MatTransColoringApplyDenToSp_SeqAIJ(MatTransposeColoring matcoloring, Mat Cden, Mat Csp)
    1774: {
    1775:   Mat_SeqAIJ        *csp = (Mat_SeqAIJ *)Csp->data;
    1776:   const PetscScalar *ca_den, *ca_den_ptr;
    1777:   PetscScalar       *ca = csp->a;
    1778:   PetscInt           k, l, m = Cden->rmap->n, ncolors = matcoloring->ncolors;
    1779:   PetscInt           brows = matcoloring->brows, *den2sp = matcoloring->den2sp;
    1780:   PetscInt           nrows, *row, *idx;
    1781:   PetscInt          *rows = matcoloring->rows, *colorforrow = matcoloring->colorforrow;

    1783:   PetscFunctionBegin;
    1784:   PetscCall(MatDenseGetArrayRead(Cden, &ca_den));

    1786:   if (brows > 0) {
    1787:     PetscInt *lstart, row_end, row_start;
    1788:     lstart = matcoloring->lstart;
    1789:     PetscCall(PetscArrayzero(lstart, ncolors));

    1791:     row_end = brows;
    1792:     if (row_end > m) row_end = m;
    1793:     for (row_start = 0; row_start < m; row_start += brows) { /* loop over row blocks of Csp */
    1794:       ca_den_ptr = ca_den;
    1795:       for (k = 0; k < ncolors; k++) { /* loop over colors (columns of Cden) */
    1796:         nrows = matcoloring->nrows[k];
    1797:         row   = rows + colorforrow[k];
    1798:         idx   = den2sp + colorforrow[k];
    1799:         for (l = lstart[k]; l < nrows; l++) {
    1800:           if (row[l] >= row_end) {
    1801:             lstart[k] = l;
    1802:             break;
    1803:           } else {
    1804:             ca[idx[l]] = ca_den_ptr[row[l]];
    1805:           }
    1806:         }
    1807:         ca_den_ptr += m;
    1808:       }
    1809:       row_end += brows;
    1810:       if (row_end > m) row_end = m;
    1811:     }
    1812:   } else { /* non-blocked impl: loop over columns of Csp - slow if Csp is large */
    1813:     ca_den_ptr = ca_den;
    1814:     for (k = 0; k < ncolors; k++) {
    1815:       nrows = matcoloring->nrows[k];
    1816:       row   = rows + colorforrow[k];
    1817:       idx   = den2sp + colorforrow[k];
    1818:       for (l = 0; l < nrows; l++) ca[idx[l]] = ca_den_ptr[row[l]];
    1819:       ca_den_ptr += m;
    1820:     }
    1821:   }

    1823:   PetscCall(MatDenseRestoreArrayRead(Cden, &ca_den));
    1824: #if defined(PETSC_USE_INFO)
    1825:   if (matcoloring->brows > 0) {
    1826:     PetscCall(PetscInfo(Csp, "Loop over %" PetscInt_FMT " row blocks for den2sp\n", brows));
    1827:   } else {
    1828:     PetscCall(PetscInfo(Csp, "Loop over colors/columns of Cden, inefficient for large sparse matrix product \n"));
    1829:   }
    1830: #endif
    1831:   PetscFunctionReturn(PETSC_SUCCESS);
    1832: }

    1834: PetscErrorCode MatTransposeColoringCreate_SeqAIJ(Mat mat, ISColoring iscoloring, MatTransposeColoring c)
    1835: {
    1836:   PetscInt        i, n, nrows, Nbs, j, k, m, ncols, col, cm;
    1837:   const PetscInt *is, *ci, *cj, *row_idx;
    1838:   PetscInt        nis = iscoloring->n, *rowhit, bs = 1;
    1839:   IS             *isa;
    1840:   Mat_SeqAIJ     *csp = (Mat_SeqAIJ *)mat->data;
    1841:   PetscInt       *colorforrow, *rows, *rows_i, *idxhit, *spidx, *den2sp, *den2sp_i;
    1842:   PetscInt       *colorforcol, *columns, *columns_i, brows;
    1843:   PetscBool       flg;

    1845:   PetscFunctionBegin;
    1846:   PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, PETSC_IGNORE, &isa));

    1848:   /* bs > 1 is not being tested yet! */
    1849:   Nbs       = mat->cmap->N / bs;
    1850:   c->M      = mat->rmap->N / bs; /* set total rows, columns and local rows */
    1851:   c->N      = Nbs;
    1852:   c->m      = c->M;
    1853:   c->rstart = 0;
    1854:   c->brows  = 100;

    1856:   c->ncolors = nis;
    1857:   PetscCall(PetscMalloc3(nis, &c->ncolumns, nis, &c->nrows, nis + 1, &colorforrow));
    1858:   PetscCall(PetscMalloc1(csp->nz + 1, &rows));
    1859:   PetscCall(PetscMalloc1(csp->nz + 1, &den2sp));

    1861:   brows = c->brows;
    1862:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-matden2sp_brows", &brows, &flg));
    1863:   if (flg) c->brows = brows;
    1864:   if (brows > 0) PetscCall(PetscMalloc1(nis + 1, &c->lstart));

    1866:   colorforrow[0] = 0;
    1867:   rows_i         = rows;
    1868:   den2sp_i       = den2sp;

    1870:   PetscCall(PetscMalloc1(nis + 1, &colorforcol));
    1871:   PetscCall(PetscMalloc1(Nbs + 1, &columns));

    1873:   colorforcol[0] = 0;
    1874:   columns_i      = columns;

    1876:   /* get column-wise storage of mat */
    1877:   PetscCall(MatGetColumnIJ_SeqAIJ_Color(mat, 0, PETSC_FALSE, PETSC_FALSE, &ncols, &ci, &cj, &spidx, NULL));

    1879:   cm = c->m;
    1880:   PetscCall(PetscMalloc1(cm + 1, &rowhit));
    1881:   PetscCall(PetscMalloc1(cm + 1, &idxhit));
    1882:   for (i = 0; i < nis; i++) { /* loop over color */
    1883:     PetscCall(ISGetLocalSize(isa[i], &n));
    1884:     PetscCall(ISGetIndices(isa[i], &is));

    1886:     c->ncolumns[i] = n;
    1887:     if (n) PetscCall(PetscArraycpy(columns_i, is, n));
    1888:     colorforcol[i + 1] = colorforcol[i] + n;
    1889:     columns_i += n;

    1891:     /* fast, crude version requires O(N*N) work */
    1892:     PetscCall(PetscArrayzero(rowhit, cm));

    1894:     for (j = 0; j < n; j++) { /* loop over columns*/
    1895:       col     = is[j];
    1896:       row_idx = cj + ci[col];
    1897:       m       = ci[col + 1] - ci[col];
    1898:       for (k = 0; k < m; k++) { /* loop over columns marking them in rowhit */
    1899:         idxhit[*row_idx]   = spidx[ci[col] + k];
    1900:         rowhit[*row_idx++] = col + 1;
    1901:       }
    1902:     }
    1903:     /* count the number of hits */
    1904:     nrows = 0;
    1905:     for (j = 0; j < cm; j++) {
    1906:       if (rowhit[j]) nrows++;
    1907:     }
    1908:     c->nrows[i]        = nrows;
    1909:     colorforrow[i + 1] = colorforrow[i] + nrows;

    1911:     nrows = 0;
    1912:     for (j = 0; j < cm; j++) { /* loop over rows */
    1913:       if (rowhit[j]) {
    1914:         rows_i[nrows]   = j;
    1915:         den2sp_i[nrows] = idxhit[j];
    1916:         nrows++;
    1917:       }
    1918:     }
    1919:     den2sp_i += nrows;

    1921:     PetscCall(ISRestoreIndices(isa[i], &is));
    1922:     rows_i += nrows;
    1923:   }
    1924:   PetscCall(MatRestoreColumnIJ_SeqAIJ_Color(mat, 0, PETSC_FALSE, PETSC_FALSE, &ncols, &ci, &cj, &spidx, NULL));
    1925:   PetscCall(PetscFree(rowhit));
    1926:   PetscCall(ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &isa));
    1927:   PetscCheck(csp->nz == colorforrow[nis], PETSC_COMM_SELF, PETSC_ERR_PLIB, "csp->nz %" PetscInt_FMT " != colorforrow[nis] %" PetscInt_FMT, csp->nz, colorforrow[nis]);

    1929:   c->colorforrow = colorforrow;
    1930:   c->rows        = rows;
    1931:   c->den2sp      = den2sp;
    1932:   c->colorforcol = colorforcol;
    1933:   c->columns     = columns;

    1935:   PetscCall(PetscFree(idxhit));
    1936:   PetscFunctionReturn(PETSC_SUCCESS);
    1937: }

    1939: static PetscErrorCode MatProductNumeric_AtB_SeqAIJ_SeqAIJ(Mat C)
    1940: {
    1941:   Mat_Product *product = C->product;
    1942:   Mat          A = product->A, B = product->B;

    1944:   PetscFunctionBegin;
    1945:   if (C->ops->mattransposemultnumeric) {
    1946:     /* Alg: "outerproduct" */
    1947:     PetscCall((*C->ops->mattransposemultnumeric)(A, B, C));
    1948:   } else {
    1949:     /* Alg: "matmatmult" -- C = At*B */
    1950:     Mat_MatTransMatMult *atb = (Mat_MatTransMatMult *)product->data;

    1952:     PetscCheck(atb, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct");
    1953:     if (atb->At) {
    1954:       /* At is computed in MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ();
    1955:          user may have called MatProductReplaceMats() to get this A=product->A */
    1956:       PetscCall(MatTransposeSetPrecursor(A, atb->At));
    1957:       PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &atb->At));
    1958:     }
    1959:     PetscCall(MatMatMultNumeric_SeqAIJ_SeqAIJ(atb->At ? atb->At : A, B, C));
    1960:   }
    1961:   PetscFunctionReturn(PETSC_SUCCESS);
    1962: }

    1964: static PetscErrorCode MatProductSymbolic_AtB_SeqAIJ_SeqAIJ(Mat C)
    1965: {
    1966:   Mat_Product *product = C->product;
    1967:   Mat          A = product->A, B = product->B;
    1968:   PetscReal    fill = product->fill;

    1970:   PetscFunctionBegin;
    1971:   PetscCall(MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(A, B, fill, C));

    1973:   C->ops->productnumeric = MatProductNumeric_AtB_SeqAIJ_SeqAIJ;
    1974:   PetscFunctionReturn(PETSC_SUCCESS);
    1975: }

    1977: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_AB(Mat C)
    1978: {
    1979:   Mat_Product *product = C->product;
    1980:   PetscInt     alg     = 0; /* default algorithm */
    1981:   PetscBool    flg     = PETSC_FALSE;
    1982: #if !defined(PETSC_HAVE_HYPRE)
    1983:   const char *algTypes[7] = {"sorted", "scalable", "scalable_fast", "heap", "btheap", "llcondensed", "rowmerge"};
    1984:   PetscInt    nalg        = 7;
    1985: #else
    1986:   const char *algTypes[8] = {"sorted", "scalable", "scalable_fast", "heap", "btheap", "llcondensed", "rowmerge", "hypre"};
    1987:   PetscInt    nalg        = 8;
    1988: #endif

    1990:   PetscFunctionBegin;
    1991:   /* Set default algorithm */
    1992:   PetscCall(PetscStrcmp(C->product->alg, "default", &flg));
    1993:   if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg]));

    1995:   /* Get runtime option */
    1996:   if (product->api_user) {
    1997:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatMatMult", "Mat");
    1998:     PetscCall(PetscOptionsEList("-matmatmult_via", "Algorithmic approach", "MatMatMult", algTypes, nalg, algTypes[0], &alg, &flg));
    1999:     PetscOptionsEnd();
    2000:   } else {
    2001:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_AB", "Mat");
    2002:     PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_AB", algTypes, nalg, algTypes[0], &alg, &flg));
    2003:     PetscOptionsEnd();
    2004:   }
    2005:   if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg]));

    2007:   C->ops->productsymbolic = MatProductSymbolic_AB;
    2008:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqAIJ_SeqAIJ;
    2009:   PetscFunctionReturn(PETSC_SUCCESS);
    2010: }

    2012: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_AtB(Mat C)
    2013: {
    2014:   Mat_Product *product     = C->product;
    2015:   PetscInt     alg         = 0; /* default algorithm */
    2016:   PetscBool    flg         = PETSC_FALSE;
    2017:   const char  *algTypes[3] = {"default", "at*b", "outerproduct"};
    2018:   PetscInt     nalg        = 3;

    2020:   PetscFunctionBegin;
    2021:   /* Get runtime option */
    2022:   if (product->api_user) {
    2023:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatTransposeMatMult", "Mat");
    2024:     PetscCall(PetscOptionsEList("-mattransposematmult_via", "Algorithmic approach", "MatTransposeMatMult", algTypes, nalg, algTypes[alg], &alg, &flg));
    2025:     PetscOptionsEnd();
    2026:   } else {
    2027:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_AtB", "Mat");
    2028:     PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_AtB", algTypes, nalg, algTypes[alg], &alg, &flg));
    2029:     PetscOptionsEnd();
    2030:   }
    2031:   if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg]));

    2033:   C->ops->productsymbolic = MatProductSymbolic_AtB_SeqAIJ_SeqAIJ;
    2034:   PetscFunctionReturn(PETSC_SUCCESS);
    2035: }

    2037: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_ABt(Mat C)
    2038: {
    2039:   Mat_Product *product     = C->product;
    2040:   PetscInt     alg         = 0; /* default algorithm */
    2041:   PetscBool    flg         = PETSC_FALSE;
    2042:   const char  *algTypes[2] = {"default", "color"};
    2043:   PetscInt     nalg        = 2;

    2045:   PetscFunctionBegin;
    2046:   /* Set default algorithm */
    2047:   PetscCall(PetscStrcmp(C->product->alg, "default", &flg));
    2048:   if (!flg) {
    2049:     alg = 1;
    2050:     PetscCall(MatProductSetAlgorithm(C, algTypes[alg]));
    2051:   }

    2053:   /* Get runtime option */
    2054:   if (product->api_user) {
    2055:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatMatTransposeMult", "Mat");
    2056:     PetscCall(PetscOptionsEList("-matmattransmult_via", "Algorithmic approach", "MatMatTransposeMult", algTypes, nalg, algTypes[alg], &alg, &flg));
    2057:     PetscOptionsEnd();
    2058:   } else {
    2059:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_ABt", "Mat");
    2060:     PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_ABt", algTypes, nalg, algTypes[alg], &alg, &flg));
    2061:     PetscOptionsEnd();
    2062:   }
    2063:   if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg]));

    2065:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ;
    2066:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
    2067:   PetscFunctionReturn(PETSC_SUCCESS);
    2068: }

    2070: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_PtAP(Mat C)
    2071: {
    2072:   Mat_Product *product = C->product;
    2073:   PetscBool    flg     = PETSC_FALSE;
    2074:   PetscInt     alg     = 0; /* default algorithm -- alg=1 should be default!!! */
    2075: #if !defined(PETSC_HAVE_HYPRE)
    2076:   const char *algTypes[2] = {"scalable", "rap"};
    2077:   PetscInt    nalg        = 2;
    2078: #else
    2079:   const char *algTypes[3] = {"scalable", "rap", "hypre"};
    2080:   PetscInt    nalg        = 3;
    2081: #endif

    2083:   PetscFunctionBegin;
    2084:   /* Set default algorithm */
    2085:   PetscCall(PetscStrcmp(product->alg, "default", &flg));
    2086:   if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg]));

    2088:   /* Get runtime option */
    2089:   if (product->api_user) {
    2090:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatPtAP", "Mat");
    2091:     PetscCall(PetscOptionsEList("-matptap_via", "Algorithmic approach", "MatPtAP", algTypes, nalg, algTypes[0], &alg, &flg));
    2092:     PetscOptionsEnd();
    2093:   } else {
    2094:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_PtAP", "Mat");
    2095:     PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_PtAP", algTypes, nalg, algTypes[0], &alg, &flg));
    2096:     PetscOptionsEnd();
    2097:   }
    2098:   if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg]));

    2100:   C->ops->productsymbolic = MatProductSymbolic_PtAP_SeqAIJ_SeqAIJ;
    2101:   PetscFunctionReturn(PETSC_SUCCESS);
    2102: }

    2104: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_RARt(Mat C)
    2105: {
    2106:   Mat_Product *product     = C->product;
    2107:   PetscBool    flg         = PETSC_FALSE;
    2108:   PetscInt     alg         = 0; /* default algorithm */
    2109:   const char  *algTypes[3] = {"r*a*rt", "r*art", "coloring_rart"};
    2110:   PetscInt     nalg        = 3;

    2112:   PetscFunctionBegin;
    2113:   /* Set default algorithm */
    2114:   PetscCall(PetscStrcmp(product->alg, "default", &flg));
    2115:   if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg]));

    2117:   /* Get runtime option */
    2118:   if (product->api_user) {
    2119:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatRARt", "Mat");
    2120:     PetscCall(PetscOptionsEList("-matrart_via", "Algorithmic approach", "MatRARt", algTypes, nalg, algTypes[0], &alg, &flg));
    2121:     PetscOptionsEnd();
    2122:   } else {
    2123:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_RARt", "Mat");
    2124:     PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_RARt", algTypes, nalg, algTypes[0], &alg, &flg));
    2125:     PetscOptionsEnd();
    2126:   }
    2127:   if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg]));

    2129:   C->ops->productsymbolic = MatProductSymbolic_RARt_SeqAIJ_SeqAIJ;
    2130:   PetscFunctionReturn(PETSC_SUCCESS);
    2131: }

    2133: /* ABC = A*B*C = A*(B*C); ABC's algorithm must be chosen from AB's algorithm */
    2134: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_ABC(Mat C)
    2135: {
    2136:   Mat_Product *product     = C->product;
    2137:   PetscInt     alg         = 0; /* default algorithm */
    2138:   PetscBool    flg         = PETSC_FALSE;
    2139:   const char  *algTypes[7] = {"sorted", "scalable", "scalable_fast", "heap", "btheap", "llcondensed", "rowmerge"};
    2140:   PetscInt     nalg        = 7;

    2142:   PetscFunctionBegin;
    2143:   /* Set default algorithm */
    2144:   PetscCall(PetscStrcmp(product->alg, "default", &flg));
    2145:   if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg]));

    2147:   /* Get runtime option */
    2148:   if (product->api_user) {
    2149:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatMatMatMult", "Mat");
    2150:     PetscCall(PetscOptionsEList("-matmatmatmult_via", "Algorithmic approach", "MatMatMatMult", algTypes, nalg, algTypes[alg], &alg, &flg));
    2151:     PetscOptionsEnd();
    2152:   } else {
    2153:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_ABC", "Mat");
    2154:     PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_ABC", algTypes, nalg, algTypes[alg], &alg, &flg));
    2155:     PetscOptionsEnd();
    2156:   }
    2157:   if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg]));

    2159:   C->ops->matmatmultsymbolic = MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ;
    2160:   C->ops->productsymbolic    = MatProductSymbolic_ABC;
    2161:   PetscFunctionReturn(PETSC_SUCCESS);
    2162: }

    2164: PetscErrorCode MatProductSetFromOptions_SeqAIJ(Mat C)
    2165: {
    2166:   Mat_Product *product = C->product;

    2168:   PetscFunctionBegin;
    2169:   switch (product->type) {
    2170:   case MATPRODUCT_AB:
    2171:     PetscCall(MatProductSetFromOptions_SeqAIJ_AB(C));
    2172:     break;
    2173:   case MATPRODUCT_AtB:
    2174:     PetscCall(MatProductSetFromOptions_SeqAIJ_AtB(C));
    2175:     break;
    2176:   case MATPRODUCT_ABt:
    2177:     PetscCall(MatProductSetFromOptions_SeqAIJ_ABt(C));
    2178:     break;
    2179:   case MATPRODUCT_PtAP:
    2180:     PetscCall(MatProductSetFromOptions_SeqAIJ_PtAP(C));
    2181:     break;
    2182:   case MATPRODUCT_RARt:
    2183:     PetscCall(MatProductSetFromOptions_SeqAIJ_RARt(C));
    2184:     break;
    2185:   case MATPRODUCT_ABC:
    2186:     PetscCall(MatProductSetFromOptions_SeqAIJ_ABC(C));
    2187:     break;
    2188:   default:
    2189:     break;
    2190:   }
    2191:   PetscFunctionReturn(PETSC_SUCCESS);
    2192: }