Actual source code: ex62.c
petsc-3.14.6 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,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;
42: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
43: MPI_Comm_size(PETSC_COMM_WORLD,&size);
44: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
46: /* Load the matrices A_save and B */
47: PetscOptionsGetString(NULL,NULL,"-fA",file[0],sizeof(file[0]),&flg);
48: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix A with the -fA option.");
49: PetscOptionsGetString(NULL,NULL,"-fB",file[1],sizeof(file[1]),&flg);
50: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix B with the -fB option.");
52: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);
53: MatCreate(PETSC_COMM_WORLD,&A_save);
54: MatSetFromOptions(A_save);
55: MatLoad(A_save,viewer);
56: PetscViewerDestroy(&viewer);
58: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);
59: MatCreate(PETSC_COMM_WORLD,&B);
60: MatSetFromOptions(B);
61: MatLoad(B,viewer);
62: PetscViewerDestroy(&viewer);
64: MatGetType(B,&mattype);
66: MatGetSize(B,&M,&N);
67: nzp = PetscMax((PetscInt)(0.1*M),5);
68: PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
69: a = (PetscScalar*)(idxn + nzp);
71: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
72: PetscRandomSetFromOptions(rdm);
74: /* 1) MatMatMult() */
75: /* ----------------*/
76: if (Test_MatMatMult) {
77: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
79: /* (1.1) Test developer API */
80: MatProductCreate(A,B,NULL,&C);
81: MatSetOptionsPrefix(C,"AB_");
82: MatProductSetType(C,MATPRODUCT_AB);
83: MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);
84: MatProductSetFill(C,PETSC_DEFAULT);
85: MatProductSetFromOptions(C);
86: /* we can inquire about MATOP_PRODUCTSYMBOLIC even if the destination matrix type has not been set yet */
87: MatHasOperation(C,MATOP_PRODUCTSYMBOLIC,&flg);
88: MatProductSymbolic(C);
89: MatProductNumeric(C);
91: /* Test reuse symbolic C */
92: alpha = 0.9;
93: MatScale(A,alpha);
94: MatProductNumeric(C);
96: MatMatMultEqual(A,B,C,10,&flg);
97: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B");
98: MatDestroy(&C);
100: /* (1.2) Test user driver */
101: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
103: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
104: alpha = 1.0;
105: for (i=0; i<2; i++) {
106: alpha -= 0.1;
107: MatScale(A,alpha);
108: MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
109: }
110: MatMatMultEqual(A,B,C,10,&flg);
111: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()");
112: MatDestroy(&A);
114: /* Test MatProductClear() */
115: MatProductClear(C);
116: MatDestroy(&C);
118: /* Test MatMatMult() for dense and aij matrices */
119: MatConvert(A_save,MATDENSE,MAT_INITIAL_MATRIX,&A);
120: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
121: MatDestroy(&C);
122: MatDestroy(&A);
124: }
126: /* Create P and R = P^T */
127: /* --------------------- */
128: PN = M/2;
129: nzp = 5; /* num of nonzeros in each row of P */
130: MatCreate(PETSC_COMM_WORLD,&P);
131: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
132: MatSetType(P,mattype);
133: MatSeqAIJSetPreallocation(P,nzp,NULL);
134: MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
135: MatGetOwnershipRange(P,&rstart,&rend);
136: for (i=0; i<nzp; i++) {
137: PetscRandomGetValue(rdm,&a[i]);
138: }
139: for (i=rstart; i<rend; i++) {
140: for (j=0; j<nzp; j++) {
141: PetscRandomGetValue(rdm,&rval);
142: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
143: }
144: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
145: }
146: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
147: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
149: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
151: /* 2) MatTransposeMatMult() */
152: /* ------------------------ */
153: if (Test_MatTrMat) {
154: /* (2.1) Test developer driver C = P^T*B */
155: MatProductCreate(P,B,NULL,&C);
156: MatSetOptionsPrefix(C,"AtB_");
157: MatProductSetType(C,MATPRODUCT_AtB);
158: MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);
159: MatProductSetFill(C,PETSC_DEFAULT);
160: MatProductSetFromOptions(C);
161: MatProductSymbolic(C); /* equivalent to MatSetUp() */
162: MatSetOption(C,MAT_USE_INODES,PETSC_FALSE); /* illustrate how to call MatSetOption() */
163: MatProductNumeric(C);
164: MatProductNumeric(C); /* test reuse symbolic C */
166: MatTransposeMatMultEqual(P,B,C,10,&flg);
167: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error: developer driver C = P^T*B");
168: MatDestroy(&C);
170: /* (2.2) Test user driver C = P^T*B */
171: MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
172: MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
173: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
174: MatProductClear(C);
176: /* Compare P^T*B and R*B */
177: MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1);
178: MatNormDifference(C,C1,&norm);
179: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatTransposeMatMult(): %g",(double)norm);
180: MatDestroy(&C1);
182: /* Test MatDuplicate() of C=P^T*B */
183: MatDuplicate(C,MAT_COPY_VALUES,&C1);
184: MatDestroy(&C1);
185: MatDestroy(&C);
186: }
188: /* 3) MatMatTransposeMult() */
189: /* ------------------------ */
190: if (Test_MatMatTr) {
191: /* C = B*R^T */
192: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);
193: if (size == 1 && seqaij) {
194: MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
195: MatSetOptionsPrefix(C,"ABt_"); /* enable '-ABt_' for matrix C */
196: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
198: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
199: MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
201: /* Check */
202: MatMatMult(B,P,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 MatMatTransposeMult() %g",(double)norm);
205: MatDestroy(&C1);
206: MatDestroy(&C);
207: }
208: }
210: /* 4) Test MatPtAP() */
211: /*-------------------*/
212: if (Test_MatPtAP) {
213: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
215: /* (4.1) Test developer API */
216: MatProductCreate(A,P,NULL,&C);
217: MatSetOptionsPrefix(C,"PtAP_");
218: MatProductSetType(C,MATPRODUCT_PtAP);
219: MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);
220: MatProductSetFill(C,PETSC_DEFAULT);
221: MatProductSetFromOptions(C);
222: MatProductSymbolic(C);
223: MatProductNumeric(C);
224: MatProductNumeric(C); /* reuse symbolic C */
226: MatPtAPMultEqual(A,P,C,10,&flg);
227: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
228: MatDestroy(&C);
230: /* (4.2) Test user driver */
231: MatPtAP(A,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
233: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
234: alpha=1.0;
235: for (i=0; i<2; i++) {
236: alpha -= 0.1;
237: MatScale(A,alpha);
238: MatPtAP(A,P,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
239: }
241: MatPtAPMultEqual(A,P,C,10,&flg);
242: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP");
244: /* 5) Test MatRARt() */
245: /* ----------------- */
246: if (Test_MatRARt) {
247: Mat RARt;
248: MatTranspose(P,MAT_REUSE_MATRIX,&R);
250: /* (5.1) Test developer driver RARt = R*A*Rt */
251: MatProductCreate(A,R,NULL,&RARt);
252: MatSetOptionsPrefix(RARt,"RARt_");
253: MatProductSetType(RARt,MATPRODUCT_RARt);
254: MatProductSetAlgorithm(RARt,MATPRODUCTALGORITHM_DEFAULT);
255: MatProductSetFill(RARt,PETSC_DEFAULT);
256: MatProductSetFromOptions(RARt);
257: MatProductSymbolic(RARt); /* equivalent to MatSetUp() */
258: MatSetOption(RARt,MAT_USE_INODES,PETSC_FALSE); /* illustrate how to call MatSetOption() */
259: MatProductNumeric(RARt);
260: MatProductNumeric(RARt); /* test reuse symbolic RARt */
261: MatDestroy(&RARt);
263: /* (2.2) Test user driver RARt = R*A*Rt */
264: MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);
265: MatRARt(A,R,MAT_REUSE_MATRIX,2.0,&RARt);
267: MatNormDifference(C,RARt,&norm);
268: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"|PtAP - RARt| = %g",(double)norm);
269: MatDestroy(&RARt);
270: }
272: MatDestroy(&A);
273: MatDestroy(&C);
274: }
276: /* Destroy objects */
277: PetscRandomDestroy(&rdm);
278: PetscFree(idxn);
280: MatDestroy(&A_save);
281: MatDestroy(&B);
282: MatDestroy(&P);
283: MatDestroy(&R);
285: PetscFinalize();
286: return ierr;
287: }
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
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
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
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
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
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
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
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
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
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
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
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
366: TEST*/