Actual source code: matmatmult.c

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

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

     13: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
     14: {

     18:   if (C->ops->matmultnumeric) {
     19:     if (C->ops->matmultnumeric == MatMatMultNumeric_SeqAIJ_SeqAIJ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Recursive call");
     20:     (*C->ops->matmultnumeric)(A,B,C);
     21:   } else {
     22:     MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(A,B,C);
     23:   }
     24:   return(0);
     25: }

     27: /* Modified from MatCreateSeqAIJWithArrays() */
     28: PETSC_INTERN PetscErrorCode MatSetSeqAIJWithArrays_private(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],MatType mtype,Mat mat)
     29: {
     31:   PetscInt       ii;
     32:   Mat_SeqAIJ     *aij;
     33:   PetscBool      isseqaij;

     36:   if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0");
     37:   MatSetSizes(mat,m,n,m,n);

     39:   if (!mtype) {
     40:     PetscObjectBaseTypeCompare((PetscObject)mat,MATSEQAIJ,&isseqaij);
     41:     if (!isseqaij) { MatSetType(mat,MATSEQAIJ); }
     42:   } else {
     43:     MatSetType(mat,mtype);
     44:   }
     45:   MatSeqAIJSetPreallocation_SeqAIJ(mat,MAT_SKIP_ALLOCATION,NULL);
     46:   aij  = (Mat_SeqAIJ*)(mat)->data;
     47:   PetscMalloc1(m,&aij->imax);
     48:   PetscMalloc1(m,&aij->ilen);

     50:   aij->i            = i;
     51:   aij->j            = j;
     52:   aij->a            = a;
     53:   aij->singlemalloc = PETSC_FALSE;
     54:   aij->nonew        = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
     55:   aij->free_a       = PETSC_FALSE;
     56:   aij->free_ij      = PETSC_FALSE;

     58:   for (ii=0; ii<m; ii++) {
     59:     aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii];
     60:   }

     62:   return(0);
     63: }

     65: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C)
     66: {
     67:   PetscErrorCode      ierr;
     68:   Mat_Product         *product = C->product;
     69:   MatProductAlgorithm alg;
     70:   PetscBool           flg;

     73:   if (product) {
     74:     alg = product->alg;
     75:   } else {
     76:     alg = "sorted";
     77:   }
     78:   /* sorted */
     79:   PetscStrcmp(alg,"sorted",&flg);
     80:   if (flg) {
     81:     MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(A,B,fill,C);
     82:     return(0);
     83:   }

     85:   /* scalable */
     86:   PetscStrcmp(alg,"scalable",&flg);
     87:   if (flg) {
     88:     MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A,B,fill,C);
     89:     return(0);
     90:   }

     92:   /* scalable_fast */
     93:   PetscStrcmp(alg,"scalable_fast",&flg);
     94:   if (flg) {
     95:     MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(A,B,fill,C);
     96:     return(0);
     97:   }

     99:   /* heap */
    100:   PetscStrcmp(alg,"heap",&flg);
    101:   if (flg) {
    102:     MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(A,B,fill,C);
    103:     return(0);
    104:   }

    106:   /* btheap */
    107:   PetscStrcmp(alg,"btheap",&flg);
    108:   if (flg) {
    109:     MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(A,B,fill,C);
    110:     return(0);
    111:   }

    113:   /* llcondensed */
    114:   PetscStrcmp(alg,"llcondensed",&flg);
    115:   if (flg) {
    116:     MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(A,B,fill,C);
    117:     return(0);
    118:   }

    120:   /* rowmerge */
    121:   PetscStrcmp(alg,"rowmerge",&flg);
    122:   if (flg) {
    123:     MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(A,B,fill,C);
    124:     return(0);
    125:   }

    127: #if defined(PETSC_HAVE_HYPRE)
    128:   PetscStrcmp(alg,"hypre",&flg);
    129:   if (flg) {
    130:     MatMatMultSymbolic_AIJ_AIJ_wHYPRE(A,B,fill,C);
    131:     return(0);
    132:   }
    133: #endif

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

    138: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(Mat A,Mat B,PetscReal fill,Mat C)
    139: {
    140:   PetscErrorCode     ierr;
    141:   Mat_SeqAIJ         *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
    142:   PetscInt           *ai=a->i,*bi=b->i,*ci,*cj;
    143:   PetscInt           am =A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
    144:   PetscReal          afill;
    145:   PetscInt           i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;
    146:   PetscTable         ta;
    147:   PetscBT            lnkbt;
    148:   PetscFreeSpaceList free_space=NULL,current_space=NULL;

    151:   /* Get ci and cj */
    152:   /*---------------*/
    153:   /* Allocate ci array, arrays for fill computation and */
    154:   /* free space for accumulating nonzero column info */
    155:   PetscMalloc1(am+2,&ci);
    156:   ci[0] = 0;

    158:   /* create and initialize a linked list */
    159:   PetscTableCreate(bn,bn,&ta);
    160:   MatRowMergeMax_SeqAIJ(b,bm,ta);
    161:   PetscTableGetCount(ta,&Crmax);
    162:   PetscTableDestroy(&ta);

    164:   PetscLLCondensedCreate(Crmax,bn,&lnk,&lnkbt);

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

    169:   current_space = free_space;

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

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

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

    198:     current_space->array           += cnzi;
    199:     current_space->local_used      += cnzi;
    200:     current_space->local_remaining -= cnzi;

    202:     ci[i+1] = ci[i] + cnzi;
    203:   }

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

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

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

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

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

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

    246: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(Mat A,Mat B,Mat C)
    247: {
    248:   PetscErrorCode    ierr;
    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;

    262:   MatSeqAIJGetArrayRead(A,&aa);
    263:   MatSeqAIJGetArrayRead(B,&ba);
    264:   if (!c->a) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
    265:     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:   PetscObjectQuery((PetscObject)C,"__PETSc__ab_dense",(PetscObject*)&cab_dense);
    274:   if (!cab_dense) {
    275:     PetscMalloc1(B->cmap->N,&ab_dense);
    276:     PetscContainerCreate(PETSC_COMM_SELF,&cab_dense);
    277:     PetscContainerSetPointer(cab_dense,ab_dense);
    278:     PetscContainerSetUserDestroy(cab_dense,PetscContainerUserDestroyDefault);
    279:     PetscObjectCompose((PetscObject)C,"__PETSc__ab_dense",(PetscObject)cab_dense);
    280:     PetscObjectDereference((PetscObject)cab_dense);
    281:   }
    282:   PetscContainerGetPointer(cab_dense,(void**)&ab_dense);
    283:   PetscArrayzero(ab_dense,B->cmap->N);

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

    306:     cnzi = ci[i+1] - ci[i];
    307:     for (k=0; k<cnzi; k++) {
    308:       ca[k]          += ab_dense[cj[k]];
    309:       ab_dense[cj[k]] = 0.0; /* zero ab_dense */
    310:     }
    311:     flops += cnzi;
    312:     cj    += cnzi; ca += cnzi;
    313:   }
    314: #if defined(PETSC_HAVE_DEVICE)
    315:   if (C->offloadmask != PETSC_OFFLOAD_UNALLOCATED) C->offloadmask = PETSC_OFFLOAD_CPU;
    316: #endif
    317:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    318:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
    319:   PetscLogFlops(flops);
    320:   MatSeqAIJRestoreArrayRead(A,&aa);
    321:   MatSeqAIJRestoreArrayRead(B,&ba);
    322:   return(0);
    323: }

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

    340:   MatSeqAIJGetArrayRead(A,&aa);
    341:   MatSeqAIJGetArrayRead(B,&ba);
    342:   if (!ca) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
    343:     PetscMalloc1(ci[cm]+1,&ca);
    344:     c->a      = ca;
    345:     c->free_a = PETSC_TRUE;
    346:   }

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

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

    398:   /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_fast() */
    399:   /*-----------------------------------------------------------------------------------------*/
    400:   /* Allocate arrays for fill computation and free space for accumulating nonzero column */
    401:   PetscMalloc1(am+2,&ci);
    402:   ci[0] = 0;

    404:   /* create and initialize a linked list */
    405:   PetscTableCreate(bn,bn,&ta);
    406:   MatRowMergeMax_SeqAIJ(b,bm,ta);
    407:   PetscTableGetCount(ta,&Crmax);
    408:   PetscTableDestroy(&ta);

    410:   PetscLLCondensedCreate_fast(Crmax,&lnk);

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

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

    433:     /* If free space is not available, make more free space */
    434:     /* Double the amount of total space in the list */
    435:     if (current_space->local_remaining<cnzi) {
    436:       PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),&current_space);
    437:       ndouble++;
    438:     }

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

    443:     current_space->array           += cnzi;
    444:     current_space->local_used      += cnzi;
    445:     current_space->local_remaining -= cnzi;

    447:     ci[i+1] = ci[i] + cnzi;
    448:   }

    450:   /* Column indices are in the list of free space */
    451:   /* Allocate space for cj, initialize cj, and */
    452:   /* destroy list of free space and other temporary array(s) */
    453:   PetscMalloc1(ci[am]+1,&cj);
    454:   PetscFreeSpaceContiguous(&free_space,cj);
    455:   PetscLLCondensedDestroy_fast(lnk);

    457:   /* Allocate space for ca */
    458:   PetscCalloc1(ci[am]+1,&ca);

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

    464:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    465:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    466:   c          = (Mat_SeqAIJ*)(C->data);
    467:   c->free_a  = PETSC_TRUE;
    468:   c->free_ij = PETSC_TRUE;
    469:   c->nonew   = 0;

    471:   /* slower, less memory */
    472:   C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable;

    474:   /* set MatInfo */
    475:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
    476:   if (afill < 1.0) afill = 1.0;
    477:   c->maxnz                     = ci[am];
    478:   c->nz                        = ci[am];
    479:   C->info.mallocs           = ndouble;
    480:   C->info.fill_ratio_given  = fill;
    481:   C->info.fill_ratio_needed = afill;

    483: #if defined(PETSC_USE_INFO)
    484:   if (ci[am]) {
    485:     PetscInfo3(C,"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
    486:     PetscInfo1(C,"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
    487:   } else {
    488:     PetscInfo(C,"Empty matrix product\n");
    489:   }
    490: #endif
    491:   return(0);
    492: }

    494: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(Mat A,Mat B,PetscReal fill,Mat C)
    495: {
    496:   PetscErrorCode     ierr;
    497:   Mat_SeqAIJ         *a  = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
    498:   PetscInt           *ai = a->i,*bi=b->i,*ci,*cj;
    499:   PetscInt           am  = A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
    500:   MatScalar          *ca;
    501:   PetscReal          afill;
    502:   PetscInt           i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;
    503:   PetscTable         ta;
    504:   PetscFreeSpaceList free_space=NULL,current_space=NULL;

    507:   /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_Scalalbe() */
    508:   /*---------------------------------------------------------------------------------------------*/
    509:   /* Allocate arrays for fill computation and free space for accumulating nonzero column */
    510:   PetscMalloc1(am+2,&ci);
    511:   ci[0] = 0;

    513:   /* create and initialize a linked list */
    514:   PetscTableCreate(bn,bn,&ta);
    515:   MatRowMergeMax_SeqAIJ(b,bm,ta);
    516:   PetscTableGetCount(ta,&Crmax);
    517:   PetscTableDestroy(&ta);
    518:   PetscLLCondensedCreate_Scalable(Crmax,&lnk);

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

    524:   /* Determine ci and cj */
    525:   for (i=0; i<am; i++) {
    526:     anzi = ai[i+1] - ai[i];
    527:     aj   = a->j + ai[i];
    528:     for (j=0; j<anzi; j++) {
    529:       brow = aj[j];
    530:       bnzj = bi[brow+1] - bi[brow];
    531:       bj   = b->j + bi[brow];
    532:       /* add non-zero cols of B into the sorted linked list lnk */
    533:       PetscLLCondensedAddSorted_Scalable(bnzj,bj,lnk);
    534:     }
    535:     /* add possible missing diagonal entry */
    536:     if (C->force_diagonals) {
    537:       PetscLLCondensedAddSorted_Scalable(1,&i,lnk);
    538:     }

    540:     cnzi = lnk[0];

    542:     /* If free space is not available, make more free space */
    543:     /* Double the amount of total space in the list */
    544:     if (current_space->local_remaining<cnzi) {
    545:       PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),&current_space);
    546:       ndouble++;
    547:     }

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

    552:     current_space->array           += cnzi;
    553:     current_space->local_used      += cnzi;
    554:     current_space->local_remaining -= cnzi;

    556:     ci[i+1] = ci[i] + cnzi;
    557:   }

    559:   /* Column indices are in the list of free space */
    560:   /* Allocate space for cj, initialize cj, and */
    561:   /* destroy list of free space and other temporary array(s) */
    562:   PetscMalloc1(ci[am]+1,&cj);
    563:   PetscFreeSpaceContiguous(&free_space,cj);
    564:   PetscLLCondensedDestroy_Scalable(lnk);

    566:   /* Allocate space for ca */
    567:   /*-----------------------*/
    568:   PetscCalloc1(ci[am]+1,&ca);

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

    574:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    575:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    576:   c          = (Mat_SeqAIJ*)(C->data);
    577:   c->free_a  = PETSC_TRUE;
    578:   c->free_ij = PETSC_TRUE;
    579:   c->nonew   = 0;

    581:   /* slower, less memory */
    582:   C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable;

    584:   /* set MatInfo */
    585:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
    586:   if (afill < 1.0) afill = 1.0;
    587:   c->maxnz                     = ci[am];
    588:   c->nz                        = ci[am];
    589:   C->info.mallocs           = ndouble;
    590:   C->info.fill_ratio_given  = fill;
    591:   C->info.fill_ratio_needed = afill;

    593: #if defined(PETSC_USE_INFO)
    594:   if (ci[am]) {
    595:     PetscInfo3(C,"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
    596:     PetscInfo1(C,"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
    597:   } else {
    598:     PetscInfo(C,"Empty matrix product\n");
    599:   }
    600: #endif
    601:   return(0);
    602: }

    604: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(Mat A,Mat B,PetscReal fill,Mat C)
    605: {
    606:   PetscErrorCode     ierr;
    607:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
    608:   const PetscInt     *ai=a->i,*bi=b->i,*aj=a->j,*bj=b->j;
    609:   PetscInt           *ci,*cj,*bb;
    610:   PetscInt           am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
    611:   PetscReal          afill;
    612:   PetscInt           i,j,col,ndouble = 0;
    613:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
    614:   PetscHeap          h;

    617:   /* Get ci and cj - by merging sorted rows using a heap */
    618:   /*---------------------------------------------------------------------------------------------*/
    619:   /* Allocate arrays for fill computation and free space for accumulating nonzero column */
    620:   PetscMalloc1(am+2,&ci);
    621:   ci[0] = 0;

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

    627:   PetscHeapCreate(a->rmax,&h);
    628:   PetscMalloc1(a->rmax,&bb);

    630:   /* Determine ci and cj */
    631:   for (i=0; i<am; i++) {
    632:     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 */
    633:     const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
    634:     ci[i+1] = ci[i];
    635:     /* Populate the min heap */
    636:     for (j=0; j<anzi; j++) {
    637:       bb[j] = bi[acol[j]];         /* bb points at the start of the row */
    638:       if (bb[j] < bi[acol[j]+1]) { /* Add if row is nonempty */
    639:         PetscHeapAdd(h,j,bj[bb[j]++]);
    640:       }
    641:     }
    642:     /* Pick off the min element, adding it to free space */
    643:     PetscHeapPop(h,&j,&col);
    644:     while (j >= 0) {
    645:       if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
    646:         PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2,current_space->total_array_size),16 << 20),&current_space);
    647:         ndouble++;
    648:       }
    649:       *(current_space->array++) = col;
    650:       current_space->local_used++;
    651:       current_space->local_remaining--;
    652:       ci[i+1]++;

    654:       /* stash if anything else remains in this row of B */
    655:       if (bb[j] < bi[acol[j]+1]) {PetscHeapStash(h,j,bj[bb[j]++]);}
    656:       while (1) {               /* pop and stash any other rows of B that also had an entry in this column */
    657:         PetscInt j2,col2;
    658:         PetscHeapPeek(h,&j2,&col2);
    659:         if (col2 != col) break;
    660:         PetscHeapPop(h,&j2,&col2);
    661:         if (bb[j2] < bi[acol[j2]+1]) {PetscHeapStash(h,j2,bj[bb[j2]++]);}
    662:       }
    663:       /* Put any stashed elements back into the min heap */
    664:       PetscHeapUnstash(h);
    665:       PetscHeapPop(h,&j,&col);
    666:     }
    667:   }
    668:   PetscFree(bb);
    669:   PetscHeapDestroy(&h);

    671:   /* Column indices are in the list of free space */
    672:   /* Allocate space for cj, initialize cj, and */
    673:   /* destroy list of free space and other temporary array(s) */
    674:   PetscMalloc1(ci[am],&cj);
    675:   PetscFreeSpaceContiguous(&free_space,cj);

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

    681:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    682:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    683:   c          = (Mat_SeqAIJ*)(C->data);
    684:   c->free_a  = PETSC_TRUE;
    685:   c->free_ij = PETSC_TRUE;
    686:   c->nonew   = 0;

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

    690:   /* set MatInfo */
    691:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
    692:   if (afill < 1.0) afill = 1.0;
    693:   c->maxnz                     = ci[am];
    694:   c->nz                        = ci[am];
    695:   C->info.mallocs           = ndouble;
    696:   C->info.fill_ratio_given  = fill;
    697:   C->info.fill_ratio_needed = afill;

    699: #if defined(PETSC_USE_INFO)
    700:   if (ci[am]) {
    701:     PetscInfo3(C,"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
    702:     PetscInfo1(C,"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
    703:   } else {
    704:     PetscInfo(C,"Empty matrix product\n");
    705:   }
    706: #endif
    707:   return(0);
    708: }

    710: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(Mat A,Mat B,PetscReal fill,Mat C)
    711: {
    712:   PetscErrorCode     ierr;
    713:   Mat_SeqAIJ         *a  = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
    714:   const PetscInt     *ai = a->i,*bi=b->i,*aj=a->j,*bj=b->j;
    715:   PetscInt           *ci,*cj,*bb;
    716:   PetscInt           am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
    717:   PetscReal          afill;
    718:   PetscInt           i,j,col,ndouble = 0;
    719:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
    720:   PetscHeap          h;
    721:   PetscBT            bt;

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

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

    733:   current_space = free_space;

    735:   PetscHeapCreate(a->rmax,&h);
    736:   PetscMalloc1(a->rmax,&bb);
    737:   PetscBTCreate(bn,&bt);

    739:   /* Determine ci and cj */
    740:   for (i=0; i<am; i++) {
    741:     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 */
    742:     const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
    743:     const PetscInt *fptr = current_space->array; /* Save beginning of the row so we can clear the BT later */
    744:     ci[i+1] = ci[i];
    745:     /* Populate the min heap */
    746:     for (j=0; j<anzi; j++) {
    747:       PetscInt brow = acol[j];
    748:       for (bb[j] = bi[brow]; bb[j] < bi[brow+1]; bb[j]++) {
    749:         PetscInt bcol = bj[bb[j]];
    750:         if (!PetscBTLookupSet(bt,bcol)) { /* new entry */
    751:           PetscHeapAdd(h,j,bcol);
    752:           bb[j]++;
    753:           break;
    754:         }
    755:       }
    756:     }
    757:     /* Pick off the min element, adding it to free space */
    758:     PetscHeapPop(h,&j,&col);
    759:     while (j >= 0) {
    760:       if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
    761:         fptr = NULL;                      /* need PetscBTMemzero */
    762:         PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2,current_space->total_array_size),16 << 20),&current_space);
    763:         ndouble++;
    764:       }
    765:       *(current_space->array++) = col;
    766:       current_space->local_used++;
    767:       current_space->local_remaining--;
    768:       ci[i+1]++;

    770:       /* stash if anything else remains in this row of B */
    771:       for (; bb[j] < bi[acol[j]+1]; bb[j]++) {
    772:         PetscInt bcol = bj[bb[j]];
    773:         if (!PetscBTLookupSet(bt,bcol)) { /* new entry */
    774:           PetscHeapAdd(h,j,bcol);
    775:           bb[j]++;
    776:           break;
    777:         }
    778:       }
    779:       PetscHeapPop(h,&j,&col);
    780:     }
    781:     if (fptr) {                 /* Clear the bits for this row */
    782:       for (; fptr<current_space->array; fptr++) {PetscBTClear(bt,*fptr);}
    783:     } else {                    /* We reallocated so we don't remember (easily) how to clear only the bits we changed */
    784:       PetscBTMemzero(bn,bt);
    785:     }
    786:   }
    787:   PetscFree(bb);
    788:   PetscHeapDestroy(&h);
    789:   PetscBTDestroy(&bt);

    791:   /* Column indices are in the list of free space */
    792:   /* Allocate space for cj, initialize cj, and */
    793:   /* destroy list of free space and other temporary array(s) */
    794:   PetscMalloc1(ci[am],&cj);
    795:   PetscFreeSpaceContiguous(&free_space,cj);

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

    801:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    802:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    803:   c          = (Mat_SeqAIJ*)(C->data);
    804:   c->free_a  = PETSC_TRUE;
    805:   c->free_ij = PETSC_TRUE;
    806:   c->nonew   = 0;

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

    810:   /* set MatInfo */
    811:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
    812:   if (afill < 1.0) afill = 1.0;
    813:   c->maxnz                     = ci[am];
    814:   c->nz                        = ci[am];
    815:   C->info.mallocs           = ndouble;
    816:   C->info.fill_ratio_given  = fill;
    817:   C->info.fill_ratio_needed = afill;

    819: #if defined(PETSC_USE_INFO)
    820:   if (ci[am]) {
    821:     PetscInfo3(C,"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
    822:     PetscInfo1(C,"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
    823:   } else {
    824:     PetscInfo(C,"Empty matrix product\n");
    825:   }
    826: #endif
    827:   return(0);
    828: }

    831: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(Mat A,Mat B,PetscReal fill,Mat C)
    832: {
    833:   PetscErrorCode     ierr;
    834:   Mat_SeqAIJ         *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
    835:   const PetscInt     *ai=a->i,*bi=b->i,*aj=a->j,*bj=b->j,*inputi,*inputj,*inputcol,*inputcol_L1;
    836:   PetscInt           *ci,*cj,*outputj,worki_L1[9],worki_L2[9];
    837:   PetscInt           c_maxmem,a_maxrownnz=0,a_rownnz;
    838:   const PetscInt     workcol[8]={0,1,2,3,4,5,6,7};
    839:   const PetscInt     am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
    840:   const PetscInt     *brow_ptr[8],*brow_end[8];
    841:   PetscInt           window[8];
    842:   PetscInt           window_min,old_window_min,ci_nnz,outputi_nnz=0,L1_nrows,L2_nrows;
    843:   PetscInt           i,k,ndouble=0,L1_rowsleft,rowsleft;
    844:   PetscReal          afill;
    845:   PetscInt           *workj_L1,*workj_L2,*workj_L3;
    846:   PetscInt           L1_nnz,L2_nnz;

    848:   /* Step 1: Get upper bound on memory required for allocation.
    849:              Because of the way virtual memory works,
    850:              only the memory pages that are actually needed will be physically allocated. */
    852:   PetscMalloc1(am+1,&ci);
    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:     a_rownnz = 0;
    857:     for (k=0; k<anzi; ++k) {
    858:       a_rownnz += bi[acol[k]+1] - bi[acol[k]];
    859:       if (a_rownnz > bn) {
    860:         a_rownnz = bn;
    861:         break;
    862:       }
    863:     }
    864:     a_maxrownnz = PetscMax(a_maxrownnz, a_rownnz);
    865:   }
    866:   /* temporary work areas for merging rows */
    867:   PetscMalloc1(a_maxrownnz*8,&workj_L1);
    868:   PetscMalloc1(a_maxrownnz*8,&workj_L2);
    869:   PetscMalloc1(a_maxrownnz,&workj_L3);

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

    876:   ci_nnz       = 0;
    877:   ci[0]        = 0;
    878:   worki_L1[0]  = 0;
    879:   worki_L2[0]  = 0;
    880:   for (i=0; i<am; i++) {
    881:     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 */
    882:     const PetscInt *acol = aj + ai[i];      /* column indices of nonzero entries in this row */
    883:     rowsleft             = anzi;
    884:     inputcol_L1          = acol;
    885:     L2_nnz               = 0;
    886:     L2_nrows             = 1;  /* Number of rows to be merged on Level 3. output of L3 already exists -> initial value 1   */
    887:     worki_L2[1]          = 0;
    888:     outputi_nnz          = 0;

    890:     /* If the number of indices in C so far + the max number of columns in the next row > c_maxmem  -> allocate more memory */
    891:     while (ci_nnz+a_maxrownnz > c_maxmem) {
    892:       c_maxmem *= 2;
    893:       ndouble++;
    894:       PetscRealloc(sizeof(PetscInt)*c_maxmem,&cj);
    895:     }

    897:     while (rowsleft) {
    898:       L1_rowsleft = PetscMin(64, rowsleft); /* In the inner loop max 64 rows of B can be merged */
    899:       L1_nrows    = 0;
    900:       L1_nnz      = 0;
    901:       inputcol    = inputcol_L1;
    902:       inputi      = bi;
    903:       inputj      = bj;

    905:       /* The following macro is used to specialize for small rows in A.
    906:          This helps with compiler unrolling, improving performance substantially.
    907:           Input:  inputj   inputi  inputcol  bn
    908:           Output: outputj  outputi_nnz                       */
    909:        #define MatMatMultSymbolic_RowMergeMacro(ANNZ)                        \
    910:          window_min  = bn;                                                   \
    911:          outputi_nnz = 0;                                                    \
    912:          for (k=0; k<ANNZ; ++k) {                                            \
    913:            brow_ptr[k] = inputj + inputi[inputcol[k]];                       \
    914:            brow_end[k] = inputj + inputi[inputcol[k]+1];                     \
    915:            window[k]   = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn;   \
    916:            window_min  = PetscMin(window[k], window_min);                    \
    917:          }                                                                   \
    918:          while (window_min < bn) {                                           \
    919:            outputj[outputi_nnz++] = window_min;                              \
    920:            /* advance front and compute new minimum */                       \
    921:            old_window_min = window_min;                                      \
    922:            window_min = bn;                                                  \
    923:            for (k=0; k<ANNZ; ++k) {                                          \
    924:              if (window[k] == old_window_min) {                              \
    925:                brow_ptr[k]++;                                                \
    926:                window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \
    927:              }                                                               \
    928:              window_min = PetscMin(window[k], window_min);                   \
    929:            }                                                                 \
    930:          }

    932:       /************** L E V E L  1 ***************/
    933:       /* Merge up to 8 rows of B to L1 work array*/
    934:       while (L1_rowsleft) {
    935:         outputi_nnz = 0;
    936:         if (anzi > 8)  outputj = workj_L1 + L1_nnz;     /* Level 1 rowmerge*/
    937:         else           outputj = cj + ci_nnz;           /* Merge directly to C */

    939:         switch (L1_rowsleft) {
    940:         case 1:  brow_ptr[0] = inputj + inputi[inputcol[0]];
    941:                  brow_end[0] = inputj + inputi[inputcol[0]+1];
    942:                  for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
    943:                  inputcol    += L1_rowsleft;
    944:                  rowsleft    -= L1_rowsleft;
    945:                  L1_rowsleft  = 0;
    946:                  break;
    947:         case 2:  MatMatMultSymbolic_RowMergeMacro(2);
    948:                  inputcol    += L1_rowsleft;
    949:                  rowsleft    -= L1_rowsleft;
    950:                  L1_rowsleft  = 0;
    951:                  break;
    952:         case 3: MatMatMultSymbolic_RowMergeMacro(3);
    953:                  inputcol    += L1_rowsleft;
    954:                  rowsleft    -= L1_rowsleft;
    955:                  L1_rowsleft  = 0;
    956:                  break;
    957:         case 4:  MatMatMultSymbolic_RowMergeMacro(4);
    958:                  inputcol    += L1_rowsleft;
    959:                  rowsleft    -= L1_rowsleft;
    960:                  L1_rowsleft  = 0;
    961:                  break;
    962:         case 5:  MatMatMultSymbolic_RowMergeMacro(5);
    963:                  inputcol    += L1_rowsleft;
    964:                  rowsleft    -= L1_rowsleft;
    965:                  L1_rowsleft  = 0;
    966:                  break;
    967:         case 6:  MatMatMultSymbolic_RowMergeMacro(6);
    968:                  inputcol    += L1_rowsleft;
    969:                  rowsleft    -= L1_rowsleft;
    970:                  L1_rowsleft  = 0;
    971:                  break;
    972:         case 7:  MatMatMultSymbolic_RowMergeMacro(7);
    973:                  inputcol    += L1_rowsleft;
    974:                  rowsleft    -= L1_rowsleft;
    975:                  L1_rowsleft  = 0;
    976:                  break;
    977:         default: MatMatMultSymbolic_RowMergeMacro(8);
    978:                  inputcol    += 8;
    979:                  rowsleft    -= 8;
    980:                  L1_rowsleft -= 8;
    981:                  break;
    982:         }
    983:         inputcol_L1           = inputcol;
    984:         L1_nnz               += outputi_nnz;
    985:         worki_L1[++L1_nrows]  = L1_nnz;
    986:       }

    988:       /********************** L E V E L  2 ************************/
    989:       /* Merge from L1 work array to either C or to L2 work array */
    990:       if (anzi > 8) {
    991:         inputi      = worki_L1;
    992:         inputj      = workj_L1;
    993:         inputcol    = workcol;
    994:         outputi_nnz = 0;

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

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

    1016:         /************************ L E V E L  3 **********************/
    1017:         /* Merge from L2 work array to either C or to L2 work array */
    1018:         if (anzi > 64 && (L2_nrows == 8 || rowsleft == 0)) {
    1019:           inputi      = worki_L2;
    1020:           inputj      = workj_L2;
    1021:           inputcol    = workcol;
    1022:           outputi_nnz = 0;
    1023:           if (rowsleft) outputj = workj_L3;
    1024:           else          outputj = cj + ci_nnz;
    1025:           switch (L2_nrows) {
    1026:           case 1:  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:  MatMatMultSymbolic_RowMergeMacro(2); break;
    1031:           case 3:  MatMatMultSymbolic_RowMergeMacro(3); break;
    1032:           case 4:  MatMatMultSymbolic_RowMergeMacro(4); break;
    1033:           case 5:  MatMatMultSymbolic_RowMergeMacro(5); break;
    1034:           case 6:  MatMatMultSymbolic_RowMergeMacro(6); break;
    1035:           case 7:  MatMatMultSymbolic_RowMergeMacro(7); break;
    1036:           case 8:  MatMatMultSymbolic_RowMergeMacro(8); break;
    1037:           default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatMult logic error: Not merging 1-8 rows from L2 work array!");
    1038:           }
    1039:           L2_nrows    = 1;
    1040:           L2_nnz      = outputi_nnz;
    1041:           worki_L2[1] = outputi_nnz;
    1042:           /* Copy to workj_L2 */
    1043:           if (rowsleft) {
    1044:             for (k=0; k<outputi_nnz; ++k)  workj_L2[k] = outputj[k];
    1045:           }
    1046:         }
    1047:       }
    1048:     }  /* while (rowsleft) */
    1049: #undef MatMatMultSymbolic_RowMergeMacro

    1051:     /* terminate current row */
    1052:     ci_nnz += outputi_nnz;
    1053:     ci[i+1] = ci_nnz;
    1054:   }

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

    1060:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    1061:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    1062:   c          = (Mat_SeqAIJ*)(C->data);
    1063:   c->free_a  = PETSC_TRUE;
    1064:   c->free_ij = PETSC_TRUE;
    1065:   c->nonew   = 0;

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

    1069:   /* set MatInfo */
    1070:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
    1071:   if (afill < 1.0) afill = 1.0;
    1072:   c->maxnz                     = ci[am];
    1073:   c->nz                        = ci[am];
    1074:   C->info.mallocs           = ndouble;
    1075:   C->info.fill_ratio_given  = fill;
    1076:   C->info.fill_ratio_needed = afill;

    1078: #if defined(PETSC_USE_INFO)
    1079:   if (ci[am]) {
    1080:     PetscInfo3(C,"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
    1081:     PetscInfo1(C,"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
    1082:   } else {
    1083:     PetscInfo(C,"Empty matrix product\n");
    1084:   }
    1085: #endif

    1087:   /* Step 4: Free temporary work areas */
    1088:   PetscFree(workj_L1);
    1089:   PetscFree(workj_L2);
    1090:   PetscFree(workj_L3);
    1091:   return(0);
    1092: }

    1094: /* concatenate unique entries and then sort */
    1095: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(Mat A,Mat B,PetscReal fill,Mat C)
    1096: {
    1098:   Mat_SeqAIJ     *a  = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
    1099:   const PetscInt *ai = a->i,*bi=b->i,*aj=a->j,*bj=b->j;
    1100:   PetscInt       *ci,*cj,bcol;
    1101:   PetscInt       am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
    1102:   PetscReal      afill;
    1103:   PetscInt       i,j,ndouble = 0;
    1104:   PetscSegBuffer seg,segrow;
    1105:   char           *seen;

    1108:   PetscMalloc1(am+1,&ci);
    1109:   ci[0] = 0;

    1111:   /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
    1112:   PetscSegBufferCreate(sizeof(PetscInt),(PetscInt)(fill*(ai[am]+bi[bm])),&seg);
    1113:   PetscSegBufferCreate(sizeof(PetscInt),100,&segrow);
    1114:   PetscCalloc1(bn,&seen);

    1116:   /* Determine ci and cj */
    1117:   for (i=0; i<am; i++) {
    1118:     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 */
    1119:     const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
    1120:     PetscInt packlen = 0,*PETSC_RESTRICT crow;

    1122:     /* Pack segrow */
    1123:     for (j=0; j<anzi; j++) {
    1124:       PetscInt brow = acol[j],bjstart = bi[brow],bjend = bi[brow+1],k;
    1125:       for (k=bjstart; k<bjend; k++) {
    1126:         bcol = bj[k];
    1127:         if (!seen[bcol]) { /* new entry */
    1128:           PetscInt *PETSC_RESTRICT slot;
    1129:           PetscSegBufferGetInts(segrow,1,&slot);
    1130:           *slot = bcol;
    1131:           seen[bcol] = 1;
    1132:           packlen++;
    1133:         }
    1134:       }
    1135:     }

    1137:     /* Check i-th diagonal entry */
    1138:     if (C->force_diagonals && !seen[i]) {
    1139:       PetscInt *PETSC_RESTRICT slot;
    1140:       PetscSegBufferGetInts(segrow,1,&slot);
    1141:       *slot   = i;
    1142:       seen[i] = 1;
    1143:       packlen++;
    1144:     }

    1146:     PetscSegBufferGetInts(seg,packlen,&crow);
    1147:     PetscSegBufferExtractTo(segrow,crow);
    1148:     PetscSortInt(packlen,crow);
    1149:     ci[i+1] = ci[i] + packlen;
    1150:     for (j=0; j<packlen; j++) seen[crow[j]] = 0;
    1151:   }
    1152:   PetscSegBufferDestroy(&segrow);
    1153:   PetscFree(seen);

    1155:   /* Column indices are in the segmented buffer */
    1156:   PetscSegBufferExtractAlloc(seg,&cj);
    1157:   PetscSegBufferDestroy(&seg);

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

    1163:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    1164:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    1165:   c          = (Mat_SeqAIJ*)(C->data);
    1166:   c->free_a  = PETSC_TRUE;
    1167:   c->free_ij = PETSC_TRUE;
    1168:   c->nonew   = 0;

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

    1172:   /* set MatInfo */
    1173:   afill = (PetscReal)ci[am]/PetscMax(ai[am]+bi[bm],1) + 1.e-5;
    1174:   if (afill < 1.0) afill = 1.0;
    1175:   c->maxnz                  = ci[am];
    1176:   c->nz                     = ci[am];
    1177:   C->info.mallocs           = ndouble;
    1178:   C->info.fill_ratio_given  = fill;
    1179:   C->info.fill_ratio_needed = afill;

    1181: #if defined(PETSC_USE_INFO)
    1182:   if (ci[am]) {
    1183:     PetscInfo3(C,"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
    1184:     PetscInfo1(C,"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
    1185:   } else {
    1186:     PetscInfo(C,"Empty matrix product\n");
    1187:   }
    1188: #endif
    1189:   return(0);
    1190: }

    1192: PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(void *data)
    1193: {
    1194:   PetscErrorCode      ierr;
    1195:   Mat_MatMatTransMult *abt=(Mat_MatMatTransMult *)data;

    1198:   MatTransposeColoringDestroy(&abt->matcoloring);
    1199:   MatDestroy(&abt->Bt_den);
    1200:   MatDestroy(&abt->ABt_den);
    1201:   PetscFree(abt);
    1202:   return(0);
    1203: }

    1205: PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C)
    1206: {
    1207:   PetscErrorCode      ierr;
    1208:   Mat                 Bt;
    1209:   PetscInt            *bti,*btj;
    1210:   Mat_MatMatTransMult *abt;
    1211:   Mat_Product         *product = C->product;
    1212:   char                *alg;

    1215:   if (!product) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing product struct");
    1216:   if (product->data) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Extra product struct not empty");

    1218:   /* create symbolic Bt */
    1219:   MatGetSymbolicTranspose_SeqAIJ(B,&bti,&btj);
    1220:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,B->cmap->n,B->rmap->n,bti,btj,NULL,&Bt);
    1221:   MatSetBlockSizes(Bt,PetscAbs(A->cmap->bs),PetscAbs(B->cmap->bs));
    1222:   MatSetType(Bt,((PetscObject)A)->type_name);

    1224:   /* get symbolic C=A*Bt */
    1225:   PetscStrallocpy(product->alg,&alg);
    1226:   MatProductSetAlgorithm(C,"sorted"); /* set algorithm for C = A*Bt */
    1227:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,Bt,fill,C);
    1228:   MatProductSetAlgorithm(C,alg); /* resume original algorithm for ABt product */
    1229:   PetscFree(alg);

    1231:   /* create a supporting struct for reuse intermidiate dense matrices with matcoloring */
    1232:   PetscNew(&abt);

    1234:   product->data    = abt;
    1235:   product->destroy = MatDestroy_SeqAIJ_MatMatMultTrans;

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

    1239:   abt->usecoloring = PETSC_FALSE;
    1240:   PetscStrcmp(product->alg,"color",&abt->usecoloring);
    1241:   if (abt->usecoloring) {
    1242:     /* Create MatTransposeColoring from symbolic C=A*B^T */
    1243:     MatTransposeColoring matcoloring;
    1244:     MatColoring          coloring;
    1245:     ISColoring           iscoloring;
    1246:     Mat                  Bt_dense,C_dense;

    1248:     /* inode causes memory problem */
    1249:     MatSetOption(C,MAT_USE_INODES,PETSC_FALSE);

    1251:     MatColoringCreate(C,&coloring);
    1252:     MatColoringSetDistance(coloring,2);
    1253:     MatColoringSetType(coloring,MATCOLORINGSL);
    1254:     MatColoringSetFromOptions(coloring);
    1255:     MatColoringApply(coloring,&iscoloring);
    1256:     MatColoringDestroy(&coloring);
    1257:     MatTransposeColoringCreate(C,iscoloring,&matcoloring);

    1259:     abt->matcoloring = matcoloring;

    1261:     ISColoringDestroy(&iscoloring);

    1263:     /* Create Bt_dense and C_dense = A*Bt_dense */
    1264:     MatCreate(PETSC_COMM_SELF,&Bt_dense);
    1265:     MatSetSizes(Bt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);
    1266:     MatSetType(Bt_dense,MATSEQDENSE);
    1267:     MatSeqDenseSetPreallocation(Bt_dense,NULL);

    1269:     Bt_dense->assembled = PETSC_TRUE;
    1270:     abt->Bt_den         = Bt_dense;

    1272:     MatCreate(PETSC_COMM_SELF,&C_dense);
    1273:     MatSetSizes(C_dense,A->rmap->n,matcoloring->ncolors,A->rmap->n,matcoloring->ncolors);
    1274:     MatSetType(C_dense,MATSEQDENSE);
    1275:     MatSeqDenseSetPreallocation(C_dense,NULL);

    1277:     Bt_dense->assembled = PETSC_TRUE;
    1278:     abt->ABt_den  = C_dense;

    1280: #if defined(PETSC_USE_INFO)
    1281:     {
    1282:       Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
    1283:       PetscInfo7(C,"Use coloring of C=A*B^T; B^T: %D %D, Bt_dense: %D,%D; Cnz %D / (cm*ncolors %D) = %g\n",B->cmap->n,B->rmap->n,Bt_dense->rmap->n,Bt_dense->cmap->n,c->nz,A->rmap->n*matcoloring->ncolors,(PetscReal)(c->nz)/(A->rmap->n*matcoloring->ncolors));
    1284:     }
    1285: #endif
    1286:   }
    1287:   /* clean up */
    1288:   MatDestroy(&Bt);
    1289:   MatRestoreSymbolicTranspose_SeqAIJ(B,&bti,&btj);
    1290:   return(0);
    1291: }

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

    1305:   if (!product) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing product struct");
    1306:   abt = (Mat_MatMatTransMult *)product->data;
    1307:   if (!abt) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing product struct");
    1308:   /* clear old values in C */
    1309:   if (!c->a) {
    1310:     PetscCalloc1(ci[cm]+1,&ca);
    1311:     c->a      = ca;
    1312:     c->free_a = PETSC_TRUE;
    1313:   } else {
    1314:     ca =  c->a;
    1315:     PetscArrayzero(ca,ci[cm]+1);
    1316:   }

    1318:   if (abt->usecoloring) {
    1319:     MatTransposeColoring matcoloring = abt->matcoloring;
    1320:     Mat                  Bt_dense,C_dense = abt->ABt_den;

    1322:     /* Get Bt_dense by Apply MatTransposeColoring to B */
    1323:     Bt_dense = abt->Bt_den;
    1324:     MatTransColoringApplySpToDen(matcoloring,B,Bt_dense);

    1326:     /* C_dense = A*Bt_dense */
    1327:     MatMatMultNumeric_SeqAIJ_SeqDense(A,Bt_dense,C_dense);

    1329:     /* Recover C from C_dense */
    1330:     MatTransColoringApplyDenToSp(matcoloring,C_dense,C);
    1331:     return(0);
    1332:   }

    1334:   for (i=0; i<cm; i++) {
    1335:     anzi = ai[i+1] - ai[i];
    1336:     acol = aj + ai[i];
    1337:     aval = aa + ai[i];
    1338:     cnzi = ci[i+1] - ci[i];
    1339:     ccol = cj + ci[i];
    1340:     cval = ca + ci[i];
    1341:     for (j=0; j<cnzi; j++) {
    1342:       brow = ccol[j];
    1343:       bnzj = bi[brow+1] - bi[brow];
    1344:       bcol = bj + bi[brow];
    1345:       bval = ba + bi[brow];

    1347:       /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */
    1348:       nexta = 0; nextb = 0;
    1349:       while (nexta<anzi && nextb<bnzj) {
    1350:         while (nexta < anzi && acol[nexta] < bcol[nextb]) nexta++;
    1351:         if (nexta == anzi) break;
    1352:         while (nextb < bnzj && acol[nexta] > bcol[nextb]) nextb++;
    1353:         if (nextb == bnzj) break;
    1354:         if (acol[nexta] == bcol[nextb]) {
    1355:           cval[j] += aval[nexta]*bval[nextb];
    1356:           nexta++; nextb++;
    1357:           flops += 2;
    1358:         }
    1359:       }
    1360:     }
    1361:   }
    1362:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    1363:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
    1364:   PetscLogFlops(flops);
    1365:   return(0);
    1366: }

    1368: PetscErrorCode MatDestroy_SeqAIJ_MatTransMatMult(void *data)
    1369: {
    1370:   PetscErrorCode      ierr;
    1371:   Mat_MatTransMatMult *atb = (Mat_MatTransMatMult*)data;

    1374:   MatDestroy(&atb->At);
    1375:   if (atb->destroy) {
    1376:     (*atb->destroy)(atb->data);
    1377:   }
    1378:   PetscFree(atb);
    1379:   return(0);
    1380: }

    1382: PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C)
    1383: {
    1385:   Mat            At = NULL;
    1386:   PetscInt       *ati,*atj;
    1387:   Mat_Product    *product = C->product;
    1388:   PetscBool      flg,def,square;

    1391:   MatCheckProduct(C,4);
    1392:   square = (PetscBool)(A == B && A->symmetric && A->symmetric_set);
    1393:   /* outerproduct */
    1394:   PetscStrcmp(product->alg,"outerproduct",&flg);
    1395:   if (flg) {
    1396:     /* create symbolic At */
    1397:     if (!square) {
    1398:       MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);
    1399:       MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,A->cmap->n,A->rmap->n,ati,atj,NULL,&At);
    1400:       MatSetBlockSizes(At,PetscAbs(A->cmap->bs),PetscAbs(B->cmap->bs));
    1401:       MatSetType(At,((PetscObject)A)->type_name);
    1402:     }
    1403:     /* get symbolic C=At*B */
    1404:     MatProductSetAlgorithm(C,"sorted");
    1405:     MatMatMultSymbolic_SeqAIJ_SeqAIJ(square ? A : At,B,fill,C);

    1407:     /* clean up */
    1408:     if (!square) {
    1409:       MatDestroy(&At);
    1410:       MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);
    1411:     }

    1413:     C->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ; /* outerproduct */
    1414:     MatProductSetAlgorithm(C,"outerproduct");
    1415:     return(0);
    1416:   }

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

    1424:     if (product->data) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Extra product struct not empty");
    1425:     PetscNew(&atb);
    1426:     if (!square) {
    1427:       MatTranspose_SeqAIJ(A,MAT_INITIAL_MATRIX,&At);
    1428:     }
    1429:     MatProductSetAlgorithm(C,"sorted");
    1430:     MatMatMultSymbolic_SeqAIJ_SeqAIJ(square ? A : At,B,fill,C);
    1431:     MatProductSetAlgorithm(C,"at*b");
    1432:     product->data    = atb;
    1433:     product->destroy = MatDestroy_SeqAIJ_MatTransMatMult;
    1434:     atb->At          = At;
    1435:     atb->updateAt    = PETSC_FALSE; /* because At is computed here */

    1437:     C->ops->mattransposemultnumeric = NULL; /* see MatProductNumeric_AtB_SeqAIJ_SeqAIJ */
    1438:     return(0);
    1439:   }

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

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

    1454:   if (!c->a) {
    1455:     PetscCalloc1(ci[cm]+1,&ca);

    1457:     c->a      = ca;
    1458:     c->free_a = PETSC_TRUE;
    1459:   } else {
    1460:     ca   = c->a;
    1461:     PetscArrayzero(ca,ci[cm]);
    1462:   }

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

    1487:   /* Assemble the final matrix and clean up */
    1488:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    1489:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
    1490:   PetscLogFlops(flops);
    1491:   return(0);
    1492: }

    1494: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
    1495: {

    1499:   MatMatMultSymbolic_SeqDense_SeqDense(A,B,0.0,C);
    1500:   C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
    1501:   return(0);
    1502: }

    1504: PETSC_INTERN PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A,Mat B,Mat C,const PetscBool add)
    1505: {
    1506:   Mat_SeqAIJ        *a=(Mat_SeqAIJ*)A->data;
    1507:   Mat_SeqDense      *bd=(Mat_SeqDense*)B->data;
    1508:   Mat_SeqDense      *cd=(Mat_SeqDense*)C->data;
    1509:   PetscErrorCode    ierr;
    1510:   PetscScalar       *c,r1,r2,r3,r4,*c1,*c2,*c3,*c4;
    1511:   const PetscScalar *aa,*b,*b1,*b2,*b3,*b4,*av;
    1512:   const PetscInt    *aj;
    1513:   PetscInt          cm=C->rmap->n,cn=B->cmap->n,bm=bd->lda,am=A->rmap->n;
    1514:   PetscInt          clda=cd->lda;
    1515:   PetscInt          am4=4*clda,bm4=4*bm,col,i,j,n;

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

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

    1627: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
    1628: {

    1632:   if (B->rmap->n != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in A %D not equal rows in B %D\n",A->cmap->n,B->rmap->n);
    1633:   if (A->rmap->n != C->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows in C %D not equal rows in A %D\n",C->rmap->n,A->rmap->n);
    1634:   if (B->cmap->n != C->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in B %D not equal columns in C %D\n",B->cmap->n,C->cmap->n);

    1636:   MatMatMultNumericAdd_SeqAIJ_SeqDense(A,B,C,PETSC_FALSE);
    1637:   return(0);
    1638: }

    1640: /* ------------------------------------------------------- */
    1641: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_AB(Mat C)
    1642: {
    1644:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqAIJ_SeqDense;
    1645:   C->ops->productsymbolic = MatProductSymbolic_AB;
    1646:   return(0);
    1647: }

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

    1651: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_AtB(Mat C)
    1652: {
    1654:   C->ops->transposematmultsymbolic = MatTMatTMultSymbolic_SeqAIJ_SeqDense;
    1655:   C->ops->productsymbolic          = MatProductSymbolic_AtB;
    1656:   return(0);
    1657: }

    1659: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_ABt(Mat C)
    1660: {
    1662:   C->ops->mattransposemultsymbolic = MatTMatTMultSymbolic_SeqAIJ_SeqDense;
    1663:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
    1664:   return(0);
    1665: }

    1667: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense(Mat C)
    1668: {
    1670:   Mat_Product    *product = C->product;

    1673:   switch (product->type) {
    1674:   case MATPRODUCT_AB:
    1675:     MatProductSetFromOptions_SeqAIJ_SeqDense_AB(C);
    1676:     break;
    1677:   case MATPRODUCT_AtB:
    1678:     MatProductSetFromOptions_SeqAIJ_SeqDense_AtB(C);
    1679:     break;
    1680:   case MATPRODUCT_ABt:
    1681:     MatProductSetFromOptions_SeqAIJ_SeqDense_ABt(C);
    1682:     break;
    1683:   default:
    1684:     break;
    1685:   }
    1686:   return(0);
    1687: }
    1688: /* ------------------------------------------------------- */
    1689: static PetscErrorCode MatProductSetFromOptions_SeqXBAIJ_SeqDense_AB(Mat C)
    1690: {
    1692:   Mat_Product    *product = C->product;
    1693:   Mat            A = product->A;
    1694:   PetscBool      baij;

    1697:   PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&baij);
    1698:   if (!baij) { /* A is seqsbaij */
    1699:     PetscBool sbaij;
    1700:     PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&sbaij);
    1701:     if (!sbaij) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Mat must be either seqbaij or seqsbaij format");

    1703:     C->ops->matmultsymbolic = MatMatMultSymbolic_SeqSBAIJ_SeqDense;
    1704:   } else { /* A is seqbaij */
    1705:     C->ops->matmultsymbolic = MatMatMultSymbolic_SeqBAIJ_SeqDense;
    1706:   }

    1708:   C->ops->productsymbolic = MatProductSymbolic_AB;
    1709:   return(0);
    1710: }

    1712: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqXBAIJ_SeqDense(Mat C)
    1713: {
    1715:   Mat_Product    *product = C->product;

    1718:   MatCheckProduct(C,1);
    1719:   if (!product->A) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing A");
    1720:   if (product->type == MATPRODUCT_AB || (product->type == MATPRODUCT_AtB && product->A->symmetric)) {
    1721:     MatProductSetFromOptions_SeqXBAIJ_SeqDense_AB(C);
    1722:   }
    1723:   return(0);
    1724: }

    1726: /* ------------------------------------------------------- */
    1727: static PetscErrorCode MatProductSetFromOptions_SeqDense_SeqAIJ_AB(Mat C)
    1728: {
    1730:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqAIJ;
    1731:   C->ops->productsymbolic = MatProductSymbolic_AB;
    1732:   return(0);
    1733: }

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

    1741:   if (product->type == MATPRODUCT_AB) {
    1742:     MatProductSetFromOptions_SeqDense_SeqAIJ_AB(C);
    1743:   }
    1744:   return(0);
    1745: }
    1746: /* ------------------------------------------------------- */

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

    1759:   btval_den=btdense->v;
    1760:   PetscArrayzero(btval_den,m*n);
    1761:   for (k=0; k<ncolors; k++) {
    1762:     ncolumns = coloring->ncolumns[k];
    1763:     for (l=0; l<ncolumns; l++) { /* insert a row of B to a column of Btdense */
    1764:       col   = *(columns + colorforcol[k] + l);
    1765:       btcol = bj + bi[col];
    1766:       btval = ba + bi[col];
    1767:       anz   = bi[col+1] - bi[col];
    1768:       for (j=0; j<anz; j++) {
    1769:         brow            = btcol[j];
    1770:         btval_den[brow] = btval[j];
    1771:       }
    1772:     }
    1773:     btval_den += m;
    1774:   }
    1775:   return(0);
    1776: }

    1778: PetscErrorCode MatTransColoringApplyDenToSp_SeqAIJ(MatTransposeColoring matcoloring,Mat Cden,Mat Csp)
    1779: {
    1780:   PetscErrorCode    ierr;
    1781:   Mat_SeqAIJ        *csp = (Mat_SeqAIJ*)Csp->data;
    1782:   const PetscScalar *ca_den,*ca_den_ptr;
    1783:   PetscScalar       *ca=csp->a;
    1784:   PetscInt          k,l,m=Cden->rmap->n,ncolors=matcoloring->ncolors;
    1785:   PetscInt          brows=matcoloring->brows,*den2sp=matcoloring->den2sp;
    1786:   PetscInt          nrows,*row,*idx;
    1787:   PetscInt          *rows=matcoloring->rows,*colorforrow=matcoloring->colorforrow;

    1790:   MatDenseGetArrayRead(Cden,&ca_den);

    1792:   if (brows > 0) {
    1793:     PetscInt *lstart,row_end,row_start;
    1794:     lstart = matcoloring->lstart;
    1795:     PetscArrayzero(lstart,ncolors);

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

    1831:   MatDenseRestoreArrayRead(Cden,&ca_den);
    1832: #if defined(PETSC_USE_INFO)
    1833:   if (matcoloring->brows > 0) {
    1834:     PetscInfo1(Csp,"Loop over %D row blocks for den2sp\n",brows);
    1835:   } else {
    1836:     PetscInfo(Csp,"Loop over colors/columns of Cden, inefficient for large sparse matrix product \n");
    1837:   }
    1838: #endif
    1839:   return(0);
    1840: }

    1842: PetscErrorCode MatTransposeColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatTransposeColoring c)
    1843: {
    1845:   PetscInt       i,n,nrows,Nbs,j,k,m,ncols,col,cm;
    1846:   const PetscInt *is,*ci,*cj,*row_idx;
    1847:   PetscInt       nis = iscoloring->n,*rowhit,bs = 1;
    1848:   IS             *isa;
    1849:   Mat_SeqAIJ     *csp = (Mat_SeqAIJ*)mat->data;
    1850:   PetscInt       *colorforrow,*rows,*rows_i,*idxhit,*spidx,*den2sp,*den2sp_i;
    1851:   PetscInt       *colorforcol,*columns,*columns_i,brows;
    1852:   PetscBool      flg;

    1855:   ISColoringGetIS(iscoloring,PETSC_USE_POINTER,PETSC_IGNORE,&isa);

    1857:   /* bs >1 is not being tested yet! */
    1858:   Nbs       = mat->cmap->N/bs;
    1859:   c->M      = mat->rmap->N/bs;  /* set total rows, columns and local rows */
    1860:   c->N      = Nbs;
    1861:   c->m      = c->M;
    1862:   c->rstart = 0;
    1863:   c->brows  = 100;

    1865:   c->ncolors = nis;
    1866:   PetscMalloc3(nis,&c->ncolumns,nis,&c->nrows,nis+1,&colorforrow);
    1867:   PetscMalloc1(csp->nz+1,&rows);
    1868:   PetscMalloc1(csp->nz+1,&den2sp);

    1870:   brows = c->brows;
    1871:   PetscOptionsGetInt(NULL,NULL,"-matden2sp_brows",&brows,&flg);
    1872:   if (flg) c->brows = brows;
    1873:   if (brows > 0) {
    1874:     PetscMalloc1(nis+1,&c->lstart);
    1875:   }

    1877:   colorforrow[0] = 0;
    1878:   rows_i         = rows;
    1879:   den2sp_i       = den2sp;

    1881:   PetscMalloc1(nis+1,&colorforcol);
    1882:   PetscMalloc1(Nbs+1,&columns);

    1884:   colorforcol[0] = 0;
    1885:   columns_i      = columns;

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

    1890:   cm   = c->m;
    1891:   PetscMalloc1(cm+1,&rowhit);
    1892:   PetscMalloc1(cm+1,&idxhit);
    1893:   for (i=0; i<nis; i++) { /* loop over color */
    1894:     ISGetLocalSize(isa[i],&n);
    1895:     ISGetIndices(isa[i],&is);

    1897:     c->ncolumns[i] = n;
    1898:     if (n) {
    1899:       PetscArraycpy(columns_i,is,n);
    1900:     }
    1901:     colorforcol[i+1] = colorforcol[i] + n;
    1902:     columns_i       += n;

    1904:     /* fast, crude version requires O(N*N) work */
    1905:     PetscArrayzero(rowhit,cm);

    1907:     for (j=0; j<n; j++) { /* loop over columns*/
    1908:       col     = is[j];
    1909:       row_idx = cj + ci[col];
    1910:       m       = ci[col+1] - ci[col];
    1911:       for (k=0; k<m; k++) { /* loop over columns marking them in rowhit */
    1912:         idxhit[*row_idx]   = spidx[ci[col] + k];
    1913:         rowhit[*row_idx++] = col + 1;
    1914:       }
    1915:     }
    1916:     /* count the number of hits */
    1917:     nrows = 0;
    1918:     for (j=0; j<cm; j++) {
    1919:       if (rowhit[j]) nrows++;
    1920:     }
    1921:     c->nrows[i]      = nrows;
    1922:     colorforrow[i+1] = colorforrow[i] + nrows;

    1924:     nrows = 0;
    1925:     for (j=0; j<cm; j++) { /* loop over rows */
    1926:       if (rowhit[j]) {
    1927:         rows_i[nrows]   = j;
    1928:         den2sp_i[nrows] = idxhit[j];
    1929:         nrows++;
    1930:       }
    1931:     }
    1932:     den2sp_i += nrows;

    1934:     ISRestoreIndices(isa[i],&is);
    1935:     rows_i += nrows;
    1936:   }
    1937:   MatRestoreColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);
    1938:   PetscFree(rowhit);
    1939:   ISColoringRestoreIS(iscoloring,PETSC_USE_POINTER,&isa);
    1940:   if (csp->nz != colorforrow[nis]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"csp->nz %d != colorforrow[nis] %d",csp->nz,colorforrow[nis]);

    1942:   c->colorforrow = colorforrow;
    1943:   c->rows        = rows;
    1944:   c->den2sp      = den2sp;
    1945:   c->colorforcol = colorforcol;
    1946:   c->columns     = columns;

    1948:   PetscFree(idxhit);
    1949:   return(0);
    1950: }

    1952: /* --------------------------------------------------------------- */
    1953: static PetscErrorCode MatProductNumeric_AtB_SeqAIJ_SeqAIJ(Mat C)
    1954: {
    1956:   Mat_Product    *product = C->product;
    1957:   Mat            A=product->A,B=product->B;

    1960:   if (C->ops->mattransposemultnumeric) {
    1961:     /* Alg: "outerproduct" */
    1962:     (*C->ops->mattransposemultnumeric)(A,B,C);
    1963:   } else {
    1964:     /* Alg: "matmatmult" -- C = At*B */
    1965:     Mat_MatTransMatMult *atb = (Mat_MatTransMatMult *)product->data;
    1966:     Mat                 At;

    1968:     if (!atb) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing product struct");
    1969:     At = atb->At;
    1970:     if (atb->updateAt && At) { /* At is computed in MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ() */
    1971:       MatTranspose_SeqAIJ(A,MAT_REUSE_MATRIX,&At);
    1972:     }
    1973:     MatMatMultNumeric_SeqAIJ_SeqAIJ(At ? At : A,B,C);
    1974:     atb->updateAt = PETSC_TRUE;
    1975:   }
    1976:   return(0);
    1977: }

    1979: static PetscErrorCode MatProductSymbolic_AtB_SeqAIJ_SeqAIJ(Mat C)
    1980: {
    1982:   Mat_Product    *product = C->product;
    1983:   Mat            A=product->A,B=product->B;
    1984:   PetscReal      fill=product->fill;

    1987:   MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);

    1989:   C->ops->productnumeric = MatProductNumeric_AtB_SeqAIJ_SeqAIJ;
    1990:   return(0);
    1991: }

    1993: /* --------------------------------------------------------------- */
    1994: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_AB(Mat C)
    1995: {
    1997:   Mat_Product    *product = C->product;
    1998:   PetscInt       alg = 0; /* default algorithm */
    1999:   PetscBool      flg = PETSC_FALSE;
    2000: #if !defined(PETSC_HAVE_HYPRE)
    2001:   const char     *algTypes[7] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","rowmerge"};
    2002:   PetscInt       nalg = 7;
    2003: #else
    2004:   const char     *algTypes[8] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","rowmerge","hypre"};
    2005:   PetscInt       nalg = 8;
    2006: #endif

    2009:   /* Set default algorithm */
    2010:   PetscStrcmp(C->product->alg,"default",&flg);
    2011:   if (flg) {
    2012:     MatProductSetAlgorithm(C,(MatProductAlgorithm)algTypes[alg]);
    2013:   }

    2015:   /* Get runtime option */
    2016:   if (product->api_user) {
    2017:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatMatMult","Mat");
    2018:     PetscOptionsEList("-matmatmult_via","Algorithmic approach","MatMatMult",algTypes,nalg,algTypes[0],&alg,&flg);
    2019:     PetscOptionsEnd();
    2020:   } else {
    2021:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatProduct_AB","Mat");
    2022:     PetscOptionsEList("-matproduct_ab_via","Algorithmic approach","MatProduct_AB",algTypes,nalg,algTypes[0],&alg,&flg);
    2023:     PetscOptionsEnd();
    2024:   }
    2025:   if (flg) {
    2026:     MatProductSetAlgorithm(C,(MatProductAlgorithm)algTypes[alg]);
    2027:   }

    2029:   C->ops->productsymbolic = MatProductSymbolic_AB;
    2030:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqAIJ_SeqAIJ;
    2031:   return(0);
    2032: }

    2034: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_AtB(Mat C)
    2035: {
    2037:   Mat_Product    *product = C->product;
    2038:   PetscInt       alg = 0; /* default algorithm */
    2039:   PetscBool      flg = PETSC_FALSE;
    2040:   const char     *algTypes[3] = {"default","at*b","outerproduct"};
    2041:   PetscInt       nalg = 3;

    2044:   /* Get runtime option */
    2045:   if (product->api_user) {
    2046:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatTransposeMatMult","Mat");
    2047:     PetscOptionsEList("-mattransposematmult_via","Algorithmic approach","MatTransposeMatMult",algTypes,nalg,algTypes[alg],&alg,&flg);
    2048:     PetscOptionsEnd();
    2049:   } else {
    2050:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatProduct_AtB","Mat");
    2051:     PetscOptionsEList("-matproduct_atb_via","Algorithmic approach","MatProduct_AtB",algTypes,nalg,algTypes[alg],&alg,&flg);
    2052:     PetscOptionsEnd();
    2053:   }
    2054:   if (flg) {
    2055:     MatProductSetAlgorithm(C,(MatProductAlgorithm)algTypes[alg]);
    2056:   }

    2058:   C->ops->productsymbolic = MatProductSymbolic_AtB_SeqAIJ_SeqAIJ;
    2059:   return(0);
    2060: }

    2062: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_ABt(Mat C)
    2063: {
    2065:   Mat_Product    *product = C->product;
    2066:   PetscInt       alg = 0; /* default algorithm */
    2067:   PetscBool      flg = PETSC_FALSE;
    2068:   const char     *algTypes[2] = {"default","color"};
    2069:   PetscInt       nalg = 2;

    2072:   /* Set default algorithm */
    2073:   PetscStrcmp(C->product->alg,"default",&flg);
    2074:   if (!flg) {
    2075:     alg = 1;
    2076:     MatProductSetAlgorithm(C,(MatProductAlgorithm)algTypes[alg]);
    2077:   }

    2079:   /* Get runtime option */
    2080:   if (product->api_user) {
    2081:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatMatTransposeMult","Mat");
    2082:     PetscOptionsEList("-matmattransmult_via","Algorithmic approach","MatMatTransposeMult",algTypes,nalg,algTypes[alg],&alg,&flg);
    2083:     PetscOptionsEnd();
    2084:   } else {
    2085:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatProduct_ABt","Mat");
    2086:     PetscOptionsEList("-matproduct_abt_via","Algorithmic approach","MatProduct_ABt",algTypes,nalg,algTypes[alg],&alg,&flg);
    2087:     PetscOptionsEnd();
    2088:   }
    2089:   if (flg) {
    2090:     MatProductSetAlgorithm(C,(MatProductAlgorithm)algTypes[alg]);
    2091:   }

    2093:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ;
    2094:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
    2095:   return(0);
    2096: }

    2098: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_PtAP(Mat C)
    2099: {
    2101:   Mat_Product    *product = C->product;
    2102:   PetscBool      flg = PETSC_FALSE;
    2103:   PetscInt       alg = 0; /* default algorithm -- alg=1 should be default!!! */
    2104: #if !defined(PETSC_HAVE_HYPRE)
    2105:   const char     *algTypes[2] = {"scalable","rap"};
    2106:   PetscInt       nalg = 2;
    2107: #else
    2108:   const char     *algTypes[3] = {"scalable","rap","hypre"};
    2109:   PetscInt       nalg = 3;
    2110: #endif

    2113:   /* Set default algorithm */
    2114:   PetscStrcmp(product->alg,"default",&flg);
    2115:   if (flg) {
    2116:     MatProductSetAlgorithm(C,(MatProductAlgorithm)algTypes[alg]);
    2117:   }

    2119:   /* Get runtime option */
    2120:   if (product->api_user) {
    2121:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatPtAP","Mat");
    2122:     PetscOptionsEList("-matptap_via","Algorithmic approach","MatPtAP",algTypes,nalg,algTypes[0],&alg,&flg);
    2123:     PetscOptionsEnd();
    2124:   } else {
    2125:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatProduct_PtAP","Mat");
    2126:     PetscOptionsEList("-matproduct_ptap_via","Algorithmic approach","MatProduct_PtAP",algTypes,nalg,algTypes[0],&alg,&flg);
    2127:     PetscOptionsEnd();
    2128:   }
    2129:   if (flg) {
    2130:     MatProductSetAlgorithm(C,(MatProductAlgorithm)algTypes[alg]);
    2131:   }

    2133:   C->ops->productsymbolic = MatProductSymbolic_PtAP_SeqAIJ_SeqAIJ;
    2134:   return(0);
    2135: }

    2137: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_RARt(Mat C)
    2138: {
    2140:   Mat_Product    *product = C->product;
    2141:   PetscBool      flg = PETSC_FALSE;
    2142:   PetscInt       alg = 0; /* default algorithm */
    2143:   const char     *algTypes[3] = {"r*a*rt","r*art","coloring_rart"};
    2144:   PetscInt        nalg = 3;

    2147:   /* Set default algorithm */
    2148:   PetscStrcmp(product->alg,"default",&flg);
    2149:   if (flg) {
    2150:     MatProductSetAlgorithm(C,(MatProductAlgorithm)algTypes[alg]);
    2151:   }

    2153:   /* Get runtime option */
    2154:   if (product->api_user) {
    2155:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatRARt","Mat");
    2156:     PetscOptionsEList("-matrart_via","Algorithmic approach","MatRARt",algTypes,nalg,algTypes[0],&alg,&flg);
    2157:     PetscOptionsEnd();
    2158:   } else {
    2159:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatProduct_RARt","Mat");
    2160:     PetscOptionsEList("-matproduct_rart_via","Algorithmic approach","MatProduct_RARt",algTypes,nalg,algTypes[0],&alg,&flg);
    2161:     PetscOptionsEnd();
    2162:   }
    2163:   if (flg) {
    2164:     MatProductSetAlgorithm(C,(MatProductAlgorithm)algTypes[alg]);
    2165:   }

    2167:   C->ops->productsymbolic = MatProductSymbolic_RARt_SeqAIJ_SeqAIJ;
    2168:   return(0);
    2169: }

    2171: /* ABC = A*B*C = A*(B*C); ABC's algorithm must be chosen from AB's algorithm */
    2172: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_ABC(Mat C)
    2173: {
    2175:   Mat_Product    *product = C->product;
    2176:   PetscInt       alg = 0; /* default algorithm */
    2177:   PetscBool      flg = PETSC_FALSE;
    2178:   const char     *algTypes[7] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","rowmerge"};
    2179:   PetscInt       nalg = 7;

    2182:   /* Set default algorithm */
    2183:   PetscStrcmp(product->alg,"default",&flg);
    2184:   if (flg) {
    2185:     MatProductSetAlgorithm(C,(MatProductAlgorithm)algTypes[alg]);
    2186:   }

    2188:   /* Get runtime option */
    2189:   if (product->api_user) {
    2190:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatMatMatMult","Mat");
    2191:     PetscOptionsEList("-matmatmatmult_via","Algorithmic approach","MatMatMatMult",algTypes,nalg,algTypes[alg],&alg,&flg);
    2192:     PetscOptionsEnd();
    2193:   } else {
    2194:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatProduct_ABC","Mat");
    2195:     PetscOptionsEList("-matproduct_abc_via","Algorithmic approach","MatProduct_ABC",algTypes,nalg,algTypes[alg],&alg,&flg);
    2196:     PetscOptionsEnd();
    2197:   }
    2198:   if (flg) {
    2199:     MatProductSetAlgorithm(C,(MatProductAlgorithm)algTypes[alg]);
    2200:   }

    2202:   C->ops->matmatmultsymbolic = MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ;
    2203:   C->ops->productsymbolic    = MatProductSymbolic_ABC;
    2204:   return(0);
    2205: }

    2207: PetscErrorCode MatProductSetFromOptions_SeqAIJ(Mat C)
    2208: {
    2210:   Mat_Product    *product = C->product;

    2213:   switch (product->type) {
    2214:   case MATPRODUCT_AB:
    2215:     MatProductSetFromOptions_SeqAIJ_AB(C);
    2216:     break;
    2217:   case MATPRODUCT_AtB:
    2218:     MatProductSetFromOptions_SeqAIJ_AtB(C);
    2219:     break;
    2220:   case MATPRODUCT_ABt:
    2221:     MatProductSetFromOptions_SeqAIJ_ABt(C);
    2222:     break;
    2223:   case MATPRODUCT_PtAP:
    2224:     MatProductSetFromOptions_SeqAIJ_PtAP(C);
    2225:     break;
    2226:   case MATPRODUCT_RARt:
    2227:     MatProductSetFromOptions_SeqAIJ_RARt(C);
    2228:     break;
    2229:   case MATPRODUCT_ABC:
    2230:     MatProductSetFromOptions_SeqAIJ_ABC(C);
    2231:     break;
    2232:   default:
    2233:     break;
    2234:   }
    2235:   return(0);
    2236: }