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(void *data) 12 { 13 PetscErrorCode ierr; 14 Mat_MatTransMatMult *atb = (Mat_MatTransMatMult*)data; 15 16 PetscFunctionBegin; 17 ierr = MatDestroy(&atb->mA);CHKERRQ(ierr); 18 ierr = VecDestroy(&atb->bt);CHKERRQ(ierr); 19 ierr = VecDestroy(&atb->ct);CHKERRQ(ierr); 20 ierr = PetscFree(atb);CHKERRQ(ierr); 21 PetscFunctionReturn(0); 22 } 23 24 static PetscErrorCode MatTransposeMatMultNumeric_MPIAIJ_MPIDense(Mat,Mat,Mat); 25 26 PETSC_INTERN PetscErrorCode MatTransposeMatMultSymbolic_MPIAIJ_MPIDense(Mat A,Mat B,PetscReal fill,Mat C) 27 { 28 PetscErrorCode ierr; 29 Mat_MatTransMatMult *atb; 30 PetscBool cisdense; 31 32 PetscFunctionBegin; 33 MatCheckProduct(C,4); 34 if (C->product->data) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Extra product struct not empty"); 35 36 /* create output dense matrix C = A^T*B */ 37 ierr = MatSetSizes(C,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N);CHKERRQ(ierr); 38 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATMPIDENSE,MATMPIDENSECUDA,"");CHKERRQ(ierr); 39 if (!cisdense) { 40 ierr = MatSetType(C,((PetscObject)B)->type_name);CHKERRQ(ierr); 41 } 42 ierr = MatSetUp(C);CHKERRQ(ierr); 43 44 /* create additional data structure for the product */ 45 ierr = PetscNew(&atb);CHKERRQ(ierr); 46 if (B->cmap->N) { 47 ierr = MatCreateMAIJ(A,B->cmap->N,&atb->mA);CHKERRQ(ierr); 48 if (!atb->mA->assembled) { 49 ierr = MatAssemblyBegin(atb->mA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 50 ierr = MatAssemblyEnd(atb->mA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 51 } 52 ierr = MatCreateVecs(atb->mA,&atb->ct,&atb->bt);CHKERRQ(ierr); 53 } 54 C->product->data = atb; 55 C->product->destroy = MatDestroy_MPIDense_MatTransMatMult; 56 57 C->ops->transposematmultnumeric = MatTransposeMatMultNumeric_MPIAIJ_MPIDense; 58 PetscFunctionReturn(0); 59 } 60 61 static PetscErrorCode MatTransposeMatMultNumeric_MPIAIJ_MPIDense(Mat A,Mat B,Mat C) 62 { 63 PetscErrorCode ierr; 64 const PetscScalar *Barray,*ctarray; 65 PetscScalar *Carray,*btarray; 66 Mat_MPIDense *b=(Mat_MPIDense*)B->data,*c=(Mat_MPIDense*)C->data; 67 Mat_SeqDense *bseq=(Mat_SeqDense*)(b->A)->data,*cseq=(Mat_SeqDense*)(c->A)->data; 68 PetscInt i,j,m=A->rmap->n,n=A->cmap->n,ldb=bseq->lda,BN=B->cmap->N,ldc=cseq->lda; 69 Mat_MatTransMatMult *atb; 70 Vec bt,ct; 71 72 PetscFunctionBegin; 73 MatCheckProduct(C,3); 74 atb=(Mat_MatTransMatMult *)C->product->data; 75 if (!atb) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing product struct"); 76 if (!BN) { 77 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 78 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 79 PetscFunctionReturn(0); 80 } 81 bt = atb->bt; 82 ct = atb->ct; 83 /* transpose local arry of B, then copy it to vector bt */ 84 ierr = MatDenseGetArrayRead(B,&Barray);CHKERRQ(ierr); 85 ierr = VecGetArray(bt,&btarray);CHKERRQ(ierr); 86 87 for (j=0; j<BN; j++) { 88 for (i=0; i<m; i++) btarray[i*BN + j] = Barray[j*ldb + i]; 89 } 90 ierr = VecRestoreArray(bt,&btarray);CHKERRQ(ierr); 91 ierr = MatDenseRestoreArrayRead(B,&Barray);CHKERRQ(ierr); 92 93 /* compute ct = mA^T * cb */ 94 ierr = MatMultTranspose(atb->mA,bt,ct);CHKERRQ(ierr); 95 96 /* transpose local array of ct to matrix C */ 97 ierr = MatDenseGetArray(C,&Carray);CHKERRQ(ierr); 98 ierr = VecGetArrayRead(ct,&ctarray);CHKERRQ(ierr); 99 100 for (j=0; j<BN; j++) { 101 for (i=0; i<n; i++) Carray[j*ldc + i] = ctarray[i*BN + j]; 102 } 103 ierr = VecRestoreArrayRead(ct,&ctarray);CHKERRQ(ierr); 104 ierr = MatDenseRestoreArray(C,&Carray);CHKERRQ(ierr); 105 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 106 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 107 PetscFunctionReturn(0); 108 } 109