Actual source code: ex62.c
petsc-3.13.6 2020-09-29
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,"default");
84: MatProductSetFill(C,PETSC_DEFAULT);
85: MatProductSetFromOptions(C);
86: MatProductSymbolic(C);
87: MatProductNumeric(C);
89: /* Test reuse symbolic C */
90: alpha = 0.9;
91: MatScale(A,alpha);
92: MatProductNumeric(C);
94: MatMatMultEqual(A,B,C,10,&flg);
95: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B");
96: MatDestroy(&C);
98: /* (1.2) Test user driver */
99: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
101: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
102: alpha = 1.0;
103: for (i=0; i<2; i++) {
104: alpha -= 0.1;
105: MatScale(A,alpha);
106: MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
107: }
108: MatMatMultEqual(A,B,C,10,&flg);
109: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()");
110: MatDestroy(&A);
112: /* Test MatProductClear() */
113: MatProductClear(C);
114: MatDestroy(&C);
116: /* Test MatMatMult() for dense and aij matrices */
117: MatConvert(A_save,MATDENSE,MAT_INITIAL_MATRIX,&A);
118: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
119: MatDestroy(&C);
120: MatDestroy(&A);
122: }
124: /* Create P and R = P^T */
125: /* --------------------- */
126: PN = M/2;
127: nzp = 5; /* num of nonzeros in each row of P */
128: MatCreate(PETSC_COMM_WORLD,&P);
129: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
130: MatSetType(P,mattype);
131: MatSeqAIJSetPreallocation(P,nzp,NULL);
132: MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
133: MatGetOwnershipRange(P,&rstart,&rend);
134: for (i=0; i<nzp; i++) {
135: PetscRandomGetValue(rdm,&a[i]);
136: }
137: for (i=rstart; i<rend; i++) {
138: for (j=0; j<nzp; j++) {
139: PetscRandomGetValue(rdm,&rval);
140: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
141: }
142: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
143: }
144: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
145: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
147: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
149: /* 2) MatTransposeMatMult() */
150: /* ------------------------ */
151: if (Test_MatTrMat) {
152: /* (2.1) Test developer driver C = P^T*B */
153: MatProductCreate(P,B,NULL,&C);
154: MatSetOptionsPrefix(C,"AtB_");
155: MatProductSetType(C,MATPRODUCT_AtB);
156: MatProductSetAlgorithm(C,"default");
157: MatProductSetFill(C,PETSC_DEFAULT);
158: MatProductSetFromOptions(C);
159: MatProductSymbolic(C); /* equivalent to MatSetUp() */
160: MatSetOption(C,MAT_USE_INODES,PETSC_FALSE); /* illustrate how to call MatSetOption() */
161: MatProductNumeric(C);
162: MatProductNumeric(C); /* test reuse symbolic C */
164: MatTransposeMatMultEqual(P,B,C,10,&flg);
165: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error: developer driver C = P^T*B");
166: MatDestroy(&C);
168: /* (2.2) Test user driver C = P^T*B */
169: MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
170: MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
171: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
172: MatFreeIntermediateDataStructures(C);
174: /* Compare P^T*B and R*B */
175: MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1);
176: MatNormDifference(C,C1,&norm);
177: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatTransposeMatMult(): %g",(double)norm);
178: MatDestroy(&C1);
180: /* Test MatDuplicate() of C=P^T*B */
181: MatDuplicate(C,MAT_COPY_VALUES,&C1);
182: MatDestroy(&C1);
183: MatDestroy(&C);
184: }
186: /* 3) MatMatTransposeMult() */
187: /* ------------------------ */
188: if (Test_MatMatTr) {
189: /* C = B*R^T */
190: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);
191: if (size == 1 && seqaij) {
192: MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
193: MatSetOptionsPrefix(C,"ABt_"); /* enable '-ABt_' for matrix C */
194: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
196: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
197: MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
199: /* Check */
200: MatMatMult(B,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1);
201: MatNormDifference(C,C1,&norm);
202: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatTransposeMult() %g",(double)norm);
203: MatDestroy(&C1);
204: MatDestroy(&C);
205: }
206: }
208: /* 4) Test MatPtAP() */
209: /*-------------------*/
210: if (Test_MatPtAP) {
211: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
213: /* (4.1) Test developer API */
214: MatProductCreate(A,P,NULL,&C);
215: MatSetOptionsPrefix(C,"PtAP_");
216: MatProductSetType(C,MATPRODUCT_PtAP);
217: MatProductSetAlgorithm(C,"default");
218: MatProductSetFill(C,PETSC_DEFAULT);
219: MatProductSetFromOptions(C);
220: MatProductSymbolic(C);
221: MatProductNumeric(C);
222: MatProductNumeric(C); /* reuse symbolic C */
224: MatPtAPMultEqual(A,P,C,10,&flg);
225: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
226: MatDestroy(&C);
228: /* (4.2) Test user driver */
229: MatPtAP(A,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
231: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
232: alpha=1.0;
233: for (i=0; i<2; i++) {
234: alpha -= 0.1;
235: MatScale(A,alpha);
236: MatPtAP(A,P,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
237: }
239: MatPtAPMultEqual(A,P,C,10,&flg);
240: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP");
242: /* 5) Test MatRARt() */
243: /* ----------------- */
244: if (Test_MatRARt) {
245: Mat RARt;
246: MatTranspose(P,MAT_REUSE_MATRIX,&R);
248: /* (5.1) Test developer driver RARt = R*A*Rt */
249: MatProductCreate(A,R,NULL,&RARt);
250: MatSetOptionsPrefix(RARt,"RARt_");
251: MatProductSetType(RARt,MATPRODUCT_RARt);
252: MatProductSetAlgorithm(RARt,"default");
253: MatProductSetFill(RARt,PETSC_DEFAULT);
254: MatProductSetFromOptions(RARt);
255: MatProductSymbolic(RARt); /* equivalent to MatSetUp() */
256: MatSetOption(RARt,MAT_USE_INODES,PETSC_FALSE); /* illustrate how to call MatSetOption() */
257: MatProductNumeric(RARt);
258: MatProductNumeric(RARt); /* test reuse symbolic RARt */
259: MatDestroy(&RARt);
261: /* (2.2) Test user driver RARt = R*A*Rt */
262: MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);
263: MatRARt(A,R,MAT_REUSE_MATRIX,2.0,&RARt);
265: MatNormDifference(C,RARt,&norm);
266: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"|PtAP - RARt| = %g",(double)norm);
267: MatDestroy(&RARt);
268: }
270: MatDestroy(&A);
271: MatDestroy(&C);
272: }
274: /* Destroy objects */
275: PetscRandomDestroy(&rdm);
276: PetscFree(idxn);
278: MatDestroy(&A_save);
279: MatDestroy(&B);
280: MatDestroy(&P);
281: MatDestroy(&R);
283: PetscFinalize();
284: return ierr;
285: }
287: /*TEST
288: test:
289: suffix: 1
290: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
291: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
292: output_file: output/ex62_1.out
294: test:
295: suffix: 2_ab_scalable
296: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
297: 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
298: output_file: output/ex62_1.out
300: test:
301: suffix: 3_ab_scalable_fast
302: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
303: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via scalable_fast -matmatmult_via scalable_fast -matmattransmult_via color
304: output_file: output/ex62_1.out
306: test:
307: suffix: 4_ab_heap
308: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
309: 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
310: output_file: output/ex62_1.out
312: test:
313: suffix: 5_ab_btheap
314: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
315: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via btheap -matmatmult_via btheap -matrart_via r*art
316: output_file: output/ex62_1.out
318: test:
319: suffix: 6_ab_llcondensed
320: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
321: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via llcondensed -matmatmult_via llcondensed -matrart_via coloring_rart
322: output_file: output/ex62_1.out
324: test:
325: suffix: 7_ab_rowmerge
326: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
327: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via rowmerge -matmatmult_via rowmerge
328: output_file: output/ex62_1.out
330: test:
331: suffix: 8_ab_hypre
332: requires: hypre datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
333: 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
334: output_file: output/ex62_1.out
336: test:
337: suffix: 10
338: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
339: nsize: 3
340: args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium
341: output_file: output/ex62_1.out
343: test:
344: suffix: 11_ab_scalable
345: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
346: nsize: 3
347: 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
348: output_file: output/ex62_1.out
350: test:
351: suffix: 12_ab_seqmpi
352: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
353: nsize: 3
354: 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
355: output_file: output/ex62_1.out
357: test:
358: suffix: 13_ab_hypre
359: requires: hypre datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
360: nsize: 3
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*/