Actual source code: mattransposematmult.c
petsc-3.13.6 2020-09-29
2: /*
3: Defines matrix-matrix product routines
4: C = A^T * B
5: */
7: #include <../src/mat/impls/aij/seq/aij.h>
8: #include <../src/mat/impls/dense/seq/dense.h>
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;
17: MatDestroy(&atb->mA);
18: VecDestroy(&atb->bt);
19: VecDestroy(&atb->ct);
20: (atb->destroy)(A);
21: PetscFree(atb);
22: return(0);
23: }
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;
34: PetscNew(&atb);
36: /* create output dense matrix C = A^T*B */
37: MatSetSizes(C,n,BN,n,BN);
38: MatSetType(C,MATSEQDENSE);
39: MatSetUp(C);
41: /* create vectors bt and ct to hold locally transposed arrays of B and C */
42: VecCreate(PETSC_COMM_SELF,&bt);
43: VecSetSizes(bt,m*BN,m*BN);
44: VecSetType(bt,VECSTANDARD);
45: VecCreate(PETSC_COMM_SELF,&ct);
46: VecSetSizes(ct,n*BN,n*BN);
47: VecSetType(ct,VECSTANDARD);
48: atb->bt = bt;
49: atb->ct = ct;
51: c = (Mat_SeqDense*)C->data;
52: c->atb = atb;
53: atb->destroy = C->ops->destroy;
54: C->ops->destroy = MatDestroy_SeqDense_MatTransMatMult;
55: C->ops->transposematmultnumeric = MatTransposeMatMultNumeric_SeqAIJ_SeqDense;
56: return(0);
57: }
59: PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
60: {
61: PetscErrorCode ierr;
62: PetscInt i,j,k,m=A->rmap->n,n=A->cmap->n,BN=B->cmap->N;
63: const PetscScalar *Barray,*ctarray;
64: PetscScalar *Carray,*btarray;
65: Mat_SeqDense *c=(Mat_SeqDense*)C->data;
66: Mat_MatTransMatMult *atb=c->atb;
67: Vec bt=atb->bt,ct=atb->ct;
70: /* create MAIJ matrix mA from A -- should be done in symbolic phase */
71: MatDestroy(&atb->mA);
72: MatCreateMAIJ(A,BN,&atb->mA);
74: /* transpose local arry of B, then copy it to vector bt */
75: MatDenseGetArrayRead(B,&Barray);
76: VecGetArray(bt,&btarray);
78: k=0;
79: for (j=0; j<BN; j++) {
80: for (i=0; i<m; i++) btarray[i*BN + j] = Barray[k++];
81: }
82: VecRestoreArray(bt,&btarray);
83: MatDenseRestoreArrayRead(B,&Barray);
84:
85: /* compute ct = mA^T * cb */
86: MatMultTranspose(atb->mA,bt,ct);
88: /* transpose local arry of ct to matrix C */
89: MatDenseGetArray(C,&Carray);
90: VecGetArrayRead(ct,&ctarray);
91: k = 0;
92: for (j=0; j<BN; j++) {
93: for (i=0; i<n; i++) Carray[k++] = ctarray[i*BN + j];
94: }
95: VecRestoreArrayRead(ct,&ctarray);
96: MatDenseRestoreArray(C,&Carray);
97: return(0);
98: }