Actual source code: matmatmult.c

petsc-3.11.4 2019-09-28
Report Typos and Errors

  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> /*I "petscmat.h" I*/
  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>


 14:  PETSC_INTERN PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
 15:  {

 19:    if (scall == MAT_INITIAL_MATRIX) {
 20:      PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
 21:      MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);
 22:      PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
 23:    }

 25:    PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
 26:    MatMatMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);
 27:    PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
 28:    return(0);
 29:  }

 31:  PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
 32:  {

 36:    if (C->ops->matmultnumeric) {
 37:      (*C->ops->matmultnumeric)(A,B,C);
 38:    } else {
 39:      MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(A,B,C);
 40:    }
 41:    return(0);
 42:  }

 44:  PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
 45:  {
 47:  #if !defined(PETSC_HAVE_HYPRE)
 48:    const char     *algTypes[8] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","combined","rowmerge"};
 49:    PetscInt       nalg = 8;
 50:  #else
 51:    const char     *algTypes[9] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","combined","rowmerge","hypre"};
 52:    PetscInt       nalg = 9;
 53:  #endif
 54:    PetscInt       alg = 0; /* set default algorithm */

 57:    PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatMatMult","Mat");
 58:    PetscOptionsEList("-matmatmult_via","Algorithmic approach","MatMatMult",algTypes,nalg,algTypes[0],&alg,NULL);
 59:    PetscOptionsEnd();
 60:    switch (alg) {
 61:    case 1:
 62:      MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A,B,fill,C);
 63:      break;
 64:    case 2:
 65:      MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(A,B,fill,C);
 66:      break;
 67:    case 3:
 68:      MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(A,B,fill,C);
 69:      break;
 70:    case 4:
 71:      MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(A,B,fill,C);
 72:      break;
 73:    case 5:
 74:      MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(A,B,fill,C);
 75:      break;
 76:    case 6:
 77:      MatMatMult_SeqAIJ_SeqAIJ_Combined(A,B,fill,C);
 78:      break;
 79:    case 7:
 80:      MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(A,B,fill,C);
 81:      break;
 82:  #if defined(PETSC_HAVE_HYPRE)
 83:    case 8:
 84:      MatMatMultSymbolic_AIJ_AIJ_wHYPRE(A,B,fill,C);
 85:      break;
 86:  #endif
 87:    default:
 88:      MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(A,B,fill,C);
 89:      break;
 90:    }
 91:    return(0);
 92:  }

 94:  PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(Mat A,Mat B,PetscReal fill,Mat *C)
 95:  {
 96:    PetscErrorCode     ierr;
 97:    Mat_SeqAIJ         *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
 98:    PetscInt           *ai=a->i,*bi=b->i,*ci,*cj;
 99:    PetscInt           am =A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
100:    PetscReal          afill;
101:    PetscInt           i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;
102:    PetscTable         ta;
103:    PetscBT            lnkbt;
104:    PetscFreeSpaceList free_space=NULL,current_space=NULL;

107:    /* Get ci and cj */
108:    /*---------------*/
109:    /* Allocate ci array, arrays for fill computation and */
110:    /* free space for accumulating nonzero column info */
111:    PetscMalloc1(am+2,&ci);
112:    ci[0] = 0;

114:    /* create and initialize a linked list */
115:    PetscTableCreate(bn,bn,&ta);
116:    MatRowMergeMax_SeqAIJ(b,bm,ta);
117:    PetscTableGetCount(ta,&Crmax);
118:    PetscTableDestroy(&ta);

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

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

125:    current_space = free_space;

127:    /* Determine ci and cj */
128:    for (i=0; i<am; i++) {
129:      anzi = ai[i+1] - ai[i];
130:      aj   = a->j + ai[i];
131:      for (j=0; j<anzi; j++) {
132:        brow = aj[j];
133:        bnzj = bi[brow+1] - bi[brow];
134:        bj   = b->j + bi[brow];
135:        /* add non-zero cols of B into the sorted linked list lnk */
136:        PetscLLCondensedAddSorted(bnzj,bj,lnk,lnkbt);
137:      }
138:      cnzi = lnk[0];

140:      /* If free space is not available, make more free space */
141:      /* Double the amount of total space in the list */
142:      if (current_space->local_remaining<cnzi) {
143:        PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),&current_space);
144:        ndouble++;
145:      }

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

150:      current_space->array           += cnzi;
151:      current_space->local_used      += cnzi;
152:      current_space->local_remaining -= cnzi;

154:      ci[i+1] = ci[i] + cnzi;
155:    }

157:    /* Column indices are in the list of free space */
158:    /* Allocate space for cj, initialize cj, and */
159:    /* destroy list of free space and other temporary array(s) */
160:    PetscMalloc1(ci[am]+1,&cj);
161:    PetscFreeSpaceContiguous(&free_space,cj);
162:    PetscLLCondensedDestroy(lnk,lnkbt);

164:    /* put together the new symbolic matrix */
165:    MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
166:    MatSetBlockSizesFromMats(*C,A,B);
167:    MatSetType(*C,((PetscObject)A)->type_name);

169:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
170:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
171:   c                         = (Mat_SeqAIJ*)((*C)->data);
172:   c->free_a                 = PETSC_FALSE;
173:   c->free_ij                = PETSC_TRUE;
174:   c->nonew                  = 0;
175:   (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted; /* fast, needs non-scalable O(bn) array 'abdense' */

177:   /* set MatInfo */
178:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
179:   if (afill < 1.0) afill = 1.0;
180:   c->maxnz                     = ci[am];
181:   c->nz                        = ci[am];
182:   (*C)->info.mallocs           = ndouble;
183:   (*C)->info.fill_ratio_given  = fill;
184:   (*C)->info.fill_ratio_needed = afill;

186: #if defined(PETSC_USE_INFO)
187:   if (ci[am]) {
188:     PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
189:     PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
190:   } else {
191:     PetscInfo((*C),"Empty matrix product\n");
192:   }
193: #endif
194:   return(0);
195: }

197: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(Mat A,Mat B,Mat C)
198: {
200:   PetscLogDouble flops=0.0;
201:   Mat_SeqAIJ     *a   = (Mat_SeqAIJ*)A->data;
202:   Mat_SeqAIJ     *b   = (Mat_SeqAIJ*)B->data;
203:   Mat_SeqAIJ     *c   = (Mat_SeqAIJ*)C->data;
204:   PetscInt       *ai  =a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j;
205:   PetscInt       am   =A->rmap->n,cm=C->rmap->n;
206:   PetscInt       i,j,k,anzi,bnzi,cnzi,brow;
207:   PetscScalar    *aa=a->a,*ba=b->a,*baj,*ca,valtmp;
208:   PetscScalar    *ab_dense;

211:   if (!c->a) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
212:     PetscMalloc1(ci[cm]+1,&ca);
213:     c->a      = ca;
214:     c->free_a = PETSC_TRUE;
215:   } else {
216:     ca        = c->a;
217:   }
218:   if (!c->matmult_abdense) {
219:     PetscCalloc1(B->cmap->N,&ab_dense);
220:     c->matmult_abdense = ab_dense;
221:   } else {
222:     ab_dense = c->matmult_abdense;
223:   }

225:   /* clean old values in C */
226:   PetscMemzero(ca,ci[cm]*sizeof(MatScalar));
227:   /* Traverse A row-wise. */
228:   /* Build the ith row in C by summing over nonzero columns in A, */
229:   /* the rows of B corresponding to nonzeros of A. */
230:   for (i=0; i<am; i++) {
231:     anzi = ai[i+1] - ai[i];
232:     for (j=0; j<anzi; j++) {
233:       brow = aj[j];
234:       bnzi = bi[brow+1] - bi[brow];
235:       bjj  = bj + bi[brow];
236:       baj  = ba + bi[brow];
237:       /* perform dense axpy */
238:       valtmp = aa[j];
239:       for (k=0; k<bnzi; k++) {
240:         ab_dense[bjj[k]] += valtmp*baj[k];
241:       }
242:       flops += 2*bnzi;
243:     }
244:     aj += anzi; aa += anzi;

246:     cnzi = ci[i+1] - ci[i];
247:     for (k=0; k<cnzi; k++) {
248:       ca[k]          += ab_dense[cj[k]];
249:       ab_dense[cj[k]] = 0.0; /* zero ab_dense */
250:     }
251:     flops += cnzi;
252:     cj    += cnzi; ca += cnzi;
253:   }
254:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
255:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
256:   PetscLogFlops(flops);
257:   return(0);
258: }

260: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(Mat A,Mat B,Mat C)
261: {
263:   PetscLogDouble flops=0.0;
264:   Mat_SeqAIJ     *a   = (Mat_SeqAIJ*)A->data;
265:   Mat_SeqAIJ     *b   = (Mat_SeqAIJ*)B->data;
266:   Mat_SeqAIJ     *c   = (Mat_SeqAIJ*)C->data;
267:   PetscInt       *ai  = a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j;
268:   PetscInt       am   = A->rmap->N,cm=C->rmap->N;
269:   PetscInt       i,j,k,anzi,bnzi,cnzi,brow;
270:   PetscScalar    *aa=a->a,*ba=b->a,*baj,*ca=c->a,valtmp;
271:   PetscInt       nextb;

274:   if (!ca) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
275:     PetscMalloc1(ci[cm]+1,&ca);
276:     c->a      = ca;
277:     c->free_a = PETSC_TRUE;
278:   }

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

307:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
308:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
309:   PetscLogFlops(flops);
310:   return(0);
311: }

313: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(Mat A,Mat B,PetscReal fill,Mat *C)
314: {
315:   PetscErrorCode     ierr;
316:   Mat_SeqAIJ         *a  = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
317:   PetscInt           *ai = a->i,*bi=b->i,*ci,*cj;
318:   PetscInt           am  = A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
319:   MatScalar          *ca;
320:   PetscReal          afill;
321:   PetscInt           i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;
322:   PetscTable         ta;
323:   PetscFreeSpaceList free_space=NULL,current_space=NULL;

326:   /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_fast() */
327:   /*-----------------------------------------------------------------------------------------*/
328:   /* Allocate arrays for fill computation and free space for accumulating nonzero column */
329:   PetscMalloc1(am+2,&ci);
330:   ci[0] = 0;

332:   /* create and initialize a linked list */
333:   PetscTableCreate(bn,bn,&ta);
334:   MatRowMergeMax_SeqAIJ(b,bm,ta);
335:   PetscTableGetCount(ta,&Crmax);
336:   PetscTableDestroy(&ta);

338:   PetscLLCondensedCreate_fast(Crmax,&lnk);

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

344:   /* Determine ci and cj */
345:   for (i=0; i<am; i++) {
346:     anzi = ai[i+1] - ai[i];
347:     aj   = a->j + ai[i];
348:     for (j=0; j<anzi; j++) {
349:       brow = aj[j];
350:       bnzj = bi[brow+1] - bi[brow];
351:       bj   = b->j + bi[brow];
352:       /* add non-zero cols of B into the sorted linked list lnk */
353:       PetscLLCondensedAddSorted_fast(bnzj,bj,lnk);
354:     }
355:     cnzi = lnk[1];

357:     /* If free space is not available, make more free space */
358:     /* Double the amount of total space in the list */
359:     if (current_space->local_remaining<cnzi) {
360:       PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),&current_space);
361:       ndouble++;
362:     }

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

367:     current_space->array           += cnzi;
368:     current_space->local_used      += cnzi;
369:     current_space->local_remaining -= cnzi;

371:     ci[i+1] = ci[i] + cnzi;
372:   }

374:   /* Column indices are in the list of free space */
375:   /* Allocate space for cj, initialize cj, and */
376:   /* destroy list of free space and other temporary array(s) */
377:   PetscMalloc1(ci[am]+1,&cj);
378:   PetscFreeSpaceContiguous(&free_space,cj);
379:   PetscLLCondensedDestroy_fast(lnk);

381:   /* Allocate space for ca */
382:   PetscMalloc1(ci[am]+1,&ca);
383:   PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));

385:   /* put together the new symbolic matrix */
386:   MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,ca,C);
387:   MatSetBlockSizesFromMats(*C,A,B);
388:   MatSetType(*C,((PetscObject)A)->type_name);

390:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
391:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
392:   c          = (Mat_SeqAIJ*)((*C)->data);
393:   c->free_a  = PETSC_TRUE;
394:   c->free_ij = PETSC_TRUE;
395:   c->nonew   = 0;

397:   (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* slower, less memory */

399:   /* set MatInfo */
400:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
401:   if (afill < 1.0) afill = 1.0;
402:   c->maxnz                     = ci[am];
403:   c->nz                        = ci[am];
404:   (*C)->info.mallocs           = ndouble;
405:   (*C)->info.fill_ratio_given  = fill;
406:   (*C)->info.fill_ratio_needed = afill;

408: #if defined(PETSC_USE_INFO)
409:   if (ci[am]) {
410:     PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
411:     PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
412:   } else {
413:     PetscInfo((*C),"Empty matrix product\n");
414:   }
415: #endif
416:   return(0);
417: }


420: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(Mat A,Mat B,PetscReal fill,Mat *C)
421: {
422:   PetscErrorCode     ierr;
423:   Mat_SeqAIJ         *a  = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
424:   PetscInt           *ai = a->i,*bi=b->i,*ci,*cj;
425:   PetscInt           am  = A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
426:   MatScalar          *ca;
427:   PetscReal          afill;
428:   PetscInt           i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;
429:   PetscTable         ta;
430:   PetscFreeSpaceList free_space=NULL,current_space=NULL;

433:   /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_Scalalbe() */
434:   /*---------------------------------------------------------------------------------------------*/
435:   /* Allocate arrays for fill computation and free space for accumulating nonzero column */
436:   PetscMalloc1(am+2,&ci);
437:   ci[0] = 0;

439:   /* create and initialize a linked list */
440:   PetscTableCreate(bn,bn,&ta);
441:   MatRowMergeMax_SeqAIJ(b,bm,ta);
442:   PetscTableGetCount(ta,&Crmax);
443:   PetscTableDestroy(&ta);
444:   PetscLLCondensedCreate_Scalable(Crmax,&lnk);

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

450:   /* Determine ci and cj */
451:   for (i=0; i<am; i++) {
452:     anzi = ai[i+1] - ai[i];
453:     aj   = a->j + ai[i];
454:     for (j=0; j<anzi; j++) {
455:       brow = aj[j];
456:       bnzj = bi[brow+1] - bi[brow];
457:       bj   = b->j + bi[brow];
458:       /* add non-zero cols of B into the sorted linked list lnk */
459:       PetscLLCondensedAddSorted_Scalable(bnzj,bj,lnk);
460:     }
461:     cnzi = lnk[0];

463:     /* If free space is not available, make more free space */
464:     /* Double the amount of total space in the list */
465:     if (current_space->local_remaining<cnzi) {
466:       PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),&current_space);
467:       ndouble++;
468:     }

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

473:     current_space->array           += cnzi;
474:     current_space->local_used      += cnzi;
475:     current_space->local_remaining -= cnzi;

477:     ci[i+1] = ci[i] + cnzi;
478:   }

480:   /* Column indices are in the list of free space */
481:   /* Allocate space for cj, initialize cj, and */
482:   /* destroy list of free space and other temporary array(s) */
483:   PetscMalloc1(ci[am]+1,&cj);
484:   PetscFreeSpaceContiguous(&free_space,cj);
485:   PetscLLCondensedDestroy_Scalable(lnk);

487:   /* Allocate space for ca */
488:   /*-----------------------*/
489:   PetscMalloc1(ci[am]+1,&ca);
490:   PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));

492:   /* put together the new symbolic matrix */
493:   MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,ca,C);
494:   MatSetBlockSizesFromMats(*C,A,B);
495:   MatSetType(*C,((PetscObject)A)->type_name);

497:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
498:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
499:   c          = (Mat_SeqAIJ*)((*C)->data);
500:   c->free_a  = PETSC_TRUE;
501:   c->free_ij = PETSC_TRUE;
502:   c->nonew   = 0;

504:   (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* slower, less memory */

506:   /* set MatInfo */
507:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
508:   if (afill < 1.0) afill = 1.0;
509:   c->maxnz                     = ci[am];
510:   c->nz                        = ci[am];
511:   (*C)->info.mallocs           = ndouble;
512:   (*C)->info.fill_ratio_given  = fill;
513:   (*C)->info.fill_ratio_needed = afill;

515: #if defined(PETSC_USE_INFO)
516:   if (ci[am]) {
517:     PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
518:     PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
519:   } else {
520:     PetscInfo((*C),"Empty matrix product\n");
521:   }
522: #endif
523:   return(0);
524: }

526: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(Mat A,Mat B,PetscReal fill,Mat *C)
527: {
528:   PetscErrorCode     ierr;
529:   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
530:   const PetscInt     *ai=a->i,*bi=b->i,*aj=a->j,*bj=b->j;
531:   PetscInt           *ci,*cj,*bb;
532:   PetscInt           am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
533:   PetscReal          afill;
534:   PetscInt           i,j,col,ndouble = 0;
535:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
536:   PetscHeap          h;

539:   /* Get ci and cj - by merging sorted rows using a heap */
540:   /*---------------------------------------------------------------------------------------------*/
541:   /* Allocate arrays for fill computation and free space for accumulating nonzero column */
542:   PetscMalloc1(am+2,&ci);
543:   ci[0] = 0;

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

549:   PetscHeapCreate(a->rmax,&h);
550:   PetscMalloc1(a->rmax,&bb);

552:   /* Determine ci and cj */
553:   for (i=0; i<am; i++) {
554:     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 */
555:     const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
556:     ci[i+1] = ci[i];
557:     /* Populate the min heap */
558:     for (j=0; j<anzi; j++) {
559:       bb[j] = bi[acol[j]];         /* bb points at the start of the row */
560:       if (bb[j] < bi[acol[j]+1]) { /* Add if row is nonempty */
561:         PetscHeapAdd(h,j,bj[bb[j]++]);
562:       }
563:     }
564:     /* Pick off the min element, adding it to free space */
565:     PetscHeapPop(h,&j,&col);
566:     while (j >= 0) {
567:       if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
568:         PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2,current_space->total_array_size),16 << 20),&current_space);
569:         ndouble++;
570:       }
571:       *(current_space->array++) = col;
572:       current_space->local_used++;
573:       current_space->local_remaining--;
574:       ci[i+1]++;

