Actual source code: ex62.c
petsc-3.15.0 2021-03-30
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: */
10: #include <petscmat.h>
12: /*
13: B = A - B
14: norm = norm(B)
15: */
16: PetscErrorCode MatNormDifference(Mat A,Mat B,PetscReal *norm)
17: {
21: MatAXPY(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
22: MatNorm(B,NORM_FROBENIUS,norm);
23: return(0);
24: }
26: int main(int argc,char **args)
27: {
28: Mat A,A_save,B,C,P,C1,R;
29: PetscViewer viewer;
31: PetscMPIInt size,rank;
32: PetscInt i,j,*idxn,PM,PN = PETSC_DECIDE,rstart,rend;
33: PetscReal norm;
34: PetscRandom rdm;
35: char file[2][PETSC_MAX_PATH_LEN] = { "", ""};
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,flgA,flgB;
39: MatInfo info;
40: PetscInt nzp = 5; /* num of nonzeros in each row of P */
41: MatType mattype;
42: const char *deft = MATAIJ;
43: char A_mattype[256], B_mattype[256];
44: PetscInt mcheck = 10;
46: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
47: MPI_Comm_size(PETSC_COMM_WORLD,&size);
48: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
50: /* Load the matrices A_save and B */
51: PetscOptionsBegin(PETSC_COMM_WORLD,"","","");
52: PetscOptionsInt("-PN","Number of columns of P","",PN,&PN,NULL);
53: PetscOptionsInt("-mcheck","Number of matmult checks","",mcheck,&mcheck,NULL);
54: PetscOptionsString("-fA","Path for matrix A","",file[0],file[0],sizeof(file[0]),&flg);
55: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for matrix A with the -fA option.");
56: PetscOptionsString("-fB","Path for matrix B","",file[1],file[1],sizeof(file[1]),&flg);
57: PetscOptionsFList("-A_mat_type","Matrix type","MatSetType",MatList,deft,A_mattype,256,&flgA);
58: PetscOptionsFList("-B_mat_type","Matrix type","MatSetType",MatList,deft,B_mattype,256,&flgB);
59: PetscOptionsEnd();
61: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);
62: MatCreate(PETSC_COMM_WORLD,&A_save);
63: MatLoad(A_save,viewer);
64: PetscViewerDestroy(&viewer);
66: if (flg) {
67: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);
68: MatCreate(PETSC_COMM_WORLD,&B);
69: MatLoad(B,viewer);
70: PetscViewerDestroy(&viewer);
71: } else {
72: PetscObjectReference((PetscObject)A_save);
73: B = A_save;
74: }
76: if (flgA) {
77: MatConvert(A_save,A_mattype,MAT_INPLACE_MATRIX,&A_save);
78: }
79: if (flgB) {
80: MatConvert(B,B_mattype,MAT_INPLACE_MATRIX,&B);
81: }
82: MatSetFromOptions(A_save);
83: MatSetFromOptions(B);
85: MatGetType(B,&mattype);
87: PetscMalloc(nzp*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
88: a = (PetscScalar*)(idxn + nzp);
90: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
91: PetscRandomSetFromOptions(rdm);
93: /* 1) MatMatMult() */
94: /* ----------------*/
95: if (Test_MatMatMult) {
96: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
98: /* (1.1) Test developer API */
99: MatProductCreate(A,B,NULL,&C);
100: MatSetOptionsPrefix(C,"AB_");
101: MatProductSetType(C,MATPRODUCT_AB);
102: MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);
103: MatProductSetFill(C,PETSC_DEFAULT);
104: MatProductSetFromOptions(C);
105: /* we can inquire about MATOP_PRODUCTSYMBOLIC even if the destination matrix type has not been set yet */
106: MatHasOperation(C,MATOP_PRODUCTSYMBOLIC,&flg);
107: MatProductSymbolic(C);
108: MatProductNumeric(C);
109: MatMatMultEqual(A,B,C,mcheck,&flg);
110: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B");
112: /* Test reuse symbolic C */
113: alpha = 0.9;
114: MatScale(A,alpha);
115: MatProductNumeric(C);
117: MatMatMultEqual(A,B,C,mcheck,&flg);
118: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B");
119: MatDestroy(&C);
121: /* (1.2) Test user driver */
122: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
124: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
125: alpha = 1.0;
126: for (i=0; i<2; i++) {
127: alpha -= 0.1;
128: MatScale(A,alpha);
129: MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
130: }
131: MatMatMultEqual(A,B,C,mcheck,&flg);
132: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()");
133: MatDestroy(&A);
135: /* Test MatProductClear() */
136: MatProductClear(C);
137: MatDestroy(&C);
139: /* Test MatMatMult() for dense and aij matrices */
140: PetscObjectTypeCompareAny((PetscObject)A,&flg,MATSEQAIJ,MATMPIAIJ,"");
141: if (flg) {
142: MatConvert(A_save,MATDENSE,MAT_INITIAL_MATRIX,&A);
143: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
144: MatDestroy(&C);
145: MatDestroy(&A);
146: }
147: }
149: /* Create P and R = P^T */
150: /* --------------------- */
151: MatGetSize(B,&PM,NULL);
152: if (PN < 0) PN = PM/2;
153: MatCreate(PETSC_COMM_WORLD,&P);
154: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,PM,PN);
155: MatSetType(P,mattype);
156: MatSeqAIJSetPreallocation(P,nzp,NULL);
157: MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
158: MatGetOwnershipRange(P,&rstart,&rend);
159: for (i=0; i<nzp; i++) {
160: PetscRandomGetValue(rdm,&a[i]);
161: }
162: for (i=rstart; i<rend; i++) {
163: for (j=0; j<nzp; j++) {
164: PetscRandomGetValue(rdm,&rval);
165: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
166: }
167: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
168: }
169: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
170: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
171: MatSetFromOptions(P);
173: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
174: MatSetFromOptions(R);
176: /* 2) MatTransposeMatMult() */
177: /* ------------------------ */
178: if (Test_MatTrMat) {
179: /* (2.1) Test developer driver C = P^T*B */
180: MatProductCreate(P,B,NULL,&C);
181: MatSetOptionsPrefix(C,"AtB_");
182: MatProductSetType(C,MATPRODUCT_AtB);
183: MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);
184: MatProductSetFill(C,PETSC_DEFAULT);
185: MatProductSetFromOptions(C);
186: MatProductSymbolic(C); /* equivalent to MatSetUp() */
187: MatSetOption(C,MAT_USE_INODES,PETSC_FALSE); /* illustrate how to call MatSetOption() */
188: MatProductNumeric(C);
189: MatProductNumeric(C); /* test reuse symbolic C */
191: MatTransposeMatMultEqual(P,B,C,mcheck,&flg);
192: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error: developer driver C = P^T*B");
193: MatDestroy(&C);
195: /* (2.2) Test user driver C = P^T*B */
196: MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
197: MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
198: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
199: MatProductClear(C);
201: /* Compare P^T*B and R*B */
202: MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1);
203: MatNormDifference(C,C1,&norm);
204: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatTransposeMatMult(): %g",(double)norm);
205: MatDestroy(&C1);
207: /* Test MatDuplicate() of C=P^T*B */
208: MatDuplicate(C,MAT_COPY_VALUES,&C1);
209: MatDestroy(&C1);
210: MatDestroy(&C);
211: }
213: /* 3) MatMatTransposeMult() */
214: /* ------------------------ */
215: if (Test_MatMatTr) {
216: /* C = B*R^T */
217: PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);
218: if (seqaij) {
219: MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
220: MatSetOptionsPrefix(C,"ABt_"); /* enable '-ABt_' for matrix C */
221: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
223: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
224: MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
226: /* Check */
227: MatMatMult(B,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1);
228: MatNormDifference(C,C1,&norm);
229: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatTransposeMult() %g",(double)norm);
230: MatDestroy(&C1);
231: MatDestroy(&C);
232: }
233: }
235: /* 4) Test MatPtAP() */
236: /*-------------------*/
237: if (Test_MatPtAP) {
238: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
240: /* (4.1) Test developer API */
241: MatProductCreate(A,P,NULL,&C);
242: MatSetOptionsPrefix(C,"PtAP_");
243: MatProductSetType(C,MATPRODUCT_PtAP);
244: MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);
245: MatProductSetFill(C,PETSC_DEFAULT);
246: MatProductSetFromOptions(C);
247: MatProductSymbolic(C);
248: MatProductNumeric(C);
249: MatPtAPMultEqual(A,P,C,mcheck,&flg);
250: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
251: MatProductNumeric(C); /* reuse symbolic C */
253: MatPtAPMultEqual(A,P,C,mcheck,&flg);
254: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
255: MatDestroy(&C);
257: /* (4.2) Test user driver */
258: MatPtAP(A,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
260: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
261: alpha = 1.0;
262: for (i=0; i<2; i++) {
263: alpha -= 0.1;
264: MatScale(A,alpha);
265: MatPtAP(A,P,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
266: }
267: MatPtAPMultEqual(A,P,C,mcheck,&flg);
268: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP");
270: /* 5) Test MatRARt() */
271: /* ----------------- */
272: if (Test_MatRARt) {
273: Mat RARt;
274: MatTranspose(P,MAT_REUSE_MATRIX,&R);
276: /* (5.1) Test developer driver RARt = R*A*Rt */
277: MatProductCreate(A,R,NULL,&RARt);
278: MatSetOptionsPrefix(RARt,"RARt_");
279: MatProductSetType(RARt,MATPRODUCT_RARt);
280: MatProductSetAlgorithm(RARt,MATPRODUCTALGORITHM_DEFAULT);
281: MatProductSetFill(RARt,PETSC_DEFAULT);
282: MatProductSetFromOptions(RARt);
283: MatProductSymbolic(RARt); /* equivalent to MatSetUp() */
284: MatSetOption(RARt,MAT_USE_INODES,PETSC_FALSE); /* illustrate how to call MatSetOption() */
285: MatProductNumeric(RARt);
286: MatProductNumeric(RARt); /* test reuse symbolic RARt */
287: MatDestroy(&RARt);
289: /* (2.2) Test user driver RARt = R*A*Rt */
290: MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);
291: MatRARt(A,R,MAT_REUSE_MATRIX,2.0,&RARt);
293: MatNormDifference(C,RARt,&norm);
294: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"|PtAP - RARt| = %g",(double)norm);
295: MatDestroy(&RARt);
296: }
298: MatDestroy(&A);
299: MatDestroy(&C);
300: }
302: /* Destroy objects */
303: PetscRandomDestroy(&rdm);
304: PetscFree(idxn);
306: MatDestroy(&A_save);
307: MatDestroy(&B);
308: MatDestroy(&P);
309: MatDestroy(&R);
311: PetscFinalize();
312: return ierr;
313: }
315: /*TEST
316: test:
317: suffix: 1
318: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
319: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
320: output_file: output/ex62_1.out
322: test:
323: suffix: 2_ab_scalable
324: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
325: 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
326: output_file: output/ex62_1.out
328: test:
329: suffix: 3_ab_scalable_fast
330: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
331: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via scalable_fast -matmatmult_via scalable_fast -matmattransmult_via color
332: output_file: output/ex62_1.out
334: test:
335: suffix: 4_ab_heap
336: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
337: 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
338: output_file: output/ex62_1.out
340: test:
341: suffix: 5_ab_btheap
342: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
343: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via btheap -matmatmult_via btheap -matrart_via r*art
344: output_file: output/ex62_1.out
346: test:
347: suffix: 6_ab_llcondensed
348: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
349: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via llcondensed -matmatmult_via llcondensed -matrart_via coloring_rart
350: output_file: output/ex62_1.out
352: test:
353: suffix: 7_ab_rowmerge
354: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
355: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via rowmerge -matmatmult_via rowmerge
356: output_file: output/ex62_1.out
358: test:
359: suffix: 8_ab_hypre
360: requires: hypre datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
361: 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
362: output_file: output/ex62_1.out
364: test:
365: suffix: 9_mkl
366: TODO: broken MatScale?
367: requires: mkl datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
368: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -A_mat_type aijmkl -B_mat_type aijmkl
369: output_file: output/ex62_1.out
371: test:
372: suffix: 10
373: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
374: nsize: 3
375: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
376: output_file: output/ex62_1.out
378: test:
379: suffix: 10_backend
380: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
381: nsize: 3
382: args: -fA ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via backend -matmatmult_via backend -AtB_matproduct_atb_via backend -mattransposematmult_via backend -PtAP_matproduct_ptap_via backend -matptap_via backend
383: output_file: output/ex62_1.out
385: test:
386: suffix: 11_ab_scalable
387: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
388: nsize: 3
389: 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
390: output_file: output/ex62_1.out
392: test:
393: suffix: 12_ab_seqmpi
394: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
395: nsize: 3
396: 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
397: output_file: output/ex62_1.out
399: test:
400: suffix: 13_ab_hypre
401: requires: hypre datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
402: nsize: 3
403: 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
404: output_file: output/ex62_1.out
406: test:
407: suffix: 14_seqaij
408: requires: !complex double !define(PETSC_USE_64BIT_INDICES)
409: args: -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
410: output_file: output/ex62_1.out
412: test:
413: suffix: 14_seqaijcusparse
414: requires: cuda !complex double !define(PETSC_USE_64BIT_INDICES)
415: args: -A_mat_type aijcusparse -B_mat_type aijcusparse -mat_form_explicit_transpose -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
416: output_file: output/ex62_1.out
418: test:
419: suffix: 14_mpiaijcusparse_seq
420: nsize: 1
421: requires: cuda !complex double !define(PETSC_USE_64BIT_INDICES)
422: args: -A_mat_type mpiaijcusparse -B_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
423: output_file: output/ex62_1.out
425: test:
426: suffix: 14_mpiaijcusparse
427: nsize: 3
428: requires: cuda !complex double !define(PETSC_USE_64BIT_INDICES)
429: args: -A_mat_type mpiaijcusparse -B_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
430: output_file: output/ex62_1.out
432: test:
433: suffix: 14_seqaijkokkos
434: requires: kokkos_kernels !complex double !define(PETSC_USE_64BIT_INDICES)
435: args: -A_mat_type aijkokkos -B_mat_type aijkokkos -fA ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system -fB ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
436: output_file: output/ex62_1.out
438: # these tests use matrices with many zero rows
439: test:
440: suffix: 15_seqaijcusparse
441: requires: cuda !complex double !define(PETSC_USE_64BIT_INDICES) datafilespath
442: args: -A_mat_type aijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
443: output_file: output/ex62_1.out
445: test:
446: suffix: 15_mpiaijcusparse_seq
447: nsize: 1
448: requires: cuda !complex double !define(PETSC_USE_64BIT_INDICES) datafilespath
449: args: -A_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
450: output_file: output/ex62_1.out
452: test:
453: nsize: 3
454: suffix: 15_mpiaijcusparse
455: requires: cuda !complex double !define(PETSC_USE_64BIT_INDICES) datafilespath
456: args: -A_mat_type mpiaijcusparse -mat_form_explicit_transpose -fA ${DATAFILESPATH}/matrices/matmatmult/A4.BGriffith
457: output_file: output/ex62_1.out
459: TEST*/