1
2 static char help[] = "Test Matrix products for AIJ matrices\n\
3 Input arguments are:\n\
4 -fA <input_file> -fB <input_file> -fC <input_file>: file to load\n\n";
5 /* Example of usage:
6 ./ex62 -fA <A_binary> -fB <B_binary>
7 mpiexec -n 3 ./ex62 -fA medium -fB medium
8 */
9
10 #include <petscmat.h>
11
12 /*
13 B = A - B
14 norm = norm(B)
15 */
MatNormDifference(Mat A,Mat B,PetscReal * norm)16 PetscErrorCode MatNormDifference(Mat A,Mat B,PetscReal *norm)
17 {
18 PetscErrorCode ierr;
19
20 PetscFunctionBegin;
21 ierr = MatAXPY(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
22 ierr = MatNorm(B,NORM_FROBENIUS,norm);CHKERRQ(ierr);
23 PetscFunctionReturn(0);
24 }
25
main(int argc,char ** args)26 int main(int argc,char **args)
27 {
28 Mat A,A_save,B,C,P,C1,R;
29 PetscViewer viewer;
30 PetscErrorCode ierr;
31 PetscMPIInt size,rank;
32 PetscInt i,j,*idxn,M,N,nzp,PN,rstart,rend;
33 PetscReal norm;
34 PetscRandom rdm;
35 char file[2][128];
36 PetscScalar *a,rval,alpha;
37 PetscBool Test_MatMatMult=PETSC_TRUE,Test_MatTrMat=PETSC_TRUE,Test_MatMatTr=PETSC_TRUE;
38 PetscBool Test_MatPtAP=PETSC_TRUE,Test_MatRARt=PETSC_TRUE,flg,seqaij;
39 MatInfo info;
40 MatType mattype;
41
42 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
43 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
44 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
45
46 /* Load the matrices A_save and B */
47 ierr = PetscOptionsGetString(NULL,NULL,"-fA",file[0],sizeof(file[0]),&flg);CHKERRQ(ierr);
48 if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix A with the -fA option.");
49 ierr = PetscOptionsGetString(NULL,NULL,"-fB",file[1],sizeof(file[1]),&flg);CHKERRQ(ierr);
50 if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix B with the -fB option.");
51
52 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);CHKERRQ(ierr);
53 ierr = MatCreate(PETSC_COMM_WORLD,&A_save);CHKERRQ(ierr);
54 ierr = MatSetFromOptions(A_save);CHKERRQ(ierr);
55 ierr = MatLoad(A_save,viewer);CHKERRQ(ierr);
56 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
57
58 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);CHKERRQ(ierr);
59 ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
60 ierr = MatSetFromOptions(B);CHKERRQ(ierr);
61 ierr = MatLoad(B,viewer);CHKERRQ(ierr);
62 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
63
64 ierr = MatGetType(B,&mattype);CHKERRQ(ierr);
65
66 ierr = MatGetSize(B,&M,&N);CHKERRQ(ierr);
67 nzp = PetscMax((PetscInt)(0.1*M),5);
68 ierr = PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);CHKERRQ(ierr);
69 a = (PetscScalar*)(idxn + nzp);
70
71 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rdm);CHKERRQ(ierr);
72 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);
73
74 /* 1) MatMatMult() */
75 /* ----------------*/
76 if (Test_MatMatMult) {
77 ierr = MatDuplicate(A_save,MAT_COPY_VALUES,&A);CHKERRQ(ierr);
78
79 /* (1.1) Test developer API */
80 ierr = MatProductCreate(A,B,NULL,&C);CHKERRQ(ierr);
81 ierr = MatSetOptionsPrefix(C,"AB_");CHKERRQ(ierr);
82 ierr = MatProductSetType(C,MATPRODUCT_AB);CHKERRQ(ierr);
83 ierr = MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
84 ierr = MatProductSetFill(C,PETSC_DEFAULT);CHKERRQ(ierr);
85 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr);
86 /* we can inquire about MATOP_PRODUCTSYMBOLIC even if the destination matrix type has not been set yet */
87 ierr = MatHasOperation(C,MATOP_PRODUCTSYMBOLIC,&flg);CHKERRQ(ierr);
88 ierr = MatProductSymbolic(C);CHKERRQ(ierr);
89 ierr = MatProductNumeric(C);CHKERRQ(ierr);
90
91 /* Test reuse symbolic C */
92 alpha = 0.9;
93 ierr = MatScale(A,alpha);CHKERRQ(ierr);
94 ierr = MatProductNumeric(C);CHKERRQ(ierr);
95
96 ierr = MatMatMultEqual(A,B,C,10,&flg);CHKERRQ(ierr);
97 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B");
98 ierr = MatDestroy(&C);CHKERRQ(ierr);
99
100 /* (1.2) Test user driver */
101 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
102
103 /* Test MAT_REUSE_MATRIX - reuse symbolic C */
104 alpha = 1.0;
105 for (i=0; i<2; i++) {
106 alpha -= 0.1;
107 ierr = MatScale(A,alpha);CHKERRQ(ierr);
108 ierr = MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
109 }
110 ierr = MatMatMultEqual(A,B,C,10,&flg);CHKERRQ(ierr);
111 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()");
112 ierr = MatDestroy(&A);CHKERRQ(ierr);
113
114 /* Test MatProductClear() */
115 ierr = MatProductClear(C);CHKERRQ(ierr);
116 ierr = MatDestroy(&C);CHKERRQ(ierr);
117
118 /* Test MatMatMult() for dense and aij matrices */
119 ierr = MatConvert(A_save,MATDENSE,MAT_INITIAL_MATRIX,&A);CHKERRQ(ierr);
120 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
121 ierr = MatDestroy(&C);CHKERRQ(ierr);
122 ierr = MatDestroy(&A);CHKERRQ(ierr);
123
124 }
125
126 /* Create P and R = P^T */
127 /* --------------------- */
128 PN = M/2;
129 nzp = 5; /* num of nonzeros in each row of P */
130 ierr = MatCreate(PETSC_COMM_WORLD,&P);CHKERRQ(ierr);
131 ierr = MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);CHKERRQ(ierr);
132 ierr = MatSetType(P,mattype);CHKERRQ(ierr);
133 ierr = MatSeqAIJSetPreallocation(P,nzp,NULL);CHKERRQ(ierr);
134 ierr = MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);CHKERRQ(ierr);
135 ierr = MatGetOwnershipRange(P,&rstart,&rend);CHKERRQ(ierr);
136 for (i=0; i<nzp; i++) {
137 ierr = PetscRandomGetValue(rdm,&a[i]);CHKERRQ(ierr);
138 }
139 for (i=rstart; i<rend; i++) {
140 for (j=0; j<nzp; j++) {
141 ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
142 idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
143 }
144 ierr = MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);CHKERRQ(ierr);
145 }
146 ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
147 ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
148
149 ierr = MatTranspose(P,MAT_INITIAL_MATRIX,&R);CHKERRQ(ierr);
150
151 /* 2) MatTransposeMatMult() */
152 /* ------------------------ */
153 if (Test_MatTrMat) {
154 /* (2.1) Test developer driver C = P^T*B */
155 ierr = MatProductCreate(P,B,NULL,&C);CHKERRQ(ierr);
156 ierr = MatSetOptionsPrefix(C,"AtB_");CHKERRQ(ierr);
157 ierr = MatProductSetType(C,MATPRODUCT_AtB);CHKERRQ(ierr);
158 ierr = MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
159 ierr = MatProductSetFill(C,PETSC_DEFAULT);CHKERRQ(ierr);
160 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr);
161 ierr = MatProductSymbolic(C);CHKERRQ(ierr); /* equivalent to MatSetUp() */
162 ierr = MatSetOption(C,MAT_USE_INODES,PETSC_FALSE);CHKERRQ(ierr); /* illustrate how to call MatSetOption() */
163 ierr = MatProductNumeric(C);CHKERRQ(ierr);
164 ierr = MatProductNumeric(C);CHKERRQ(ierr); /* test reuse symbolic C */
165
166 ierr = MatTransposeMatMultEqual(P,B,C,10,&flg);CHKERRQ(ierr);
167 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error: developer driver C = P^T*B");
168 ierr = MatDestroy(&C);CHKERRQ(ierr);
169
170 /* (2.2) Test user driver C = P^T*B */
171 ierr = MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
172 ierr = MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
173 ierr = MatGetInfo(C,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr);
174 ierr = MatProductClear(C);CHKERRQ(ierr);
175
176 /* Compare P^T*B and R*B */
177 ierr = MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1);CHKERRQ(ierr);
178 ierr = MatNormDifference(C,C1,&norm);CHKERRQ(ierr);
179 if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatTransposeMatMult(): %g",(double)norm);
180 ierr = MatDestroy(&C1);CHKERRQ(ierr);
181
182 /* Test MatDuplicate() of C=P^T*B */
183 ierr = MatDuplicate(C,MAT_COPY_VALUES,&C1);CHKERRQ(ierr);
184 ierr = MatDestroy(&C1);CHKERRQ(ierr);
185 ierr = MatDestroy(&C);CHKERRQ(ierr);
186 }
187
188 /* 3) MatMatTransposeMult() */
189 /* ------------------------ */
190 if (Test_MatMatTr) {
191 /* C = B*R^T */
192 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr);
193 if (size == 1 && seqaij) {
194 ierr = MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
195 ierr = MatSetOptionsPrefix(C,"ABt_");CHKERRQ(ierr); /* enable '-ABt_' for matrix C */
196 ierr = MatGetInfo(C,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr);
197
198 /* Test MAT_REUSE_MATRIX - reuse symbolic C */
199 ierr = MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
200
201 /* Check */
202 ierr = MatMatMult(B,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1);CHKERRQ(ierr);
203 ierr = MatNormDifference(C,C1,&norm);CHKERRQ(ierr);
204 if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatTransposeMult() %g",(double)norm);
205 ierr = MatDestroy(&C1);CHKERRQ(ierr);
206 ierr = MatDestroy(&C);CHKERRQ(ierr);
207 }
208 }
209
210 /* 4) Test MatPtAP() */
211 /*-------------------*/
212 if (Test_MatPtAP) {
213 ierr = MatDuplicate(A_save,MAT_COPY_VALUES,&A);CHKERRQ(ierr);
214
215 /* (4.1) Test developer API */
216 ierr = MatProductCreate(A,P,NULL,&C);CHKERRQ(ierr);
217 ierr = MatSetOptionsPrefix(C,"PtAP_");CHKERRQ(ierr);
218 ierr = MatProductSetType(C,MATPRODUCT_PtAP);CHKERRQ(ierr);
219 ierr = MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
220 ierr = MatProductSetFill(C,PETSC_DEFAULT);CHKERRQ(ierr);
221 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr);
222 ierr = MatProductSymbolic(C);CHKERRQ(ierr);
223 ierr = MatProductNumeric(C);CHKERRQ(ierr);
224 ierr = MatProductNumeric(C);CHKERRQ(ierr); /* reuse symbolic C */
225
226 ierr = MatPtAPMultEqual(A,P,C,10,&flg);CHKERRQ(ierr);
227 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
228 ierr = MatDestroy(&C);CHKERRQ(ierr);
229
230 /* (4.2) Test user driver */
231 ierr = MatPtAP(A,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
232
233 /* Test MAT_REUSE_MATRIX - reuse symbolic C */
234 alpha=1.0;
235 for (i=0; i<2; i++) {
236 alpha -= 0.1;
237 ierr = MatScale(A,alpha);CHKERRQ(ierr);
238 ierr = MatPtAP(A,P,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr);
239 }
240
241 ierr = MatPtAPMultEqual(A,P,C,10,&flg);CHKERRQ(ierr);
242 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP");
243
244 /* 5) Test MatRARt() */
245 /* ----------------- */
246 if (Test_MatRARt) {
247 Mat RARt;
248 ierr = MatTranspose(P,MAT_REUSE_MATRIX,&R);CHKERRQ(ierr);
249
250 /* (5.1) Test developer driver RARt = R*A*Rt */
251 ierr = MatProductCreate(A,R,NULL,&RARt);CHKERRQ(ierr);
252 ierr = MatSetOptionsPrefix(RARt,"RARt_");CHKERRQ(ierr);
253 ierr = MatProductSetType(RARt,MATPRODUCT_RARt);CHKERRQ(ierr);
254 ierr = MatProductSetAlgorithm(RARt,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
255 ierr = MatProductSetFill(RARt,PETSC_DEFAULT);CHKERRQ(ierr);
256 ierr = MatProductSetFromOptions(RARt);CHKERRQ(ierr);
257 ierr = MatProductSymbolic(RARt);CHKERRQ(ierr); /* equivalent to MatSetUp() */
258 ierr = MatSetOption(RARt,MAT_USE_INODES,PETSC_FALSE);CHKERRQ(ierr); /* illustrate how to call MatSetOption() */
259 ierr = MatProductNumeric(RARt);CHKERRQ(ierr);
260 ierr = MatProductNumeric(RARt);CHKERRQ(ierr); /* test reuse symbolic RARt */
261 ierr = MatDestroy(&RARt);CHKERRQ(ierr);
262
263 /* (2.2) Test user driver RARt = R*A*Rt */
264 ierr = MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);CHKERRQ(ierr);
265 ierr = MatRARt(A,R,MAT_REUSE_MATRIX,2.0,&RARt);CHKERRQ(ierr);
266
267 ierr = MatNormDifference(C,RARt,&norm);CHKERRQ(ierr);
268 if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"|PtAP - RARt| = %g",(double)norm);
269 ierr = MatDestroy(&RARt);CHKERRQ(ierr);
270 }
271
272 ierr = MatDestroy(&A);CHKERRQ(ierr);
273 ierr = MatDestroy(&C);CHKERRQ(ierr);
274 }
275
276 /* Destroy objects */
277 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr);
278 ierr = PetscFree(idxn);CHKERRQ(ierr);
279
280 ierr = MatDestroy(&A_save);CHKERRQ(ierr);
281 ierr = MatDestroy(&B);CHKERRQ(ierr);
282 ierr = MatDestroy(&P);CHKERRQ(ierr);
283 ierr = MatDestroy(&R);CHKERRQ(ierr);
284
285 PetscFinalize();
286 return ierr;
287 }
288
289 /*TEST
290 test:
291 suffix: 1
292 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
293 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
294 output_file: output/ex62_1.out
295
296 test:
297 suffix: 2_ab_scalable
298 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
299 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via scalable -matmatmult_via scalable -AtB_matproduct_atb_via outerproduct -mattransposematmult_via outerproduct
300 output_file: output/ex62_1.out
301
302 test:
303 suffix: 3_ab_scalable_fast
304 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
305 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via scalable_fast -matmatmult_via scalable_fast -matmattransmult_via color
306 output_file: output/ex62_1.out
307
308 test:
309 suffix: 4_ab_heap
310 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
311 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via heap -matmatmult_via heap -PtAP_matproduct_ptap_via rap -matptap_via rap
312 output_file: output/ex62_1.out
313
314 test:
315 suffix: 5_ab_btheap
316 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
317 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via btheap -matmatmult_via btheap -matrart_via r*art
318 output_file: output/ex62_1.out
319
320 test:
321 suffix: 6_ab_llcondensed
322 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
323 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via llcondensed -matmatmult_via llcondensed -matrart_via coloring_rart
324 output_file: output/ex62_1.out
325
326 test:
327 suffix: 7_ab_rowmerge
328 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
329 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via rowmerge -matmatmult_via rowmerge
330 output_file: output/ex62_1.out
331
332 test:
333 suffix: 8_ab_hypre
334 requires: hypre datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
335 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via hypre -matmatmult_via hypre -PtAP_matproduct_ptap_via hypre -matptap_via hypre
336 output_file: output/ex62_1.out
337
338 test:
339 suffix: 10
340 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
341 nsize: 3
342 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
343 output_file: output/ex62_1.out
344
345 test:
346 suffix: 11_ab_scalable
347 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
348 nsize: 3
349 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via scalable -matmatmult_via scalable -AtB_matproduct_atb_via scalable -mattransposematmult_via scalable
350 output_file: output/ex62_1.out
351
352 test:
353 suffix: 12_ab_seqmpi
354 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
355 nsize: 3
356 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via seqmpi -matmatmult_via seqmpi -AtB_matproduct_atb_via at*b -mattransposematmult_via at*b
357 output_file: output/ex62_1.out
358
359 test:
360 suffix: 13_ab_hypre
361 requires: hypre datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
362 nsize: 3
363 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via hypre -matmatmult_via hypre -PtAP_matproduct_ptap_via hypre -matptap_via hypre
364 output_file: output/ex62_1.out
365
366 TEST*/
367