576:       /* stash if anything else remains in this row of B */
577:       if (bb[j] < bi[acol[j]+1]) {PetscHeapStash(h,j,bj[bb[j]++]);}
578:       while (1) {               /* pop and stash any other rows of B that also had an entry in this column */
579:         PetscInt j2,col2;
580:         PetscHeapPeek(h,&j2,&col2);
581:         if (col2 != col) break;
582:         PetscHeapPop(h,&j2,&col2);
583:         if (bb[j2] < bi[acol[j2]+1]) {PetscHeapStash(h,j2,bj[bb[j2]++]);}
584:       }
585:       /* Put any stashed elements back into the min heap */
586:       PetscHeapUnstash(h);
587:       PetscHeapPop(h,&j,&col);
588:     }
589:   }
590:   PetscFree(bb);
591:   PetscHeapDestroy(&h);

593:   /* Column indices are in the list of free space */
594:   /* Allocate space for cj, initialize cj, and */
595:   /* destroy list of free space and other temporary array(s) */
596:   PetscMalloc1(ci[am],&cj);
597:   PetscFreeSpaceContiguous(&free_space,cj);

599:   /* put together the new symbolic matrix */
600:   MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
601:   MatSetBlockSizesFromMats(*C,A,B);
602:   MatSetType(*C,((PetscObject)A)->type_name);

604:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
605:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
606:   c          = (Mat_SeqAIJ*)((*C)->data);
607:   c->free_a  = PETSC_TRUE;
608:   c->free_ij = PETSC_TRUE;
609:   c->nonew   = 0;

611:   (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

613:   /* set MatInfo */
614:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
615:   if (afill < 1.0) afill = 1.0;
616:   c->maxnz                     = ci[am];
617:   c->nz                        = ci[am];
618:   (*C)->info.mallocs           = ndouble;
619:   (*C)->info.fill_ratio_given  = fill;
620:   (*C)->info.fill_ratio_needed = afill;

622: #if defined(PETSC_USE_INFO)
623:   if (ci[am]) {
624:     PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
625:     PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
626:   } else {
627:     PetscInfo((*C),"Empty matrix product\n");
628:   }
629: #endif
630:   return(0);
631: }

633: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(Mat A,Mat B,PetscReal fill,Mat *C)
634: {
635:   PetscErrorCode     ierr;
636:   Mat_SeqAIJ         *a  = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
637:   const PetscInt     *ai = a->i,*bi=b->i,*aj=a->j,*bj=b->j;
638:   PetscInt           *ci,*cj,*bb;
639:   PetscInt           am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
640:   PetscReal          afill;
641:   PetscInt           i,j,col,ndouble = 0;
642:   PetscFreeSpaceList free_space=NULL,current_space=NULL;
643:   PetscHeap          h;
644:   PetscBT            bt;

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

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

656:   current_space = free_space;

658:   PetscHeapCreate(a->rmax,&h);
659:   PetscMalloc1(a->rmax,&bb);
660:   PetscBTCreate(bn,&bt);

662:   /* Determine ci and cj */
663:   for (i=0; i<am; i++) {
664:     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 */
665:     const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
666:     const PetscInt *fptr = current_space->array; /* Save beginning of the row so we can clear the BT later */
667:     ci[i+1] = ci[i];
668:     /* Populate the min heap */
669:     for (j=0; j<anzi; j++) {
670:       PetscInt brow = acol[j];
671:       for (bb[j] = bi[brow]; bb[j] < bi[brow+1]; bb[j]++) {
672:         PetscInt bcol = bj[bb[j]];
673:         if (!PetscBTLookupSet(bt,bcol)) { /* new entry */
674:           PetscHeapAdd(h,j,bcol);
675:           bb[j]++;
676:           break;
677:         }
678:       }
679:     }
680:     /* Pick off the min element, adding it to free space */
681:     PetscHeapPop(h,&j,&col);
682:     while (j >= 0) {
683:       if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
684:         fptr = NULL;                      /* need PetscBTMemzero */
685:         PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2,current_space->total_array_size),16 << 20),&current_space);
686:         ndouble++;
687:       }
688:       *(current_space->array++) = col;
689:       current_space->local_used++;
690:       current_space->local_remaining--;
691:       ci[i+1]++;

693:       /* stash if anything else remains in this row of B */
694:       for (; bb[j] < bi[acol[j]+1]; bb[j]++) {
695:         PetscInt bcol = bj[bb[j]];
696:         if (!PetscBTLookupSet(bt,bcol)) { /* new entry */
697:           PetscHeapAdd(h,j,bcol);
698:           bb[j]++;
699:           break;
700:         }
701:       }
702:       PetscHeapPop(h,&j,&col);
703:     }
704:     if (fptr) {                 /* Clear the bits for this row */
705:       for (; fptr<current_space->array; fptr++) {PetscBTClear(bt,*fptr);}
706:     } else {                    /* We reallocated so we don't remember (easily) how to clear only the bits we changed */
707:       PetscBTMemzero(bn,bt);
708:     }
709:   }
710:   PetscFree(bb);
711:   PetscHeapDestroy(&h);
712:   PetscBTDestroy(&bt);

714:   /* Column indices are in the list of free space */
715:   /* Allocate space for cj, initialize cj, and */
716:   /* destroy list of free space and other temporary array(s) */
717:   PetscMalloc1(ci[am],&cj);
718:   PetscFreeSpaceContiguous(&free_space,cj);

720:   /* put together the new symbolic matrix */
721:   MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
722:   MatSetBlockSizesFromMats(*C,A,B);
723:   MatSetType(*C,((PetscObject)A)->type_name);

725:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
726:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
727:   c          = (Mat_SeqAIJ*)((*C)->data);
728:   c->free_a  = PETSC_TRUE;
729:   c->free_ij = PETSC_TRUE;
730:   c->nonew   = 0;

732:   (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

734:   /* set MatInfo */
735:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
736:   if (afill < 1.0) afill = 1.0;
737:   c->maxnz                     = ci[am];
738:   c->nz                        = ci[am];
739:   (*C)->info.mallocs           = ndouble;
740:   (*C)->info.fill_ratio_given  = fill;
741:   (*C)->info.fill_ratio_needed = afill;

743: #if defined(PETSC_USE_INFO)
744:   if (ci[am]) {
745:     PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
746:     PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
747:   } else {
748:     PetscInfo((*C),"Empty matrix product\n");
749:   }
750: #endif
751:   return(0);
752: }


