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 9 #undef __FUNCT__ 10 #define __FUNCT__ "MatTransposeMatMult_SeqAIJ_SeqDense" 11 PetscErrorCode MatTransposeMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 12 { 13 PetscErrorCode ierr; 14 PetscInt i,j,k,m,n,dof=B->cmap->N; 15 Mat_MatTransMatMult *atb; 16 Mat Cdense; 17 PetscScalar *Barray,*Carray,*btarray,*ctarray; 18 Vec bt,ct; 19 20 PetscFunctionBegin; 21 ierr = PetscNew(Mat_MatTransMatMult,&atb);CHKERRQ(ierr); 22 23 /* create MAIJ matrix mA from A */ 24 ierr = MatCreateMAIJ(A,dof,&atb->mA);CHKERRQ(ierr); 25 26 /* create output dense matrix C = A^T*B */ 27 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 28 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cdense);CHKERRQ(ierr); 29 ierr = MatSetSizes(Cdense,n,dof,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 30 ierr = MatSetType(Cdense,MATDENSE);CHKERRQ(ierr); 31 ierr = MatSeqDenseSetPreallocation(Cdense,NULL);CHKERRQ(ierr); 32 //ierr = MatMPIDenseSetPreallocation(Cdense,NULL);CHKERRQ(ierr); 33 34 /* create vectors bt and ct to hold locally transposed arrays of B and C */ 35 ierr = VecCreate(PETSC_COMM_SELF,&bt);CHKERRQ(ierr); 36 ierr = VecSetSizes(bt,m*dof,m*dof);CHKERRQ(ierr); 37 ierr = VecSetFromOptions(bt);CHKERRQ(ierr); 38 ierr = VecCreate(PETSC_COMM_SELF,&ct);CHKERRQ(ierr); 39 ierr = VecSetSizes(ct,n*dof,n*dof);CHKERRQ(ierr); 40 ierr = VecSetFromOptions(ct);CHKERRQ(ierr); 41 atb->bt = bt; 42 atb->ct = ct; 43 44 /* transpose local arry of B, then copy it to vector bt */ 45 ierr = MatDenseGetArray(B,&Barray);CHKERRQ(ierr); 46 ierr = VecGetArray(bt,&btarray);CHKERRQ(ierr); 47 48 k=0; 49 for (j=0; j<dof; j++) { 50 for (i=0; i<m; i++) btarray[i*dof + j] = Barray[k++]; 51 } 52 ierr = VecRestoreArray(bt,&btarray);CHKERRQ(ierr); 53 ierr = MatDenseRestoreArray(B,&Barray);CHKERRQ(ierr); 54 55 /* compute ct = mA^T * cb */ 56 ierr = MatMultTranspose(atb->mA,bt,ct);CHKERRQ(ierr); 57 58 /* transpose local arry of ct to matrix C */ 59 ierr = MatDenseGetArray(Cdense,&Carray);CHKERRQ(ierr); 60 ierr = VecGetArray(ct,&ctarray);CHKERRQ(ierr); 61 k = 0; 62 for (j=0; j<dof; j++) { 63 for (i=0; i<n; i++) Carray[k++] = ctarray[i*dof + j]; 64 } 65 ierr = VecRestoreArray(ct,&ctarray);CHKERRQ(ierr); 66 ierr = MatDenseRestoreArray(Cdense,&Carray);CHKERRQ(ierr); 67 68 *C = Cdense; 69 70 ierr = MatDestroy(&atb->mA);CHKERRQ(ierr); 71 ierr = VecDestroy(&atb->bt);CHKERRQ(ierr); 72 ierr = VecDestroy(&atb->ct);CHKERRQ(ierr); 73 ierr = PetscFree(atb);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqAIJ_SeqDense" 79 PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 80 { 81 //PetscErrorCode ierr; 82 83 PetscFunctionBegin; 84 PetscFunctionReturn(0); 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqAIJ_SeqDense" 89 PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C) 90 { 91 //PetscErrorCode ierr; 92 93 PetscFunctionBegin; 94 PetscFunctionReturn(0); 95 } 96