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 PetscInt i,j,k,m,n,BN=B->cmap->N; 33 Mat_MatTransMatMult *atb; 34 Mat Cdense; 35 PetscScalar *Barray,*Carray,*btarray,*ctarray; 36 Vec bt,ct; 37 Mat_SeqDense *c; 38 39 PetscFunctionBegin; 40 ierr = PetscNew(Mat_MatTransMatMult,&atb);CHKERRQ(ierr); 41 42 /* create MAIJ matrix mA from A */ 43 ierr = MatCreateMAIJ(A,BN,&atb->mA);CHKERRQ(ierr); 44 45 /* create output dense matrix C = A^T*B */ 46 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 47 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cdense);CHKERRQ(ierr); 48 ierr = MatSetSizes(Cdense,n,BN,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 49 ierr = MatSetType(Cdense,MATDENSE);CHKERRQ(ierr); 50 ierr = MatSeqDenseSetPreallocation(Cdense,NULL);CHKERRQ(ierr); 51 52 /* create vectors bt and ct to hold locally transposed arrays of B and C */ 53 ierr = VecCreate(PETSC_COMM_SELF,&bt);CHKERRQ(ierr); 54 ierr = VecSetSizes(bt,m*BN,m*BN);CHKERRQ(ierr); 55 ierr = VecSetFromOptions(bt);CHKERRQ(ierr); 56 ierr = VecCreate(PETSC_COMM_SELF,&ct);CHKERRQ(ierr); 57 ierr = VecSetSizes(ct,n*BN,n*BN);CHKERRQ(ierr); 58 ierr = VecSetFromOptions(ct);CHKERRQ(ierr); 59 atb->bt = bt; 60 atb->ct = ct; 61 62 /* transpose local arry of B, then copy it to vector bt */ 63 ierr = MatDenseGetArray(B,&Barray);CHKERRQ(ierr); 64 ierr = VecGetArray(bt,&btarray);CHKERRQ(ierr); 65 66 k=0; 67 for (j=0; j<BN; j++) { 68 for (i=0; i<m; i++) btarray[i*BN + j] = Barray[k++]; 69 } 70 ierr = VecRestoreArray(bt,&btarray);CHKERRQ(ierr); 71 ierr = MatDenseRestoreArray(B,&Barray);CHKERRQ(ierr); 72 73 /* compute ct = mA^T * cb */ 74 ierr = MatMultTranspose(atb->mA,bt,ct);CHKERRQ(ierr); 75 76 /* transpose local arry of ct to matrix C */ 77 ierr = MatDenseGetArray(Cdense,&Carray);CHKERRQ(ierr); 78 ierr = VecGetArray(ct,&ctarray);CHKERRQ(ierr); 79 k = 0; 80 for (j=0; j<BN; j++) { 81 for (i=0; i<n; i++) Carray[k++] = ctarray[i*BN + j]; 82 } 83 ierr = VecRestoreArray(ct,&ctarray);CHKERRQ(ierr); 84 ierr = MatDenseRestoreArray(Cdense,&Carray);CHKERRQ(ierr); 85 86 *C = Cdense; 87 c = (Mat_SeqDense*)Cdense->data; 88 c->atb = atb; 89 atb->destroy = Cdense->ops->destroy; 90 Cdense->ops->destroy = MatDestroy_SeqDense_MatTransMatMult; 91 PetscFunctionReturn(0); 92 } 93 94 #undef __FUNCT__ 95 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqAIJ_SeqDense" 96 PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 97 { 98 //PetscErrorCode ierr; 99 100 PetscFunctionBegin; 101 PetscFunctionReturn(0); 102 } 103 104 #undef __FUNCT__ 105 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqAIJ_SeqDense" 106 PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C) 107 { 108 //PetscErrorCode ierr; 109 110 PetscFunctionBegin; 111 PetscFunctionReturn(0); 112 } 113