755: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(Mat A,Mat B,PetscReal fill,Mat *C)
756: {
757:   PetscErrorCode     ierr;
758:   Mat_SeqAIJ         *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
759:   const PetscInt     *ai=a->i,*bi=b->i,*aj=a->j,*bj=b->j,*inputi,*inputj,*inputcol,*inputcol_L1;
760:   PetscInt           *ci,*cj,*outputj,worki_L1[9],worki_L2[9];
761:   PetscInt           c_maxmem,a_maxrownnz=0,a_rownnz;
762:   const PetscInt     workcol[8]={0,1,2,3,4,5,6,7};
763:   const PetscInt     am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
764:   const PetscInt     *brow_ptr[8],*brow_end[8];
765:   PetscInt           window[8];
766:   PetscInt           window_min,old_window_min,ci_nnz,outputi_nnz=0,L1_nrows,L2_nrows;
767:   PetscInt           i,k,ndouble=0,L1_rowsleft,rowsleft;
768:   PetscReal          afill;
769:   PetscInt           *workj_L1,*workj_L2,*workj_L3;
770:   PetscInt           L1_nnz,L2_nnz;

772:   /* Step 1: Get upper bound on memory required for allocation.
773:              Because of the way virtual memory works,
774:              only the memory pages that are actually needed will be physically allocated. */
776:   PetscMalloc1(am+1,&ci);
777:   for (i=0; i<am; i++) {
778:     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 */
779:     const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
780:     a_rownnz = 0;
781:     for (k=0; k<anzi; ++k) {
782:       a_rownnz += bi[acol[k]+1] - bi[acol[k]];
783:       if (a_rownnz > bn) {
784:         a_rownnz = bn;
785:         break;
786:       }
787:     }
788:     a_maxrownnz = PetscMax(a_maxrownnz, a_rownnz);
789:   }
790:   /* temporary work areas for merging rows */
791:   PetscMalloc1(a_maxrownnz*8,&workj_L1);
792:   PetscMalloc1(a_maxrownnz*8,&workj_L2);
793:   PetscMalloc1(a_maxrownnz,&workj_L3);

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

800:   ci_nnz       = 0;
801:   ci[0]        = 0;
802:   worki_L1[0]  = 0;
803:   worki_L2[0]  = 0;
804:   for (i=0; i<am; i++) {
805:     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 */
806:     const PetscInt *acol = aj + ai[i];      /* column indices of nonzero entries in this row */
807:     rowsleft             = anzi;
808:     inputcol_L1          = acol;
809:     L2_nnz               = 0;
810:     L2_nrows             = 1;  /* Number of rows to be merged on Level 3. output of L3 already exists -> initial value 1   */
811:     worki_L2[1]          = 0;
812:     outputi_nnz          = 0;

814:     /* If the number of indices in C so far + the max number of columns in the next row > c_maxmem  -> allocate more memory */
815:     while (ci_nnz+a_maxrownnz > c_maxmem) {
816:       c_maxmem *= 2;
817:       ndouble++;
818:       PetscRealloc(sizeof(PetscInt)*c_maxmem,&cj);
819:     }

821:     while (rowsleft) {
822:       L1_rowsleft = PetscMin(64, rowsleft); /* In the inner loop max 64 rows of B can be merged */
823:       L1_nrows    = 0;
824:       L1_nnz      = 0;
825:       inputcol    = inputcol_L1;
826:       inputi      = bi;
827:       inputj      = bj;

829:       /* The following macro is used to specialize for small rows in A.
830:          This helps with compiler unrolling, improving performance substantially.
831:           Input:  inputj   inputi  inputcol  bn
832:           Output: outputj  outputi_nnz                       */
833:        #define MatMatMultSymbolic_RowMergeMacro(ANNZ)                        \
834:          window_min  = bn;                                                   \
835:          outputi_nnz = 0;                                                    \
836:          for (k=0; k<ANNZ; ++k) {                                            \
837:            brow_ptr[k] = inputj + inputi[inputcol[k]];                       \
838:            brow_end[k] = inputj + inputi[inputcol[k]+1];                     \
839:            window[k]   = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn;   \
840:            window_min  = PetscMin(window[k], window_min);                    \
841:          }                                                                   \
842:          while (window_min < bn) {                                           \
843:            outputj[outputi_nnz++] = window_min;                              \
844:            /* advance front and compute new minimum */                       \
845:            old_window_min = window_min;                                      \
846:            window_min = bn;                                                  \
847:            for (k=0; k<ANNZ; ++k) {                                          \
848:              if (window[k] == old_window_min) {                              \
849:                brow_ptr[k]++;                                                \
850:                window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \
851:              }                                                               \
852:              window_min = PetscMin(window[k], window_min);                   \
853:            }                                                                 \
854:          }

856:       /************** L E V E L  1 ***************/
857:       /* Merge up to 8 rows of B to L1 work array*/
858:       while (L1_rowsleft) {
859:         outputi_nnz = 0;
860:         if (anzi > 8)  outputj = workj_L1 + L1_nnz;     /* Level 1 rowmerge*/
861:         else           outputj = cj + ci_nnz;           /* Merge directly to C */

863:         switch (L1_rowsleft) {
864:         case 1:  brow_ptr[0] = inputj + inputi[inputcol[0]];
865:                  brow_end[0] = inputj + inputi[inputcol[0]+1];
866:                  for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
867:                  inputcol    += L1_rowsleft;
868:                  rowsleft    -= L1_rowsleft;
869:                  L1_rowsleft  = 0;
870:                  break;
871:         case 2:  MatMatMultSymbolic_RowMergeMacro(2);
872:                  inputcol    += L1_rowsleft;
873:                  rowsleft    -= L1_rowsleft;
874:                  L1_rowsleft  = 0;
875:                  break;
876:         case 3: MatMatMultSymbolic_RowMergeMacro(3);
877:                  inputcol    += L1_rowsleft;
878:                  rowsleft    -= L1_rowsleft;
879:                  L1_rowsleft  = 0;
880:                  break;
881:         case 4:  MatMatMultSymbolic_RowMergeMacro(4);
882:                  inputcol    += L1_rowsleft;
883:                  rowsleft    -= L1_rowsleft;
884:                  L1_rowsleft  = 0;
885:                  break;
886:         case 5:  MatMatMultSymbolic_RowMergeMacro(5);
887:                  inputcol    += L1_rowsleft;
888:                  rowsleft    -= L1_rowsleft;
889:                  L1_rowsleft  = 0;
890:                  break;
891:         case 6:  MatMatMultSymbolic_RowMergeMacro(6);
892:                  inputcol    += L1_rowsleft;
893:                  rowsleft    -= L1_rowsleft;
894:                  L1_rowsleft  = 0;
895:                  break;
896:         case 7:  MatMatMultSymbolic_RowMergeMacro(7);
897:                  inputcol    += L1_rowsleft;
898:                  rowsleft    -= L1_rowsleft;
899:                  L1_rowsleft  = 0;
900:                  break;
901:         default: MatMatMultSymbolic_RowMergeMacro(8);
902:                  inputcol    += 8;
903:                  rowsleft    -= 8;
904:                  L1_rowsleft -= 8;
905:                  break;
906:         }
907:         inputcol_L1           = inputcol;
908:         L1_nnz               += outputi_nnz;
909:         worki_L1[++L1_nrows]  = L1_nnz;
910:       }

912:       /********************** L E V E L  2 ************************/
913:       /* Merge from L1 work array to either C or to L2 work array */
914:       if (anzi > 8) {
915:         inputi      = worki_L1;
916:         inputj      = workj_L1;
917:         inputcol    = workcol;
918:         outputi_nnz = 0;

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

923:         switch (L1_nrows) {
924:         case 1:  brow_ptr[0] = inputj + inputi[inputcol[0]];
925:                  brow_end[0] = inputj + inputi[inputcol[0]+1];
926:                  for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
927:                  break;
928:         case 2:  MatMatMultSymbolic_RowMergeMacro(2); break;
929:         case 3:  MatMatMultSymbolic_RowMergeMacro(3); break;
930:         case 4:  MatMatMultSymbolic_RowMergeMacro(4); break;
931:         case 5:  MatMatMultSymbolic_RowMergeMacro(5); break;
932:         case 6:  MatMatMultSymbolic_RowMergeMacro(6); break;
933:         case 7:  MatMatMultSymbolic_RowMergeMacro(7); break;
934:         case 8:  MatMatMultSymbolic_RowMergeMacro(8); break;
935:         default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatMult logic error: Not merging 1-8 rows from L1 work array!");
936:         }
937:         L2_nnz               += outputi_nnz;
938:         worki_L2[++L2_nrows]  = L2_nnz;

940:         /************************ L E V E L  3 **********************/
941:         /* Merge from L2 work array to either C or to L2 work array */
942:         if (anzi > 64 && (L2_nrows == 8 || rowsleft == 0)) {
943:           inputi      = worki_L2;
944:           inputj      = workj_L2;
945:           inputcol    = workcol;
946:           outputi_nnz = 0;
947:           if (rowsleft) outputj = workj_L3;
948:           else          outputj = cj + ci_nnz;
949:           switch (L2_nrows) {
950:           case 1:  brow_ptr[0] = inputj + inputi[inputcol[0]];
951:                    brow_end[0] = inputj + inputi[inputcol[0]+1];
952:                    for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
953:                    break;
954:           case 2:  MatMatMultSymbolic_RowMergeMacro(2); break;
955:           case 3:  MatMatMultSymbolic_RowMergeMacro(3); break;
956:           case 4:  MatMatMultSymbolic_RowMergeMacro(4); break;
957:           case 5:  MatMatMultSymbolic_RowMergeMacro(5); break;
958:           case 6:  MatMatMultSymbolic_RowMergeMacro(6); break;
959:           case 7:  MatMatMultSymbolic_RowMergeMacro(7); break;
960:           case 8:  MatMatMultSymbolic_RowMergeMacro(8); break;
961:           default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatMult logic error: Not merging 1-8 rows from L2 work array!");
962:           }
963:           L2_nrows    = 1;
964:           L2_nnz      = outputi_nnz;
965:           worki_L2[1] = outputi_nnz;
966:           /* Copy to workj_L2 */
967:           if (rowsleft) {
968:             for (k=0; k<outputi_nnz; ++k)  workj_L2[k] = outputj[k];
969:           }
970:         }
971:       }
972:     }  /* while (rowsleft) */
973: #undef MatMatMultSymbolic_RowMergeMacro

975:     /* terminate current row */
976:     ci_nnz += outputi_nnz;
977:     ci[i+1] = ci_nnz;
978:   }

980:   /* Step 3: Create the new symbolic matrix */
981:   MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
982:   MatSetBlockSizesFromMats(*C,A,B);
983:   MatSetType(*C,((PetscObject)A)->type_name);

985:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
986:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
987:   c          = (Mat_SeqAIJ*)((*C)->data);
988:   c->free_a  = PETSC_TRUE;
989:   c->free_ij = PETSC_TRUE;
990:   c->nonew   = 0;

992:   (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

994:   /* set MatInfo */
995:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
996:   if (afill < 1.0) afill = 1.0;
997:   c->maxnz                     = ci[am];
998:   c->nz                        = ci[am];
999:   (*C)->info.mallocs           = ndouble;
1000:   (*C)->info.fill_ratio_given  = fill;
1001:   (*C)->info.fill_ratio_needed = afill;

1003: #if defined(PETSC_USE_INFO)
1004:   if (ci[am]) {
1005:     PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
1006:     PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
1007:   } else {
1008:     PetscInfo((*C),"Empty matrix product\n");
1009:   }
1010: #endif

1012:   /* Step 4: Free temporary work areas */
1013:   PetscFree(workj_L1);
1014:   PetscFree(workj_L2);
1015:   PetscFree(workj_L3);
1016:   return(0);
1017: }

1019: /* concatenate unique entries and then sort */
1020: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(Mat A,Mat B,PetscReal fill,Mat *C)
1021: {
1022:   PetscErrorCode     ierr;
1023:   Mat_SeqAIJ         *a  = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
1024:   const PetscInt     *ai = a->i,*bi=b->i,*aj=a->j,*bj=b->j;
1025:   PetscInt           *ci,*cj;
1026:   PetscInt           am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
1027:   PetscReal          afill;
1028:   PetscInt           i,j,ndouble = 0;
1029:   PetscSegBuffer     seg,segrow;
1030:   char               *seen;

1033:   PetscMalloc1(am+1,&ci);
1034:   ci[0] = 0;

1036:   /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
1037:   PetscSegBufferCreate(sizeof(PetscInt),(PetscInt)(fill*(ai[am]+bi[bm])),&seg);
1038:   PetscSegBufferCreate(sizeof(PetscInt),100,&segrow);
1039:   PetscMalloc1(bn,&seen);
1040:   PetscMemzero(seen,bn*sizeof(char));

1042:   /* Determine ci and cj */
1043:   for (i=0; i<am; i++) {
1044:     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 */
1045:     const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
1046:     PetscInt packlen = 0,*PETSC_RESTRICT crow;
1047:     /* Pack segrow */
1048:     for (j=0; j<anzi; j++) {
1049:       PetscInt brow = acol[j],bjstart = bi[brow],bjend = bi[brow+1],k;
1050:       for (k=bjstart; k<bjend; k++) {
1051:         PetscInt bcol = bj[k];
1052:         if (!seen[bcol]) { /* new entry */
1053:           PetscInt *PETSC_RESTRICT slot;
1054:           PetscSegBufferGetInts(segrow,1,&slot);
1055:           *slot = bcol;
1056:           seen[bcol] = 1;
1057:           packlen++;
1058:         }
1059:       }
1060:     }
1061:     PetscSegBufferGetInts(seg,packlen,&crow);
1062:     PetscSegBufferExtractTo(segrow,crow);
1063:     PetscSortInt(packlen,crow);
1064:     ci[i+1] = ci[i] + packlen;
1065:     for (j=0; j<packlen; j++) seen[crow[j]] = 0;
1066:   }
1067:   PetscSegBufferDestroy(&segrow);
1068:   PetscFree(seen);

1070:   /* Column indices are in the segmented buffer */
1071:   PetscSegBufferExtractAlloc(seg,&cj);
1072:   PetscSegBufferDestroy(&seg);

1074:   /* put together the new symbolic matrix */
1075:   MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
1076:   MatSetBlockSizesFromMats(*C,A,B);
1077:   MatSetType(*C,((PetscObject)A)->type_name);

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

1086:   (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

1088:   /* set MatInfo */
1089:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
1090:   if (afill < 1.0) afill = 1.0;
1091:   c->maxnz                     = ci[am];
1092:   c->nz                        = ci[am];
1093:   (*C)->info.mallocs           = ndouble;
1094:   (*C)->info.fill_ratio_given  = fill;
1095:   (*C)->info.fill_ratio_needed = afill;

1097: #if defined(PETSC_USE_INFO)
1098:   if (ci[am]) {
1099:     PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
1100:     PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
1101:   } else {
1102:     PetscInfo((*C),"Empty matrix product\n");
1103:   }
1104: #endif
1105:   return(0);
1106: }

1108: /* This routine is not used. Should be removed! */
1109: PetscErrorCode MatMatTransposeMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1110: {

1114:   if (scall == MAT_INITIAL_MATRIX) {
1115:     PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
1116:     MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);
1117:     PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
1118:   }
1119:   PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
1120:   MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);
1121:   PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
1122:   return(0);
1123: }

