2: /*
  3:   Defines matrix-matrix product routines for pairs of MPIAIJ matrices
  4:           C = A^T * B
  5:   The routines are slightly modified from MatTransposeMatMultxxx_SeqAIJ_SeqDense().
  6: */
  7:  #include <../src/mat/impls/aij/seq/aij.h>
  8:  #include <../src/mat/impls/aij/mpi/mpiaij.h>
  9:  #include <../src/mat/impls/dense/mpi/mpidense.h>

 11: PetscErrorCode MatDestroy_MPIDense_MatTransMatMult(Mat A)
 12: {
 13:   PetscErrorCode      ierr;
 14:   Mat_MPIDense        *a = (Mat_MPIDense*)A->data;
 15:   Mat_MatTransMatMult *atb = a->atb;

 18:   MatDestroy(&atb->mA);
 19:   VecDestroy(&atb->bt);
 20:   VecDestroy(&atb->ct);
 21:   (atb->destroy)(A);
 22:   PetscFree(atb);
 23:   return(0);
 24: }

 26: PetscErrorCode MatTransposeMatMult_MPIAIJ_MPIDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
 27: {

 31:   if (scall == MAT_INITIAL_MATRIX) {
 32:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
 33:     MatTransposeMatMultSymbolic_MPIAIJ_MPIDense(A,B,fill,C);
 34:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
 35:   }
 36:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
 37:   MatTransposeMatMultNumeric_MPIAIJ_MPIDense(A,B,*C);
 38:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
 39:   return(0);
 40: }

 42: PetscErrorCode MatTransposeMatMultSymbolic_MPIAIJ_MPIDense(Mat A,Mat B,PetscReal fill,Mat *C)
 43: {
 44:   PetscErrorCode      ierr;
 45:   PetscInt            m=A->rmap->n,n=A->cmap->n,BN=B->cmap->N;
 46:   Mat_MatTransMatMult *atb;
 47:   Mat                 Cdense;
 48:   Vec                 bt,ct;
 49:   Mat_MPIDense        *c;

 52:   PetscNew(&atb);

 54:   /* create output dense matrix C = A^T*B */
 55:   MatCreate(PetscObjectComm((PetscObject)A),&Cdense);
 56:   MatSetSizes(Cdense,n,PETSC_DECIDE,PETSC_DECIDE,BN);
 57:   MatSetType(Cdense,MATMPIDENSE);
 58:   MatMPIDenseSetPreallocation(Cdense,NULL);

 60:   /* create vectors bt and ct to hold locally transposed arrays of B and C */
 61:   VecCreate(PetscObjectComm((PetscObject)A),&bt);
 62:   VecSetSizes(bt,m*BN,PETSC_DECIDE);
 63:   VecSetType(bt,VECSTANDARD);
 64:   VecCreate(PetscObjectComm((PetscObject)A),&ct);
 65:   VecSetSizes(ct,n*BN,PETSC_DECIDE);
 66:   VecSetType(ct,VECSTANDARD);
 67:   atb->bt = bt;
 68:   atb->ct = ct;

 70:   *C = Cdense;
 71:   c                    = (Mat_MPIDense*)Cdense->data;
 72:   c->atb               = atb;
 73:   atb->destroy         = Cdense->ops->destroy;
 74:   Cdense->ops->destroy = MatDestroy_MPIDense_MatTransMatMult;
 75:   return(0);
 76: }

 78: PetscErrorCode MatTransposeMatMultNumeric_MPIAIJ_MPIDense(Mat A,Mat B,Mat C)
 79: {
 80:   PetscErrorCode      ierr;
 81:   const PetscScalar   *Barray,*ctarray;
 82:   PetscScalar         *Carray,*btarray;
 83:   Mat_MPIDense        *b=(Mat_MPIDense*)B->data,*c=(Mat_MPIDense*)C->data;
 84:   Mat_SeqDense        *bseq=(Mat_SeqDense*)(b->A)->data,*cseq=(Mat_SeqDense*)(c->A)->data;
 85:   PetscInt            i,j,m=A->rmap->n,n=A->cmap->n,ldb=bseq->lda,BN=B->cmap->N,ldc=cseq->lda;
 86:   Mat_MatTransMatMult *atb=c->atb;
 87:   Vec                 bt=atb->bt,ct=atb->ct;

 90:   /* create MAIJ matrix mA from A -- should be done in symbolic phase */
 91:   MatDestroy(&atb->mA);
 92:   MatCreateMAIJ(A,BN,&atb->mA);

 94:   /* transpose local arry of B, then copy it to vector bt */
 95:   MatDenseGetArrayRead(B,&Barray);
 96:   VecGetArray(bt,&btarray);

 98:   for (j=0; j<BN; j++) {
 99:     for (i=0; i<m; i++) btarray[i*BN + j] = Barray[j*ldb + i];
100:   }
101:   VecRestoreArray(bt,&btarray);
102:   MatDenseRestoreArrayRead(B,&Barray);

104:   /* compute ct = mA^T * cb */
105:   MatMultTranspose(atb->mA,bt,ct);

107:   /* transpose local array of ct to matrix C */
108:   MatDenseGetArray(C,&Carray);
109:   VecGetArrayRead(ct,&ctarray);

111:   for (j=0; j<BN; j++) {
112:     for (i=0; i<n; i++) Carray[j*ldc + i] = ctarray[i*BN + j];
113:   }
114:   VecRestoreArrayRead(ct,&ctarray);
115:   MatDenseRestoreArray(C,&Carray);
116:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
117:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
118:   return(0);
119: }