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 #undef __FUNCT__ 12 #define __FUNCT__ "MatDestroy_MPIDense_MatTransMatMult" 13 PetscErrorCode MatDestroy_MPIDense_MatTransMatMult(Mat A) 14 { 15 PetscErrorCode ierr; 16 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 17 Mat_MatTransMatMult *atb = a->atb; 18 19 PetscFunctionBegin; 20 ierr = MatDestroy(&atb->mA);CHKERRQ(ierr); 21 ierr = VecDestroy(&atb->bt);CHKERRQ(ierr); 22 ierr = VecDestroy(&atb->ct);CHKERRQ(ierr); 23 ierr = (atb->destroy)(A);CHKERRQ(ierr); 24 ierr = PetscFree(atb);CHKERRQ(ierr); 25 PetscFunctionReturn(0); 26 } 27 28 #undef __FUNCT__ 29 #define __FUNCT__ "MatTransposeMatMult_MPIAIJ_MPIDense" 30 PetscErrorCode MatTransposeMatMult_MPIAIJ_MPIDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 31 { 32 PetscErrorCode ierr; 33 34 PetscFunctionBegin; 35 if (scall == MAT_INITIAL_MATRIX) { 36 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 37 ierr = MatTransposeMatMultSymbolic_MPIAIJ_MPIDense(A,B,fill,C);CHKERRQ(ierr); 38 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 39 } 40 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 41 ierr = MatTransposeMatMultNumeric_MPIAIJ_MPIDense(A,B,*C);CHKERRQ(ierr); 42 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 46 #undef __FUNCT__ 47 #define __FUNCT__ "MatTransposeMatMultSymbolic_MPIAIJ_MPIDense" 48 PetscErrorCode MatTransposeMatMultSymbolic_MPIAIJ_MPIDense(Mat A,Mat B,PetscReal fill,Mat *C) 49 { 50 PetscErrorCode ierr; 51 PetscInt m=A->rmap->n,n=A->cmap->n,BN=B->cmap->N; 52 Mat_MatTransMatMult *atb; 53 Mat Cdense; 54 Vec bt,ct; 55 Mat_MPIDense *c; 56 57 PetscFunctionBegin; 58 ierr = PetscNew(&atb);CHKERRQ(ierr); 59 60 /* create output dense matrix C = A^T*B */ 61 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cdense);CHKERRQ(ierr); 62 ierr = MatSetSizes(Cdense,n,PETSC_DECIDE,PETSC_DECIDE,BN);CHKERRQ(ierr); 63 ierr = MatSetType(Cdense,MATMPIDENSE);CHKERRQ(ierr); 64 ierr = MatMPIDenseSetPreallocation(Cdense,NULL);CHKERRQ(ierr); 65 66 /* create vectors bt and ct to hold locally transposed arrays of B and C */ 67 ierr = VecCreate(PetscObjectComm((PetscObject)A),&bt);CHKERRQ(ierr); 68 ierr = VecSetSizes(bt,m*BN,PETSC_DECIDE);CHKERRQ(ierr); 69 ierr = VecSetType(bt,VECSTANDARD);CHKERRQ(ierr); 70 ierr = VecCreate(PetscObjectComm((PetscObject)A),&ct);CHKERRQ(ierr); 71 ierr = VecSetSizes(ct,n*BN,PETSC_DECIDE);CHKERRQ(ierr); 72 ierr = VecSetType(ct,VECSTANDARD);CHKERRQ(ierr); 73 atb->bt = bt; 74 atb->ct = ct; 75 76 *C = Cdense; 77 c = (Mat_MPIDense*)Cdense->data; 78 c->atb = atb; 79 atb->destroy = Cdense->ops->destroy; 80 Cdense->ops->destroy = MatDestroy_MPIDense_MatTransMatMult; 81 PetscFunctionReturn(0); 82 } 83 84 #undef __FUNCT__ 85 #define __FUNCT__ "MatTransposeMatMultNumeric_MPIAIJ_MPIDense" 86 PetscErrorCode MatTransposeMatMultNumeric_MPIAIJ_MPIDense(Mat A,Mat B,Mat C) 87 { 88 PetscErrorCode ierr; 89 PetscInt i,j,k,m=A->rmap->n,n=A->cmap->n,BN=B->cmap->N; 90 PetscScalar *Barray,*Carray,*btarray,*ctarray; 91 Mat_MPIDense *c=(Mat_MPIDense*)C->data; 92 Mat_MatTransMatMult *atb=c->atb; 93 Vec bt=atb->bt,ct=atb->ct; 94 95 PetscFunctionBegin; 96 /* create MAIJ matrix mA from A -- should be done in symbolic phase */ 97 ierr = MatDestroy(&atb->mA);CHKERRQ(ierr); 98 ierr = MatCreateMAIJ(A,BN,&atb->mA);CHKERRQ(ierr); 99 100 /* transpose local arry of B, then copy it to vector bt */ 101 ierr = MatDenseGetArray(B,&Barray);CHKERRQ(ierr); 102 ierr = VecGetArray(bt,&btarray);CHKERRQ(ierr); 103 104 k=0; 105 for (j=0; j<BN; j++) { 106 for (i=0; i<m; i++) btarray[i*BN + j] = Barray[k++]; 107 } 108 ierr = VecRestoreArray(bt,&btarray);CHKERRQ(ierr); 109 ierr = MatDenseRestoreArray(B,&Barray);CHKERRQ(ierr); 110 111 /* compute ct = mA^T * cb */ 112 ierr = MatMultTranspose(atb->mA,bt,ct);CHKERRQ(ierr); 113 114 /* transpose local arry of ct to matrix C */ 115 ierr = MatDenseGetArray(C,&Carray);CHKERRQ(ierr); 116 ierr = VecGetArray(ct,&ctarray);CHKERRQ(ierr); 117 k = 0; 118 for (j=0; j<BN; j++) { 119 for (i=0; i<n; i++) Carray[k++] = ctarray[i*BN + j]; 120 } 121 ierr = VecRestoreArray(ct,&ctarray);CHKERRQ(ierr); 122 ierr = MatDenseRestoreArray(C,&Carray);CHKERRQ(ierr); 123 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 124 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 125 PetscFunctionReturn(0); 126 } 127