1125: PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(Mat A)
1126: {
1127:   PetscErrorCode      ierr;
1128:   Mat_SeqAIJ          *a=(Mat_SeqAIJ*)A->data;
1129:   Mat_MatMatTransMult *abt=a->abt;

1132:   (abt->destroy)(A);
1133:   MatTransposeColoringDestroy(&abt->matcoloring);
1134:   MatDestroy(&abt->Bt_den);
1135:   MatDestroy(&abt->ABt_den);
1136:   PetscFree(abt);
1137:   return(0);
1138: }

1140: PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
1141: {
1142:   PetscErrorCode      ierr;
1143:   Mat                 Bt;
1144:   PetscInt            *bti,*btj;
1145:   Mat_MatMatTransMult *abt;
1146:   Mat_SeqAIJ          *c;

1149:   /* create symbolic Bt */
1150:   MatGetSymbolicTranspose_SeqAIJ(B,&bti,&btj);
1151:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,B->cmap->n,B->rmap->n,bti,btj,NULL,&Bt);
1152:   MatSetBlockSizes(Bt,PetscAbs(A->cmap->bs),PetscAbs(B->cmap->bs));
1153:   MatSetType(Bt,((PetscObject)A)->type_name);

1155:   /* get symbolic C=A*Bt */
1156:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,Bt,fill,C);

1158:   /* create a supporting struct for reuse intermidiate dense matrices with matcoloring */
1159:   PetscNew(&abt);
1160:   c      = (Mat_SeqAIJ*)(*C)->data;
1161:   c->abt = abt;

1163:   abt->usecoloring = PETSC_FALSE;
1164:   abt->destroy     = (*C)->ops->destroy;
1165:   (*C)->ops->destroy     = MatDestroy_SeqAIJ_MatMatMultTrans;

1167:   PetscOptionsGetBool(((PetscObject)A)->options,NULL,"-matmattransmult_color",&abt->usecoloring,NULL);
1168:   if (abt->usecoloring) {
1169:     /* Create MatTransposeColoring from symbolic C=A*B^T */
1170:     MatTransposeColoring matcoloring;
1171:     MatColoring          coloring;
1172:     ISColoring           iscoloring;
1173:     Mat                  Bt_dense,C_dense;
1174:     Mat_SeqAIJ           *c=(Mat_SeqAIJ*)(*C)->data;
1175:     /* inode causes memory problem, don't know why */
1176:     if (c->inode.use) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MAT_USE_INODES is not supported. Use '-mat_no_inode'");

1178:     MatColoringCreate(*C,&coloring);
1179:     MatColoringSetDistance(coloring,2);
1180:     MatColoringSetType(coloring,MATCOLORINGSL);
1181:     MatColoringSetFromOptions(coloring);
1182:     MatColoringApply(coloring,&iscoloring);
1183:     MatColoringDestroy(&coloring);
1184:     MatTransposeColoringCreate(*C,iscoloring,&matcoloring);

1186:     abt->matcoloring = matcoloring;

1188:     ISColoringDestroy(&iscoloring);

1190:     /* Create Bt_dense and C_dense = A*Bt_dense */
1191:     MatCreate(PETSC_COMM_SELF,&Bt_dense);
1192:     MatSetSizes(Bt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);
1193:     MatSetType(Bt_dense,MATSEQDENSE);
1194:     MatSeqDenseSetPreallocation(Bt_dense,NULL);

1196:     Bt_dense->assembled = PETSC_TRUE;
1197:     abt->Bt_den   = Bt_dense;

1199:     MatCreate(PETSC_COMM_SELF,&C_dense);
1200:     MatSetSizes(C_dense,A->rmap->n,matcoloring->ncolors,A->rmap->n,matcoloring->ncolors);
1201:     MatSetType(C_dense,MATSEQDENSE);
1202:     MatSeqDenseSetPreallocation(C_dense,NULL);

1204:     Bt_dense->assembled = PETSC_TRUE;
1205:     abt->ABt_den  = C_dense;

1207: #if defined(PETSC_USE_INFO)
1208:     {
1209:       Mat_SeqAIJ *c = (Mat_SeqAIJ*)(*C)->data;
1210:       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));
1211:     }
1212: #endif
1213:   }
1214:   /* clean up */
1215:   MatDestroy(&Bt);
1216:   MatRestoreSymbolicTranspose_SeqAIJ(B,&bti,&btj);
1217:   return(0);
1218: }

1220: PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
1221: {
1222:   PetscErrorCode      ierr;
1223:   Mat_SeqAIJ          *a   =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;
1224:   PetscInt            *ai  =a->i,*aj=a->j,*bi=b->i,*bj=b->j,anzi,bnzj,nexta,nextb,*acol,*bcol,brow;
1225:   PetscInt            cm   =C->rmap->n,*ci=c->i,*cj=c->j,i,j,cnzi,*ccol;
1226:   PetscLogDouble      flops=0.0;
1227:   MatScalar           *aa  =a->a,*aval,*ba=b->a,*bval,*ca,*cval;
1228:   Mat_MatMatTransMult *abt = c->abt;

1231:   /* clear old values in C */
1232:   if (!c->a) {
1233:     PetscMalloc1(ci[cm]+1,&ca);
1234:     c->a      = ca;
1235:     c->free_a = PETSC_TRUE;
1236:   } else {
1237:     ca =  c->a;
1238:   }
1239:   PetscMemzero(ca,ci[cm]*sizeof(MatScalar));

1241:   if (abt->usecoloring) {
1242:     MatTransposeColoring matcoloring = abt->matcoloring;
1243:     Mat                  Bt_dense,C_dense = abt->ABt_den;

1245:     /* Get Bt_dense by Apply MatTransposeColoring to B */
1246:     Bt_dense = abt->Bt_den;
1247:     MatTransColoringApplySpToDen(matcoloring,B,Bt_dense);

1249:     /* C_dense = A*Bt_dense */
1250:     MatMatMultNumeric_SeqAIJ_SeqDense(A,Bt_dense,C_dense);

1252:     /* Recover C from C_dense */
1253:     MatTransColoringApplyDenToSp(matcoloring,C_dense,C);
1254:     return(0);
1255:   }

1257:   for (i=0; i<cm; i++) {
1258:     anzi = ai[i+1] - ai[i];
1259:     acol = aj + ai[i];
1260:     aval = aa + ai[i];
1261:     cnzi = ci[i+1] - ci[i];
1262:     ccol = cj + ci[i];
1263:     cval = ca + ci[i];
1264:     for (j=0; j<cnzi; j++) {
1265:       brow = ccol[j];
1266:       bnzj = bi[brow+1] - bi[brow];
1267:       bcol = bj + bi[brow];
1268:       bval = ba + bi[brow];

1270:       /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */
1271:       nexta = 0; nextb = 0;
1272:       while (nexta<anzi && nextb<bnzj) {
1273:         while (nexta < anzi && acol[nexta] < bcol[nextb]) nexta++;
1274:         if (nexta == anzi) break;
1275:         while (nextb < bnzj && acol[nexta] > bcol[nextb]) nextb++;
1276:         if (nextb == bnzj) break;
1277:         if (acol[nexta] == bcol[nextb]) {
1278:           cval[j] += aval[nexta]*bval[nextb];
1279:           nexta++; nextb++;
1280:           flops += 2;
1281:         }
1282:       }
1283:     }
1284:   }
1285:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1286:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1287:   PetscLogFlops(flops);
1288:   return(0);
1289: }

