Actual source code: matmatmult.c

      1: #define PETSCMAT_DLL

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

     8:  #include src/mat/impls/aij/seq/aij.h
     9:  #include src/mat/utils/freespace.h
     10:  #include petscbt.h

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

     20:   if (scall == MAT_INITIAL_MATRIX){
     21:     MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);
     22:   }
     23:   MatMatMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);
     24:   return(0);
     25: }

     30: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
     31: {
     33:   FreeSpaceList  free_space=PETSC_NULL,current_space=PETSC_NULL;
     34:   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
     35:   PetscInt       *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci,*cj;
     36:   PetscInt       am=A->M,bn=B->N,bm=B->M;
     37:   PetscInt       i,j,anzi,brow,bnzj,cnzi,nlnk,*lnk,nspacedouble=0;
     38:   MatScalar      *ca;
     39:   PetscBT        lnkbt;

     42:   /* Set up */
     43:   /* Allocate ci array, arrays for fill computation and */
     44:   /* free space for accumulating nonzero column info */
     45:   PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);
     46:   ci[0] = 0;
     47:
     48:   /* create and initialize a linked list */
     49:   nlnk = bn+1;
     50:   PetscLLCreate(bn,bn,nlnk,lnk,lnkbt);

     52:   /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
     53:   GetMoreSpace((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);
     54:   current_space = free_space;

     56:   /* Determine symbolic info for each row of the product: */
     57:   for (i=0;i<am;i++) {
     58:     anzi = ai[i+1] - ai[i];
     59:     cnzi = 0;
     60:     j    = anzi;
     61:     aj   = a->j + ai[i];
     62:     while (j){/* assume cols are almost in increasing order, starting from its end saves computation */
     63:       j--;
     64:       brow = *(aj + j);
     65:       bnzj = bi[brow+1] - bi[brow];
     66:       bjj  = bj + bi[brow];
     67:       /* add non-zero cols of B into the sorted linked list lnk */
     68:       PetscLLAdd(bnzj,bjj,bn,nlnk,lnk,lnkbt);
     69:       cnzi += nlnk;
     70:     }

     72:     /* If free space is not available, make more free space */
     73:     /* Double the amount of total space in the list */
     74:     if (current_space->local_remaining<cnzi) {
     75:       GetMoreSpace(current_space->total_array_size,&current_space);
     76:       nspacedouble++;
     77:     }

     79:     /* Copy data into free space, then initialize lnk */
     80:     PetscLLClean(bn,bn,cnzi,lnk,current_space->array,lnkbt);
     81:     current_space->array           += cnzi;
     82:     current_space->local_used      += cnzi;
     83:     current_space->local_remaining -= cnzi;

     85:     ci[i+1] = ci[i] + cnzi;
     86:   }

     88:   /* Column indices are in the list of free space */
     89:   /* Allocate space for cj, initialize cj, and */
     90:   /* destroy list of free space and other temporary array(s) */
     91:   PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);
     92:   MakeSpaceContiguous(&free_space,cj);
     93:   PetscLLDestroy(lnk,lnkbt);
     94:
     95:   /* Allocate space for ca */
     96:   PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);
     97:   PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));
     98:
     99:   /* put together the new symbolic matrix */
    100:   MatCreateSeqAIJWithArrays(A->comm,am,bn,ci,cj,ca,C);

    102:   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    103:   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
    104:   c = (Mat_SeqAIJ *)((*C)->data);
    105:   c->freedata = PETSC_TRUE;
    106:   c->nonew    = 0;

    108:   if (nspacedouble){
    109:     PetscLogInfo(((PetscObject)(*C),"MatMatMultSymbolic_SeqAIJ_SeqAIJ: nspacedouble:%D, nnz(A):%D, nnz(B):%D, fill:%g, nnz(C):%D\n",nspacedouble,ai[am],bi[bm],fill,ci[am]));
    110:   }
    111:   return(0);
    112: }

    117: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
    118: {
    120:   PetscInt       flops=0;
    121:   Mat_SeqAIJ     *a = (Mat_SeqAIJ *)A->data;
    122:   Mat_SeqAIJ     *b = (Mat_SeqAIJ *)B->data;
    123:   Mat_SeqAIJ     *c = (Mat_SeqAIJ *)C->data;
    124:   PetscInt       *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j;
    125:   PetscInt       am=A->M,cm=C->M;
    126:   PetscInt       i,j,k,anzi,bnzi,cnzi,brow,nextb;
    127:   MatScalar      *aa=a->a,*ba=b->a,*baj,*ca=c->a;

    130:   /* clean old values in C */
    131:   PetscMemzero(ca,ci[cm]*sizeof(MatScalar));
    132:   /* Traverse A row-wise. */
    133:   /* Build the ith row in C by summing over nonzero columns in A, */
    134:   /* the rows of B corresponding to nonzeros of A. */
    135:   for (i=0;i<am;i++) {
    136:     anzi = ai[i+1] - ai[i];
    137:     for (j=0;j<anzi;j++) {
    138:       brow = *aj++;
    139:       bnzi = bi[brow+1] - bi[brow];
    140:       bjj  = bj + bi[brow];
    141:       baj  = ba + bi[brow];
    142:       nextb = 0;
    143:       for (k=0; nextb<bnzi; k++) {
    144:         if (cj[k] == bjj[nextb]){ /* ccol == bcol */
    145:           ca[k] += (*aa)*baj[nextb++];
    146:         }
    147:       }
    148:       flops += 2*bnzi;
    149:       aa++;
    150:     }
    151:     cnzi = ci[i+1] - ci[i];
    152:     ca += cnzi;
    153:     cj += cnzi;
    154:   }
    155:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    156:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

    158:   PetscLogFlops(flops);
    159:   return(0);
    160: }

    165: PetscErrorCode MatMatMultTranspose_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) {

    169:   if (scall == MAT_INITIAL_MATRIX){
    170:     MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);
    171:   }
    172:   MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ(A,B,*C);
    173:   return(0);
    174: }

    178: PetscErrorCode MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
    179: {
    181:   Mat            At;
    182:   PetscInt       *ati,*atj;

    185:   /* create symbolic At */
    186:   MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);
    187:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,A->n,A->m,ati,atj,PETSC_NULL,&At);

    189:   /* get symbolic C=At*B */
    190:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);

    192:   /* clean up */
    193:   MatDestroy(At);
    194:   MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);
    195:
    196:   return(0);
    197: }

    201: PetscErrorCode MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
    202: {
    204:   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;
    205:   PetscInt       am=A->m,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb;
    206:   PetscInt       cm=C->m,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k,flops=0;
    207:   MatScalar      *aa=a->a,*ba,*ca=c->a,*caj;
    208:
    210:   /* clear old values in C */
    211:   PetscMemzero(ca,ci[cm]*sizeof(MatScalar));

    213:   /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */
    214:   for (i=0;i<am;i++) {
    215:     bj   = b->j + bi[i];
    216:     ba   = b->a + bi[i];
    217:     bnzi = bi[i+1] - bi[i];
    218:     anzi = ai[i+1] - ai[i];
    219:     for (j=0; j<anzi; j++) {
    220:       nextb = 0;
    221:       crow  = *aj++;
    222:       cjj   = cj + ci[crow];
    223:       caj   = ca + ci[crow];
    224:       /* perform sparse axpy operation.  Note cjj includes bj. */
    225:       for (k=0; nextb<bnzi; k++) {
    226:         if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */
    227:           caj[k] += (*aa)*(*(ba+nextb));
    228:           nextb++;
    229:         }
    230:       }
    231:       flops += 2*bnzi;
    232:       aa++;
    233:     }
    234:   }

    236:   /* Assemble the final matrix and clean up */
    237:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    238:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
    239:   PetscLogFlops(flops);
    240:   return(0);
    241: }