xref: /petsc/src/mat/impls/aij/mpi/mpimatmatmult.c (revision 4ae313f42201736c22765e1d4ea08232bbac4af6) 
1 #define PETSCMAT_DLL
2 
3 /*
4   Defines matrix-matrix product routines for pairs of MPIAIJ matrices
5           C = A * B
6 */
7 #include "src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/
8 #include "src/mat/utils/freespace.h"
9 #include "src/mat/impls/aij/mpi/mpiaij.h"
10 #include "petscbt.h"
11 #include "src/mat/impls/dense/mpi/mpidense.h"
12 
13 #undef __FUNCT__
14 #define __FUNCT__ "MatMatMult_MPIAIJ_MPIAIJ"
15 PetscErrorCode MatMatMult_MPIAIJ_MPIAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill, Mat *C)
16 {
17   PetscErrorCode ierr;
18 
19   PetscFunctionBegin;
20   if (scall == MAT_INITIAL_MATRIX){
21     ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,B,fill,C);CHKERRQ(ierr);/* numeric product is computed as well */
22   } else if (scall == MAT_REUSE_MATRIX){
23     ierr = MatMatMultNumeric_MPIAIJ_MPIAIJ(A,B,*C);CHKERRQ(ierr);
24   } else {
25     SETERRQ1(PETSC_ERR_ARG_WRONG,"Invalid MatReuse %d",(int)scall);
26   }
27   PetscFunctionReturn(0);
28 }
29 
30 #undef __FUNCT__
31 #define __FUNCT__ "PetscObjectContainerDestroy_Mat_MatMatMultMPI"
32 PetscErrorCode PetscObjectContainerDestroy_Mat_MatMatMultMPI(void *ptr)
33 {
34   PetscErrorCode       ierr;
35   Mat_MatMatMultMPI    *mult=(Mat_MatMatMultMPI*)ptr;
36 
37   PetscFunctionBegin;
38   ierr = PetscFree(mult->startsj);CHKERRQ(ierr);
39   ierr = PetscFree(mult->bufa);CHKERRQ(ierr);
40   if (mult->isrowa){ierr = ISDestroy(mult->isrowa);CHKERRQ(ierr);}
41   if (mult->isrowb){ierr = ISDestroy(mult->isrowb);CHKERRQ(ierr);}
42   if (mult->iscolb){ierr = ISDestroy(mult->iscolb);CHKERRQ(ierr);}
43   if (mult->C_seq){ierr = MatDestroy(mult->C_seq);CHKERRQ(ierr);}
44   if (mult->A_loc){ierr = MatDestroy(mult->A_loc);CHKERRQ(ierr); }
45   if (mult->B_seq){ierr = MatDestroy(mult->B_seq);CHKERRQ(ierr);}
46   if (mult->B_loc){ierr = MatDestroy(mult->B_loc);CHKERRQ(ierr);}
47   if (mult->B_oth){ierr = MatDestroy(mult->B_oth);CHKERRQ(ierr);}
48   ierr = PetscFree(mult->abi);CHKERRQ(ierr);
49   ierr = PetscFree(mult->abj);CHKERRQ(ierr);
50   ierr = PetscFree(mult);CHKERRQ(ierr);
51   PetscFunctionReturn(0);
52 }
53 
54 EXTERN PetscErrorCode MatDestroy_AIJ(Mat);
55 #undef __FUNCT__
56 #define __FUNCT__ "MatDestroy_MPIAIJ_MatMatMult"
57 PetscErrorCode MatDestroy_MPIAIJ_MatMatMult(Mat A)
58 {
59   PetscErrorCode       ierr;
60   PetscObjectContainer container;
61   Mat_MatMatMultMPI    *mult=PETSC_NULL;
62 
63   PetscFunctionBegin;
64   ierr = PetscObjectQuery((PetscObject)A,"Mat_MatMatMultMPI",(PetscObject *)&container);CHKERRQ(ierr);
65   if (container) {
66     ierr = PetscObjectContainerGetPointer(container,(void **)&mult);CHKERRQ(ierr);
67   } else {
68     SETERRQ(PETSC_ERR_PLIB,"Container does not exit");
69   }
70   A->ops->destroy = mult->MatDestroy;
71   ierr = PetscObjectCompose((PetscObject)A,"Mat_MatMatMultMPI",0);CHKERRQ(ierr);
72   ierr = (*A->ops->destroy)(A);CHKERRQ(ierr);
73   ierr = PetscObjectContainerDestroy(container);CHKERRQ(ierr);
74   PetscFunctionReturn(0);
75 }
76 
77 #undef __FUNCT__
78 #define __FUNCT__ "MatDuplicate_MPIAIJ_MatMatMult"
79 PetscErrorCode MatDuplicate_MPIAIJ_MatMatMult(Mat A, MatDuplicateOption op, Mat *M) {
80   PetscErrorCode       ierr;
81   Mat_MatMatMultMPI    *mult;
82   PetscObjectContainer container;
83 
84   PetscFunctionBegin;
85   ierr = PetscObjectQuery((PetscObject)A,"Mat_MatMatMultMPI",(PetscObject *)&container);CHKERRQ(ierr);
86   if (container) {
87     ierr  = PetscObjectContainerGetPointer(container,(void **)&mult);CHKERRQ(ierr);
88   } else {
89     SETERRQ(PETSC_ERR_PLIB,"Container does not exit");
90   }
91   ierr = (*mult->MatDuplicate)(A,op,M);CHKERRQ(ierr);
92   (*M)->ops->destroy   = mult->MatDestroy;   /* =MatDestroy_MPIAIJ, *M doesn't duplicate A's container! */
93   (*M)->ops->duplicate = mult->MatDuplicate; /* =MatDuplicate_ MPIAIJ */
94   PetscFunctionReturn(0);
95 }
96 
97 #undef __FUNCT__
98 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIAIJ"
99 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
100 {
101   PetscErrorCode       ierr;
102   PetscInt             start,end;
103   Mat_MatMatMultMPI    *mult;
104   PetscObjectContainer container;
105 
106   PetscFunctionBegin;
107   if (A->cmap.rstart != B->rmap.rstart || A->cmap.rend != B->rmap.rend){
108     SETERRQ4(PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, (%D, %D) != (%D,%D)",A->cmap.rstart,A->cmap.rend,B->rmap.rstart,B->rmap.rend);
109   }
110   ierr = PetscNew(Mat_MatMatMultMPI,&mult);CHKERRQ(ierr);
111 
112   /* create a seq matrix B_seq = submatrix of B by taking rows of B that equal to nonzero col of A */
113   ierr = MatGetBrowsOfAcols(A,B,MAT_INITIAL_MATRIX,&mult->isrowb,&mult->iscolb,&mult->brstart,&mult->B_seq);CHKERRQ(ierr);
114 
115   /*  create a seq matrix A_seq = submatrix of A by taking all local rows of A */
116   start = A->rmap.rstart; end = A->rmap.rend;
117   ierr = ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&mult->isrowa);CHKERRQ(ierr);
118   ierr = MatGetLocalMatCondensed(A,MAT_INITIAL_MATRIX,&mult->isrowa,&mult->isrowb,&mult->A_loc);CHKERRQ(ierr);
119 
120   /* compute C_seq = A_seq * B_seq */
121   ierr = MatMatMult_SeqAIJ_SeqAIJ(mult->A_loc,mult->B_seq,MAT_INITIAL_MATRIX,fill,&mult->C_seq);CHKERRQ(ierr);
122 
123   /* create mpi matrix C by concatinating C_seq */
124   ierr = PetscObjectReference((PetscObject)mult->C_seq);CHKERRQ(ierr); /* prevent C_seq being destroyed by MatMerge() */
125   ierr = MatMerge(A->comm,mult->C_seq,B->cmap.n,MAT_INITIAL_MATRIX,C);CHKERRQ(ierr);
126 
127   /* attach the supporting struct to C for reuse of symbolic C */
128   ierr = PetscObjectContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr);
129   ierr = PetscObjectContainerSetPointer(container,mult);CHKERRQ(ierr);
130   ierr = PetscObjectCompose((PetscObject)(*C),"Mat_MatMatMultMPI",(PetscObject)container);CHKERRQ(ierr);
131   ierr = PetscObjectContainerSetUserDestroy(container,PetscObjectContainerDestroy_Mat_MatMatMultMPI);CHKERRQ(ierr);
132   mult->MatDestroy   = (*C)->ops->destroy;
133   mult->MatDuplicate = (*C)->ops->duplicate;
134 
135   (*C)->ops->destroy   = MatDestroy_MPIAIJ_MatMatMult;
136   (*C)->ops->duplicate = MatDuplicate_MPIAIJ_MatMatMult;
137   PetscFunctionReturn(0);
138 }
139 
140 /* This routine is called ONLY in the case of reusing previously computed symbolic C */
141 #undef __FUNCT__
142 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIAIJ"
143 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C)
144 {
145   PetscErrorCode       ierr;
146   Mat                  *seq;
147   Mat_MatMatMultMPI    *mult;
148   PetscObjectContainer container;
149 
150   PetscFunctionBegin;
151   ierr = PetscObjectQuery((PetscObject)C,"Mat_MatMatMultMPI",(PetscObject *)&container);CHKERRQ(ierr);
152   if (container) {
153     ierr  = PetscObjectContainerGetPointer(container,(void **)&mult);CHKERRQ(ierr);
154   } else {
155     SETERRQ(PETSC_ERR_PLIB,"Container does not exit");
156   }
157 
158   seq = &mult->B_seq;
159   ierr = MatGetSubMatrices(B,1,&mult->isrowb,&mult->iscolb,MAT_REUSE_MATRIX,&seq);CHKERRQ(ierr);
160   mult->B_seq = *seq;
161 
162   seq = &mult->A_loc;
163   ierr = MatGetSubMatrices(A,1,&mult->isrowa,&mult->isrowb,MAT_REUSE_MATRIX,&seq);CHKERRQ(ierr);
164   mult->A_loc = *seq;
165 
166   ierr = MatMatMult_SeqAIJ_SeqAIJ(mult->A_loc,mult->B_seq,MAT_REUSE_MATRIX,0.0,&mult->C_seq);CHKERRQ(ierr);
167 
168   ierr = PetscObjectReference((PetscObject)mult->C_seq);CHKERRQ(ierr);
169   ierr = MatMerge(A->comm,mult->C_seq,B->cmap.n,MAT_REUSE_MATRIX,&C);CHKERRQ(ierr);
170   PetscFunctionReturn(0);
171 }
172 
173 #undef __FUNCT__
174 #define __FUNCT__ "MatMatMult_MPIAIJ_MPIDense"
175 PetscErrorCode MatMatMult_MPIAIJ_MPIDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
176 {
177   PetscErrorCode ierr;
178 
179   PetscFunctionBegin;
180   if (scall == MAT_INITIAL_MATRIX){
181     ierr = MatMatMultSymbolic_MPIAIJ_MPIDense(A,B,fill,C);CHKERRQ(ierr);
182   }
183   ierr = MatMatMultNumeric_MPIAIJ_MPIDense(A,B,*C);CHKERRQ(ierr);
184   PetscFunctionReturn(0);
185 }
186 
187 #undef __FUNCT__
188 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIDense"
189 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIDense(Mat A,Mat B,PetscReal fill,Mat *C)
190 {
191   PetscErrorCode ierr;
192 
193   PetscFunctionBegin;
194   ierr = MatMatMultSymbolic_MPIDense_MPIDense(A,B,0.0,C);
195   PetscFunctionReturn(0);
196 }
197 
198 #undef __FUNCT__
199 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIDense"
200 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIDense(Mat A,Mat B,Mat C)
201 {
202   PetscErrorCode ierr;
203   PetscScalar    *b,*c,*barray,*carray;
204   PetscInt       cm=C->rmap.n, bN=B->cmap.N, bm=B->rmap.n, col;
205   Vec            vb,vc;
206 
207   PetscFunctionBegin;
208   if (!C->rmap.N || !bN) PetscFunctionReturn(0);
209   ierr = VecCreateMPI(B->comm,bm,PETSC_DECIDE,&vb);CHKERRQ(ierr);
210   ierr = VecCreateMPI(A->comm,cm,PETSC_DECIDE,&vc);CHKERRQ(ierr);
211 
212   ierr = MatGetArray(B,&barray);CHKERRQ(ierr);
213   ierr = MatGetArray(C,&carray);CHKERRQ(ierr);
214   b = barray; c = carray;
215   for (col=0; col<bN; col++){
216     ierr = VecPlaceArray(vb,b);CHKERRQ(ierr);
217     ierr = VecPlaceArray(vc,c);CHKERRQ(ierr);
218     ierr = MatMult(A,vb,vc);CHKERRQ(ierr);
219     b += bm; c += cm;
220     ierr = VecResetArray(vb);CHKERRQ(ierr);
221     ierr = VecResetArray(vc);CHKERRQ(ierr);
222   }
223   ierr = MatRestoreArray(B,&barray);CHKERRQ(ierr);
224   ierr = MatRestoreArray(C,&carray);CHKERRQ(ierr);
225   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
226   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
227   ierr = VecDestroy(vb);CHKERRQ(ierr);
228   ierr = VecDestroy(vc);CHKERRQ(ierr);
229   PetscFunctionReturn(0);
230 }
231 
232