2: /*
  3:   Defines matrix-matrix product routines 
  4:           C = A^T * B
  5: */

  7:  #include <../src/mat/impls/aij/seq/aij.h>
  8:  #include <../src/mat/impls/dense/seq/dense.h>

 10: PetscErrorCode MatDestroy_SeqDense_MatTransMatMult(Mat A)
 11: {
 12:   PetscErrorCode      ierr;
 13:   Mat_SeqDense        *a = (Mat_SeqDense*)A->data;
 14:   Mat_MatTransMatMult *atb = a->atb;

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

 25: PetscErrorCode MatTransposeMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
 26: {
 28: 
 30:   if (scall == MAT_INITIAL_MATRIX) {
 31:     PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
 32:     MatTransposeMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);
 33:     PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
 34:   }
 35:   PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
 36:   MatTransposeMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);
 37:   PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
 38:   return(0);
 39: }

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

 51:   PetscNew(&atb);

 53:   /* create output dense matrix C = A^T*B */
 54:   MatCreate(PETSC_COMM_SELF,&Cdense);
 55:   MatSetSizes(Cdense,n,BN,n,BN);
 56:   MatSetType(Cdense,MATSEQDENSE);
 57:   MatSeqDenseSetPreallocation(Cdense,NULL);

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

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

 77: PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
 78: {
 79:   PetscErrorCode      ierr;
 80:   PetscInt            i,j,k,m=A->rmap->n,n=A->cmap->n,BN=B->cmap->N;
 81:   const PetscScalar   *Barray,*ctarray;
 82:   PetscScalar         *Carray,*btarray;
 83:   Mat_SeqDense        *c=(Mat_SeqDense*)C->data;
 84:   Mat_MatTransMatMult *atb=c->atb;
 85:   Vec                 bt=atb->bt,ct=atb->ct;

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

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

 96:   k=0;
 97:   for (j=0; j<BN; j++) {
 98:     for (i=0; i<m; i++) btarray[i*BN + j] = Barray[k++];
 99:   }
100:   VecRestoreArray(bt,&btarray);
101:   MatDenseRestoreArrayRead(B,&Barray);
102: 
103:   /* compute ct = mA^T * cb */
104:   MatMultTranspose(atb->mA,bt,ct);

106:   /* transpose local arry of ct to matrix C */
107:   MatDenseGetArray(C,&Carray);
108:   VecGetArrayRead(ct,&ctarray);
109:   k = 0;
110:   for (j=0; j<BN; j++) {
111:     for (i=0; i<n; i++) Carray[k++] = ctarray[i*BN + j];
112:   }
113:   VecRestoreArrayRead(ct,&ctarray);
114:   MatDenseRestoreArray(C,&Carray);
115:   return(0);
116: }