1291: PetscErrorCode MatDestroy_SeqAIJ_MatTransMatMult(Mat A)
1292: {
1293:   PetscErrorCode      ierr;
1294:   Mat_SeqAIJ          *a = (Mat_SeqAIJ*)A->data;
1295:   Mat_MatTransMatMult *atb = a->atb;

1298:   if (atb) {
1299:     MatDestroy(&atb->At);
1300:     (*atb->destroy)(A);
1301:   }
1302:   PetscFree(atb);
1303:   return(0);
1304: }

1306: PetscErrorCode MatTransposeMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1307: {
1308:   PetscErrorCode      ierr;
1309:   const char          *algTypes[2] = {"matmatmult","outerproduct"};
1310:   PetscInt            alg=0; /* set default algorithm */
1311:   Mat                 At;
1312:   Mat_MatTransMatMult *atb;
1313:   Mat_SeqAIJ          *c;

1316:   if (scall == MAT_INITIAL_MATRIX) {
1317:     PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatTransposeMatMult","Mat");
1318:     PetscOptionsEList("-mattransposematmult_via","Algorithmic approach","MatTransposeMatMult",algTypes,2,algTypes[0],&alg,NULL);
1319:     PetscOptionsEnd();

1321:     switch (alg) {
1322:     case 1:
1323:       MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);
1324:       break;
1325:     default:
1326:       PetscNew(&atb);
1327:       MatTranspose_SeqAIJ(A,MAT_INITIAL_MATRIX,&At);
1328:       MatMatMult_SeqAIJ_SeqAIJ(At,B,MAT_INITIAL_MATRIX,fill,C);

1330:       c                  = (Mat_SeqAIJ*)(*C)->data;
1331:       c->atb             = atb;
1332:       atb->At            = At;
1333:       atb->destroy       = (*C)->ops->destroy;
1334:       (*C)->ops->destroy = MatDestroy_SeqAIJ_MatTransMatMult;

1336:       break;
1337:     }
1338:   }
1339:   if (alg) {
1340:     (*(*C)->ops->mattransposemultnumeric)(A,B,*C);
1341:   } else if (!alg && scall == MAT_REUSE_MATRIX) {
1342:     c   = (Mat_SeqAIJ*)(*C)->data;
1343:     atb = c->atb;
1344:     At  = atb->At;
1345:     MatTranspose_SeqAIJ(A,MAT_REUSE_MATRIX,&At);
1346:     MatMatMult_SeqAIJ_SeqAIJ(At,B,MAT_REUSE_MATRIX,fill,C);
1347:   }
1348:   return(0);
1349: }

1351: PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
1352: {
1354:   Mat            At;
1355:   PetscInt       *ati,*atj;

1358:   /* create symbolic At */
1359:   MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);
1360:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,A->cmap->n,A->rmap->n,ati,atj,NULL,&At);
1361:   MatSetBlockSizes(At,PetscAbs(A->cmap->bs),PetscAbs(B->cmap->bs));
1362:   MatSetType(At,((PetscObject)A)->type_name);

1364:   /* get symbolic C=At*B */
1365:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);

1367:   /* clean up */
1368:   MatDestroy(&At);
1369:   MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);

1371:   (*C)->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ;
1372:   return(0);
1373: }

1375: PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
1376: {
1378:   Mat_SeqAIJ     *a   =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;
1379:   PetscInt       am   =A->rmap->n,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb;
1380:   PetscInt       cm   =C->rmap->n,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k;
1381:   PetscLogDouble flops=0.0;
1382:   MatScalar      *aa  =a->a,*ba,*ca,*caj;

1385:   if (!c->a) {
1386:     PetscMalloc1(ci[cm]+1,&ca);

1388:     c->a      = ca;
1389:     c->free_a = PETSC_TRUE;
1390:   } else {
1391:     ca = c->a;
1392:   }
1393:   /* clear old values in C */
1394:   PetscMemzero(ca,ci[cm]*sizeof(MatScalar));

1396:   /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */
1397:   for (i=0; i<am; i++) {
1398:     bj   = b->j + bi[i];
1399:     ba   = b->a + bi[i];
1400:     bnzi = bi[i+1] - bi[i];
1401:     anzi = ai[i+1] - ai[i];
1402:     for (j=0; j<anzi; j++) {
1403:       nextb = 0;
1404:       crow  = *aj++;
1405:       cjj   = cj + ci[crow];
1406:       caj   = ca + ci[crow];
1407:       /* perform sparse axpy operation.  Note cjj includes bj. */
1408:       for (k=0; nextb<bnzi; k++) {
1409:         if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */
1410:           caj[k] += (*aa)*(*(ba+nextb));
1411:           nextb++;
1412:         }
1413:       }
1414:       flops += 2*bnzi;
1415:       aa++;
1416:     }
1417:   }

1419:   /* Assemble the final matrix and clean up */
1420:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1421:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1422:   PetscLogFlops(flops);
1423:   return(0);
1424: }

1426: PETSC_INTERN PetscErrorCode MatMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1427: {

1431:   if (scall == MAT_INITIAL_MATRIX) {
1432:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1433:     MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);
1434:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1435:   }
1436:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1437:   MatMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);
1438:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1439:   return(0);
1440: }

1442: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1443: {

1447:   MatMatMultSymbolic_SeqDense_SeqDense(A,B,0.0,C);

1449:   (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
1450:   return(0);
1451: }

1453: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
1454: {
1455:   Mat_SeqAIJ        *a=(Mat_SeqAIJ*)A->data;
1456:   Mat_SeqDense      *bd = (Mat_SeqDense*)B->data;
1457:   PetscErrorCode    ierr;
1458:   PetscScalar       *c,*b,r1,r2,r3,r4,*c1,*c2,*c3,*c4,aatmp;
1459:   const PetscScalar *aa,*b1,*b2,*b3,*b4;
1460:   const PetscInt    *aj;
1461:   PetscInt          cm=C->rmap->n,cn=B->cmap->n,bm=bd->lda,am=A->rmap->n;
1462:   PetscInt          am4=4*am,bm4=4*bm,col,i,j,n,ajtmp;

1465:   if (!cm || !cn) return(0);
1466:   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);
1467:   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);
1468:   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);
1469:   b = bd->v;
1470:   MatDenseGetArray(C,&c);
1471:   b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
1472:   c1 = c; c2 = c1 + am; c3 = c2 + am; c4 = c3 + am;
1473:   for (col=0; col<cn-4; col += 4) {  /* over columns of C */
1474:     for (i=0; i<am; i++) {        /* over rows of C in those columns */
1475:       r1 = r2 = r3 = r4 = 0.0;
1476:       n  = a->i[i+1] - a->i[i];
1477:       aj = a->j + a->i[i];
1478:       aa = a->a + a->i[i];
1479:       for (j=0; j<n; j++) {
1480:         aatmp = aa[j]; ajtmp = aj[j];
1481:         r1 += aatmp*b1[ajtmp];
1482:         r2 += aatmp*b2[ajtmp];
1483:         r3 += aatmp*b3[ajtmp];
1484:         r4 += aatmp*b4[ajtmp];
1485:       }
1486:       c1[i] = r1;
1487:       c2[i] = r2;
1488:       c3[i] = r3;
1489:       c4[i] = r4;
1490:     }
1491:     b1 += bm4; b2 += bm4; b3 += bm4; b4 += bm4;
1492:     c1 += am4; c2 += am4; c3 += am4; c4 += am4;
1493:   }
1494:   for (; col<cn; col++) {   /* over extra columns of C */
1495:     for (i=0; i<am; i++) {  /* over rows of C in those columns */
1496:       r1 = 0.0;
1497:       n  = a->i[i+1] - a->i[i];
1498:       aj = a->j + a->i[i];
1499:       aa = a->a + a->i[i];
1500:       for (j=0; j<n; j++) {
1501:         r1 += aa[j]*b1[aj[j]];
1502:       }
1503:       c1[i] = r1;
1504:     }
1505:     b1 += bm;
1506:     c1 += am;
1507:   }
1508:   PetscLogFlops(cn*(2.0*a->nz));
1509:   MatDenseRestoreArray(C,&c);
1510:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1511:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1512:   return(0);
1513: }

