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