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 PetscErrorCode MatDestroy_SeqDense_MatTransMatMult(Mat A) 11 { 12 PetscErrorCode ierr; 13 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 14 Mat_MatTransMatMult *atb = a->atb; 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 = (atb->destroy)(A);CHKERRQ(ierr); 21 ierr = PetscFree(atb);CHKERRQ(ierr); 22 PetscFunctionReturn(0); 23 } 24 25 PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 26 { 27 PetscErrorCode ierr; 28 PetscInt m=A->rmap->n,n=A->cmap->n,BN=B->cmap->N; 29 Mat_MatTransMatMult *atb; 30 Vec bt,ct; 31 Mat_SeqDense *c; 32 PetscBool cisdense; 33 34 PetscFunctionBegin; 35 ierr = PetscNew(&atb);CHKERRQ(ierr); 36 37 /* create output dense matrix C = A^T*B */ 38 ierr = MatSetSizes(C,n,BN,n,BN);CHKERRQ(ierr); 39 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 40 if (!cisdense) { 41 ierr = MatSetType(C,((PetscObject)B)->type_name);CHKERRQ(ierr); 42 } 43 ierr = MatSetUp(C);CHKERRQ(ierr); 44 45 /* create vectors bt and ct to hold locally transposed arrays of B and C */ 46 ierr = VecCreate(PETSC_COMM_SELF,&bt);CHKERRQ(ierr); 47 ierr = VecSetSizes(bt,m*BN,m*BN);CHKERRQ(ierr); 48 ierr = VecSetType(bt,A->defaultvectype);CHKERRQ(ierr); 49 ierr = VecCreate(PETSC_COMM_SELF,&ct);CHKERRQ(ierr); 50 ierr = VecSetSizes(ct,n*BN,n*BN);CHKERRQ(ierr); 51 ierr = VecSetType(ct,A->defaultvectype);CHKERRQ(ierr); 52 atb->bt = bt; 53 atb->ct = ct; 54 55 c = (Mat_SeqDense*)C->data; 56 c->atb = atb; 57 atb->destroy = C->ops->destroy; 58 C->ops->destroy = MatDestroy_SeqDense_MatTransMatMult; 59 C->ops->transposematmultnumeric = MatTransposeMatMultNumeric_SeqAIJ_SeqDense; 60 PetscFunctionReturn(0); 61 } 62 63 PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C) 64 { 65 PetscErrorCode ierr; 66 PetscInt i,j,k,m=A->rmap->n,n=A->cmap->n,BN=B->cmap->N; 67 const PetscScalar *Barray,*ctarray; 68 PetscScalar *Carray,*btarray; 69 Mat_SeqDense *c=(Mat_SeqDense*)C->data; 70 Mat_MatTransMatMult *atb=c->atb; 71 Vec bt=atb->bt,ct=atb->ct; 72 73 PetscFunctionBegin; 74 /* create MAIJ matrix mA from A -- should be done in symbolic phase */ 75 ierr = MatDestroy(&atb->mA);CHKERRQ(ierr); 76 ierr = MatCreateMAIJ(A,BN,&atb->mA);CHKERRQ(ierr); 77 78 /* transpose local array of B, then copy it to vector bt */ 79 ierr = MatDenseGetArrayRead(B,&Barray);CHKERRQ(ierr); 80 ierr = VecGetArray(bt,&btarray);CHKERRQ(ierr); 81 82 k=0; 83 for (j=0; j<BN; j++) { 84 for (i=0; i<m; i++) btarray[i*BN + j] = Barray[k++]; 85 } 86 ierr = VecRestoreArray(bt,&btarray);CHKERRQ(ierr); 87 ierr = MatDenseRestoreArrayRead(B,&Barray);CHKERRQ(ierr); 88 89 /* compute ct = mA^T * cb */ 90 ierr = MatMultTranspose(atb->mA,bt,ct);CHKERRQ(ierr); 91 92 /* transpose local array of ct to matrix C */ 93 ierr = MatDenseGetArray(C,&Carray);CHKERRQ(ierr); 94 ierr = VecGetArrayRead(ct,&ctarray);CHKERRQ(ierr); 95 k = 0; 96 for (j=0; j<BN; j++) { 97 for (i=0; i<n; i++) Carray[k++] = ctarray[i*BN + j]; 98 } 99 ierr = VecRestoreArrayRead(ct,&ctarray);CHKERRQ(ierr); 100 ierr = MatDenseRestoreArray(C,&Carray);CHKERRQ(ierr); 101 PetscFunctionReturn(0); 102 } 103