1515: /*
1516:    Note very similar to MatMult_SeqAIJ(), should generate both codes from same base
1517: */
1518: PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
1519: {
1520:   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data;
1521:   Mat_SeqDense   *bd = (Mat_SeqDense*)B->data;
1523:   PetscScalar    *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4;
1524:   MatScalar      *aa;
1525:   PetscInt       cm  = C->rmap->n, cn=B->cmap->n, bm=bd->lda, col, i,j,n,*aj, am = A->rmap->n,*ii,arm;
1526:   PetscInt       am2 = 2*am, am3 = 3*am,  bm4 = 4*bm,colam,*ridx;

1529:   if (!cm || !cn) return(0);
1530:   b = bd->v;
1531:   MatDenseGetArray(C,&c);
1532:   b1   = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;

1534:   if (a->compressedrow.use) { /* use compressed row format */
1535:     for (col=0; col<cn-4; col += 4) {  /* over columns of C */
1536:       colam = col*am;
1537:       arm   = a->compressedrow.nrows;
1538:       ii    = a->compressedrow.i;
1539:       ridx  = a->compressedrow.rindex;
1540:       for (i=0; i<arm; i++) {        /* over rows of C in those columns */
1541:         r1 = r2 = r3 = r4 = 0.0;
1542:         n  = ii[i+1] - ii[i];
1543:         aj = a->j + ii[i];
1544:         aa = a->a + ii[i];
1545:         for (j=0; j<n; j++) {
1546:           r1 += (*aa)*b1[*aj];
1547:           r2 += (*aa)*b2[*aj];
1548:           r3 += (*aa)*b3[*aj];
1549:           r4 += (*aa++)*b4[*aj++];
1550:         }
1551:         c[colam       + ridx[i]] += r1;
1552:         c[colam + am  + ridx[i]] += r2;
1553:         c[colam + am2 + ridx[i]] += r3;
1554:         c[colam + am3 + ridx[i]] += r4;
1555:       }
1556:       b1 += bm4;
1557:       b2 += bm4;
1558:       b3 += bm4;
1559:       b4 += bm4;
1560:     }
1561:     for (; col<cn; col++) {     /* over extra columns of C */
1562:       colam = col*am;
1563:       arm   = a->compressedrow.nrows;
1564:       ii    = a->compressedrow.i;
1565:       ridx  = a->compressedrow.rindex;
1566:       for (i=0; i<arm; i++) {  /* over rows of C in those columns */
1567:         r1 = 0.0;
1568:         n  = ii[i+1] - ii[i];
1569:         aj = a->j + ii[i];
1570:         aa = a->a + ii[i];

1572:         for (j=0; j<n; j++) {
1573:           r1 += (*aa++)*b1[*aj++];
1574:         }
1575:         c[colam + ridx[i]] += r1;
1576:       }
1577:       b1 += bm;
1578:     }
1579:   } else {
1580:     for (col=0; col<cn-4; col += 4) {  /* over columns of C */
1581:       colam = col*am;
1582:       for (i=0; i<am; i++) {        /* over rows of C in those columns */
1583:         r1 = r2 = r3 = r4 = 0.0;
1584:         n  = a->i[i+1] - a->i[i];
1585:         aj = a->j + a->i[i];
1586:         aa = a->a + a->i[i];
1587:         for (j=0; j<n; j++) {
1588:           r1 += (*aa)*b1[*aj];
1589:           r2 += (*aa)*b2[*aj];
1590:           r3 += (*aa)*b3[*aj];
1591:           r4 += (*aa++)*b4[*aj++];
1592:         }
1593:         c[colam + i]       += r1;
1594:         c[colam + am + i]  += r2;
1595:         c[colam + am2 + i] += r3;
1596:         c[colam + am3 + i] += r4;
1597:       }
1598:       b1 += bm4;
1599:       b2 += bm4;
1600:       b3 += bm4;
1601:       b4 += bm4;
1602:     }
1603:     for (; col<cn; col++) {     /* over extra columns of C */
1604:       colam = col*am;
1605:       for (i=0; i<am; i++) {  /* over rows of C in those columns */
1606:         r1 = 0.0;
1607:         n  = a->i[i+1] - a->i[i];
1608:         aj = a->j + a->i[i];
1609:         aa = a->a + a->i[i];

1611:         for (j=0; j<n; j++) {
1612:           r1 += (*aa++)*b1[*aj++];
1613:         }
1614:         c[colam + i] += r1;
1615:       }
1616:       b1 += bm;
1617:     }
1618:   }
1619:   PetscLogFlops(cn*2.0*a->nz);
1620:   MatDenseRestoreArray(C,&c);
1621:   return(0);
1622: }

1624: PetscErrorCode  MatTransColoringApplySpToDen_SeqAIJ(MatTransposeColoring coloring,Mat B,Mat Btdense)
1625: {
1627:   Mat_SeqAIJ     *b       = (Mat_SeqAIJ*)B->data;
1628:   Mat_SeqDense   *btdense = (Mat_SeqDense*)Btdense->data;
1629:   PetscInt       *bi      = b->i,*bj=b->j;
1630:   PetscInt       m        = Btdense->rmap->n,n=Btdense->cmap->n,j,k,l,col,anz,*btcol,brow,ncolumns;
1631:   MatScalar      *btval,*btval_den,*ba=b->a;
1632:   PetscInt       *columns=coloring->columns,*colorforcol=coloring->colorforcol,ncolors=coloring->ncolors;

1635:   btval_den=btdense->v;
1636:   PetscMemzero(btval_den,(m*n)*sizeof(MatScalar));
1637:   for (k=0; k<ncolors; k++) {
1638:     ncolumns = coloring->ncolumns[k];
1639:     for (l=0; l<ncolumns; l++) { /* insert a row of B to a column of Btdense */
1640:       col   = *(columns + colorforcol[k] + l);
1641:       btcol = bj + bi[col];
1642:       btval = ba + bi[col];
1643:       anz   = bi[col+1] - bi[col];
1644:       for (j=0; j<anz; j++) {
1645:         brow            = btcol[j];
1646:         btval_den[brow] = btval[j];
1647:       }
1648:     }
1649:     btval_den += m;
1650:   }
1651:   return(0);
1652: }

1654: PetscErrorCode MatTransColoringApplyDenToSp_SeqAIJ(MatTransposeColoring matcoloring,Mat Cden,Mat Csp)
1655: {
1657:   Mat_SeqAIJ     *csp = (Mat_SeqAIJ*)Csp->data;
1658:   PetscScalar    *ca_den,*ca_den_ptr,*ca=csp->a;
1659:   PetscInt       k,l,m=Cden->rmap->n,ncolors=matcoloring->ncolors;
1660:   PetscInt       brows=matcoloring->brows,*den2sp=matcoloring->den2sp;
1661:   PetscInt       nrows,*row,*idx;
1662:   PetscInt       *rows=matcoloring->rows,*colorforrow=matcoloring->colorforrow;

1665:   MatDenseGetArray(Cden,&ca_den);

1667:   if (brows > 0) {
1668:     PetscInt *lstart,row_end,row_start;
1669:     lstart = matcoloring->lstart;
1670:     PetscMemzero(lstart,ncolors*sizeof(PetscInt));

1672:     row_end = brows;
1673:     if (row_end > m) row_end = m;
1674:     for (row_start=0; row_start<m; row_start+=brows) { /* loop over row blocks of Csp */
1675:       ca_den_ptr = ca_den;
1676:       for (k=0; k<ncolors; k++) { /* loop over colors (columns of Cden) */
1677:         nrows = matcoloring->nrows[k];
1678:         row   = rows  + colorforrow[k];
1679:         idx   = den2sp + colorforrow[k];
1680:         for (l=lstart[k]; l<nrows; l++) {
1681:           if (row[l] >= row_end) {
1682:             lstart[k] = l;
1683:             break;
1684:           } else {
1685:             ca[idx[l]] = ca_den_ptr[row[l]];
1686:           }
1687:         }
1688:         ca_den_ptr += m;
1689:       }
1690:       row_end += brows;
1691:       if (row_end > m) row_end = m;
1692:     }
1693:   } else { /* non-blocked impl: loop over columns of Csp - slow if Csp is large */
1694:     ca_den_ptr = ca_den;
1695:     for (k=0; k<ncolors; k++) {
1696:       nrows = matcoloring->nrows[k];
1697:       row   = rows  + colorforrow[k];
1698:       idx   = den2sp + colorforrow[k];
1699:       for (l=0; l<nrows; l++) {
1700:         ca[idx[l]] = ca_den_ptr[row[l]];
1701:       }
1702:       ca_den_ptr += m;
1703:     }
1704:   }

1706:   MatDenseRestoreArray(Cden,&ca_den);
1707: #if defined(PETSC_USE_INFO)
1708:   if (matcoloring->brows > 0) {
1709:     PetscInfo1(Csp,"Loop over %D row blocks for den2sp\n",brows);
1710:   } else {
1711:     PetscInfo(Csp,"Loop over colors/columns of Cden, inefficient for large sparse matrix product \n");
1712:   }
1713: #endif
1714:   return(0);
1715: }

1717: PetscErrorCode MatTransposeColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatTransposeColoring c)
1718: {
1720:   PetscInt       i,n,nrows,Nbs,j,k,m,ncols,col,cm;
1721:   const PetscInt *is,*ci,*cj,*row_idx;
1722:   PetscInt       nis = iscoloring->n,*rowhit,bs = 1;
1723:   IS             *isa;
1724:   Mat_SeqAIJ     *csp = (Mat_SeqAIJ*)mat->data;
1725:   PetscInt       *colorforrow,*rows,*rows_i,*idxhit,*spidx,*den2sp,*den2sp_i;
1726:   PetscInt       *colorforcol,*columns,*columns_i,brows;
1727:   PetscBool      flg;

1730:   ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);

1732:   /* bs >1 is not being tested yet! */
1733:   Nbs       = mat->cmap->N/bs;
1734:   c->M      = mat->rmap->N/bs;  /* set total rows, columns and local rows */
1735:   c->N      = Nbs;
1736:   c->m      = c->M;
1737:   c->rstart = 0;
1738:   c->brows  = 100;

