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