1 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> /*I "petscmat.h" I*/ 8 #include <../src/mat/impls/aij/mpi/mpiaij.h> 9 #include <../src/mat/impls/dense/mpi/mpidense.h> 10 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; 16 17 PetscFunctionBegin; 18 ierr = MatDestroy(&atb->mA);CHKERRQ(ierr); 19 ierr = VecDestroy(&atb->bt);CHKERRQ(ierr); 20 ierr = VecDestroy(&atb->ct);CHKERRQ(ierr); 21 ierr = (atb->destroy)(A);CHKERRQ(ierr); 22 ierr = PetscFree(atb);CHKERRQ(ierr); 23 PetscFunctionReturn(0); 24 } 25 26 PetscErrorCode MatTransposeMatMultSymbolic_MPIAIJ_MPIDense(Mat A,Mat B,PetscReal fill,Mat C) 27 { 28 PetscErrorCode ierr; 29 PetscInt m=A->rmap->n,n=A->cmap->n,BN=B->cmap->N; 30 Mat_MatTransMatMult *atb; 31 Vec bt,ct; 32 Mat_MPIDense *c; 33 34 PetscFunctionBegin; 35 ierr = PetscNew(&atb);CHKERRQ(ierr); 36 37 /* create output dense matrix C = A^T*B */ 38 ierr = MatSetSizes(C,n,B->cmap->n,A->cmap->N,BN);CHKERRQ(ierr); 39 ierr = MatSetType(C,MATMPIDENSE);CHKERRQ(ierr); 40 ierr = MatSetUp(C);CHKERRQ(ierr); 41 42 /* create vectors bt and ct to hold locally transposed arrays of B and C */ 43 ierr = VecCreate(PetscObjectComm((PetscObject)A),&bt);CHKERRQ(ierr); 44 ierr = VecSetSizes(bt,m*BN,PETSC_DECIDE);CHKERRQ(ierr); 45 ierr = VecSetType(bt,VECSTANDARD);CHKERRQ(ierr); 46 ierr = VecCreate(PetscObjectComm((PetscObject)A),&ct);CHKERRQ(ierr); 47 ierr = VecSetSizes(ct,n*BN,PETSC_DECIDE);CHKERRQ(ierr); 48 ierr = VecSetType(ct,VECSTANDARD);CHKERRQ(ierr); 49 atb->bt = bt; 50 atb->ct = ct; 51 52 c = (Mat_MPIDense*)C->data; 53 c->atb = atb; 54 atb->destroy = C->ops->destroy; 55 C->ops->destroy = MatDestroy_MPIDense_MatTransMatMult; 56 C->ops->transposematmultnumeric = MatTransposeMatMultNumeric_MPIAIJ_MPIDense; 57 PetscFunctionReturn(0); 58 } 59 60 PetscErrorCode MatTransposeMatMultNumeric_MPIAIJ_MPIDense(Mat A,Mat B,Mat C) 61 { 62 PetscErrorCode ierr; 63 const PetscScalar *Barray,*ctarray; 64 PetscScalar *Carray,*btarray; 65 Mat_MPIDense *b=(Mat_MPIDense*)B->data,*c=(Mat_MPIDense*)C->data; 66 Mat_SeqDense *bseq=(Mat_SeqDense*)(b->A)->data,*cseq=(Mat_SeqDense*)(c->A)->data; 67 PetscInt i,j,m=A->rmap->n,n=A->cmap->n,ldb=bseq->lda,BN=B->cmap->N,ldc=cseq->lda; 68 Mat_MatTransMatMult *atb=c->atb; 69 Vec bt=atb->bt,ct=atb->ct; 70 71 PetscFunctionBegin; 72 /* create MAIJ matrix mA from A -- should be done in symbolic phase */ 73 ierr = MatDestroy(&atb->mA);CHKERRQ(ierr); 74 ierr = MatCreateMAIJ(A,BN,&atb->mA);CHKERRQ(ierr); 75 76 /* transpose local arry of B, then copy it to vector bt */ 77 ierr = MatDenseGetArrayRead(B,&Barray);CHKERRQ(ierr); 78 ierr = VecGetArray(bt,&btarray);CHKERRQ(ierr); 79 80 for (j=0; j<BN; j++) { 81 for (i=0; i<m; i++) btarray[i*BN + j] = Barray[j*ldb + i]; 82 } 83 ierr = VecRestoreArray(bt,&btarray);CHKERRQ(ierr); 84 ierr = MatDenseRestoreArrayRead(B,&Barray);CHKERRQ(ierr); 85 86 /* compute ct = mA^T * cb */ 87 ierr = MatMultTranspose(atb->mA,bt,ct);CHKERRQ(ierr); 88 89 /* transpose local array of ct to matrix C */ 90 ierr = MatDenseGetArray(C,&Carray);CHKERRQ(ierr); 91 ierr = VecGetArrayRead(ct,&ctarray);CHKERRQ(ierr); 92 93 for (j=0; j<BN; j++) { 94 for (i=0; i<n; i++) Carray[j*ldc + i] = ctarray[i*BN + j]; 95 } 96 ierr = VecRestoreArrayRead(ct,&ctarray);CHKERRQ(ierr); 97 ierr = MatDenseRestoreArray(C,&Carray);CHKERRQ(ierr); 98 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 99 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 100 PetscFunctionReturn(0); 101 } 102