1740:   c->ncolors = nis;
1741:   PetscMalloc3(nis,&c->ncolumns,nis,&c->nrows,nis+1,&colorforrow);
1742:   PetscMalloc1(csp->nz+1,&rows);
1743:   PetscMalloc1(csp->nz+1,&den2sp);

1745:   brows = c->brows;
1746:   PetscOptionsGetInt(NULL,NULL,"-matden2sp_brows",&brows,&flg);
1747:   if (flg) c->brows = brows;
1748:   if (brows > 0) {
1749:     PetscMalloc1(nis+1,&c->lstart);
1750:   }

1752:   colorforrow[0] = 0;
1753:   rows_i         = rows;
1754:   den2sp_i       = den2sp;

1756:   PetscMalloc1(nis+1,&colorforcol);
1757:   PetscMalloc1(Nbs+1,&columns);

1759:   colorforcol[0] = 0;
1760:   columns_i      = columns;

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

1765:   cm   = c->m;
1766:   PetscMalloc1(cm+1,&rowhit);
1767:   PetscMalloc1(cm+1,&idxhit);
1768:   for (i=0; i<nis; i++) { /* loop over color */
1769:     ISGetLocalSize(isa[i],&n);
1770:     ISGetIndices(isa[i],&is);

1772:     c->ncolumns[i] = n;
1773:     if (n) {
1774:       PetscMemcpy(columns_i,is,n*sizeof(PetscInt));
1775:     }
1776:     colorforcol[i+1] = colorforcol[i] + n;
1777:     columns_i       += n;

1779:     /* fast, crude version requires O(N*N) work */
1780:     PetscMemzero(rowhit,cm*sizeof(PetscInt));

1782:     for (j=0; j<n; j++) { /* loop over columns*/
1783:       col     = is[j];
1784:       row_idx = cj + ci[col];
1785:       m       = ci[col+1] - ci[col];
1786:       for (k=0; k<m; k++) { /* loop over columns marking them in rowhit */
1787:         idxhit[*row_idx]   = spidx[ci[col] + k];
1788:         rowhit[*row_idx++] = col + 1;
1789:       }
1790:     }
1791:     /* count the number of hits */
1792:     nrows = 0;
1793:     for (j=0; j<cm; j++) {
1794:       if (rowhit[j]) nrows++;
1795:     }
1796:     c->nrows[i]      = nrows;
1797:     colorforrow[i+1] = colorforrow[i] + nrows;

1799:     nrows = 0;
1800:     for (j=0; j<cm; j++) { /* loop over rows */
1801:       if (rowhit[j]) {
1802:         rows_i[nrows]   = j;
1803:         den2sp_i[nrows] = idxhit[j];
1804:         nrows++;
1805:       }
1806:     }
1807:     den2sp_i += nrows;

1809:     ISRestoreIndices(isa[i],&is);
1810:     rows_i += nrows;
1811:   }
1812:   MatRestoreColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);
1813:   PetscFree(rowhit);
1814:   ISColoringRestoreIS(iscoloring,&isa);
1815:   if (csp->nz != colorforrow[nis]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"csp->nz %d != colorforrow[nis] %d",csp->nz,colorforrow[nis]);

1817:   c->colorforrow = colorforrow;
1818:   c->rows        = rows;
1819:   c->den2sp      = den2sp;
1820:   c->colorforcol = colorforcol;
1821:   c->columns     = columns;

1823:   PetscFree(idxhit);
1824:   return(0);
1825: }

1827: /* The combine method has done the symbolic and numeric in the first phase, and so we just return here */
1828: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Combined(Mat A,Mat B,Mat C)
1829: {
1831:   /* Empty function */
1832:   return(0);
1833: }

1835: /* This algorithm combines the symbolic and numeric phase of matrix-matrix multiplication. */
1836: PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ_Combined(Mat A,Mat B,PetscReal fill,Mat *C)
1837: {
1838:   PetscErrorCode     ierr;
1839:   PetscLogDouble     flops=0.0;
1840:   Mat_SeqAIJ         *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
1841:   const PetscInt     *ai=a->i,*bi=b->i;
1842:   PetscInt           *ci,*cj,*cj_i;
1843:   PetscScalar        *ca,*ca_i;
1844:   PetscInt           b_maxmemrow,c_maxmem,a_col;
1845:   PetscInt           am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
1846:   PetscInt           i,k,ndouble=0;
1847:   PetscReal          afill;
1848:   PetscScalar        *c_row_val_dense;
1849:   PetscBool          *c_row_idx_flags;
1850:   PetscInt           *aj_i=a->j;
1851:   PetscScalar        *aa_i=a->a;


1855:   /* Step 1: Determine upper bounds on memory for C and allocate memory */
1856:   /* This should be enough for almost all matrices. If still more memory is needed, it is reallocated later. */
1857:   c_maxmem    = 8*(ai[am]+bi[bm]);
1858:   b_maxmemrow = PetscMin(bi[bm],bn);
1859:   PetscMalloc1(am+1,&ci);
1860:   PetscMalloc1(bn,&c_row_val_dense);
1861:   PetscMalloc1(bn,&c_row_idx_flags);
1862:   PetscMalloc1(c_maxmem,&cj);
1863:   PetscMalloc1(c_maxmem,&ca);
1864:   ca_i  = ca;
1865:   cj_i  = cj;
1866:   ci[0] = 0;
1867:   PetscMemzero(c_row_val_dense,bn*sizeof(PetscScalar));
1868:   PetscMemzero(c_row_idx_flags,bn*sizeof(PetscBool));
1869:   for (i=0; i<am; i++) {
1870:     /* Step 2: Initialize the dense row vector for C  */
1871:     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 */
1872:     PetscInt       cnzi = 0;
1873:     PetscInt       *bj_i;
1874:     PetscScalar    *ba_i;
1875:     /* If the number of indices in C so far + the max number of columns in the next row > c_maxmem  -> allocate more memory
1876:        Usually, there is enough memory in the first place, so this is not executed. */
1877:     while (ci[i] + b_maxmemrow > c_maxmem) {
1878:       c_maxmem *= 2;
1879:       ndouble++;
1880:       PetscRealloc(sizeof(PetscInt)*c_maxmem,&cj);
1881:       PetscRealloc(sizeof(PetscScalar)*c_maxmem,&ca);
1882:     }

1884:     /* Step 3: Do the numerical calculations */
1885:     for (a_col=0; a_col<anzi; a_col++) {          /* iterate over all non zero values in a row of A */
1886:       PetscInt       a_col_index = aj_i[a_col];
1887:       const PetscInt bnzi        = bi[a_col_index+1] - bi[a_col_index];
1888:       flops += 2*bnzi;
1889:       bj_i   = b->j + bi[a_col_index];   /* points to the current row in bj */
1890:       ba_i   = b->a + bi[a_col_index];   /* points to the current row in ba */
1891:       for (k=0; k<bnzi; ++k) { /* iterate over all non zeros of this row in B */
1892:         if (c_row_idx_flags[bj_i[k]] == PETSC_FALSE) {
1893:           cj_i[cnzi++]             = bj_i[k];
1894:           c_row_idx_flags[bj_i[k]] = PETSC_TRUE;
1895:         }
1896:         c_row_val_dense[bj_i[k]] += aa_i[a_col] * ba_i[k];
1897:       }
1898:     }

1900:     /* Sort array */
1901:     PetscSortInt(cnzi,cj_i);
1902:     /* Step 4 */
1903:     for (k=0; k<cnzi; k++) {
1904:       ca_i[k]                  = c_row_val_dense[cj_i[k]];
1905:       c_row_val_dense[cj_i[k]] = 0.;
1906:       c_row_idx_flags[cj_i[k]] = PETSC_FALSE;
1907:     }
1908:     /* terminate current row */
1909:     aa_i   += anzi;
1910:     aj_i   += anzi;
1911:     ca_i   += cnzi;
1912:     cj_i   += cnzi;
1913:     ci[i+1] = ci[i] + cnzi;
1914:     flops  += cnzi;
1915:   }

1917:   /* Step 5 */
1918:   /* Create the new matrix */
1919:   MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
1920:   MatSetBlockSizesFromMats(*C,A,B);
1921:   MatSetType(*C,((PetscObject)A)->type_name);

1923:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
1924:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
1925:   c          = (Mat_SeqAIJ*)((*C)->data);
1926:   c->a       = ca;
1927:   c->free_a  = PETSC_TRUE;
1928:   c->free_ij = PETSC_TRUE;
1929:   c->nonew   = 0;

1931:   /* set MatInfo */
1932:   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
1933:   if (afill < 1.0) afill = 1.0;
1934:   c->maxnz                     = ci[am];
1935:   c->nz                        = ci[am];
1936:   (*C)->info.mallocs           = ndouble;
1937:   (*C)->info.fill_ratio_given  = fill;
1938:   (*C)->info.fill_ratio_needed = afill;
1939:   (*C)->ops->matmultnumeric    = MatMatMultNumeric_SeqAIJ_SeqAIJ_Combined;

1941:   PetscFree(c_row_val_dense);
1942:   PetscFree(c_row_idx_flags);

1944:   MatAssemblyBegin(*C,MAT_FINAL_ASSEMBLY);
1945:   MatAssemblyEnd(*C,MAT_FINAL_ASSEMBLY);
1946:   PetscLogFlops(flops);
1947:   return(0);
1948: }