Actual source code: ex94.c
petsc-3.9.4 2018-09-11
2: static char help[] = "Tests sequential and parallel MatMatMult() and MatPtAP(), MatTransposeMatMult(), sequential MatMatTransposeMult(), MatRARt()\n\
3: Input arguments are:\n\
4: -f0 <input_file> -f1 <input_file> -f2 <input_file> -f3 <input_file> : file to load\n\n";
5: /* Example of usage:
6: ./ex94 -f0 <A_binary> -f1 <B_binary> -matmatmult_mat_view ascii::ascii_info -matmatmulttr_mat_view
7: mpiexec -n 3 ./ex94 -f0 medium -f1 medium -f2 arco1 -f3 arco1 -matmatmult_mat_view
8: */
9:
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,P,R,C,C1;
29: Vec x,v1,v2,v3,v4;
30: PetscViewer viewer;
32: PetscMPIInt size,rank;
33: PetscInt i,m,n,j,*idxn,M,N,nzp,rstart,rend;
34: PetscReal norm,norm_abs,norm_tmp,fill=4.0;
35: PetscRandom rdm;
36: char file[4][128];
37: PetscBool flg,preload = PETSC_TRUE;
38: PetscScalar *a,rval,alpha,none = -1.0;
39: PetscBool Test_MatMatMult=PETSC_TRUE,Test_MatMatTr=PETSC_TRUE,Test_MatPtAP=PETSC_TRUE,Test_MatRARt=PETSC_TRUE,Test_MatMatMatMult=PETSC_TRUE;
40: PetscBool Test_MatAXPY=PETSC_FALSE;
41: PetscInt pm,pn,pM,pN;
42: MatInfo info;
43: PetscBool seqaij;
44: MatType mattype;
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: PetscOptionsGetReal(NULL,NULL,"-fill",&fill,NULL);
52: /* Load the matrices A_save and B */
53: PetscOptionsGetString(NULL,NULL,"-f0",file[0],sizeof(file[0]),&flg);
54: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix A with the -f0 option.");
55: PetscOptionsGetString(NULL,NULL,"-f1",file[1],sizeof(file[1]),&flg);
56: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix B with the -f1 option.");
57: PetscOptionsGetString(NULL,NULL,"-f2",file[2],sizeof(file[2]),&flg);
58: if (!flg) {
59: preload = PETSC_FALSE;
60: } else {
61: PetscOptionsGetString(NULL,NULL,"-f3",file[3],128,&flg);
62: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for test matrix B with the -f3 option.");
63: }
65: PetscPreLoadBegin(preload,"Load system");
66: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PetscPreLoadIt],FILE_MODE_READ,&viewer);
67: MatCreate(PETSC_COMM_WORLD,&A_save);
68: MatSetFromOptions(A_save);
69: MatLoad(A_save,viewer);
70: PetscViewerDestroy(&viewer);
72: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PetscPreLoadIt+1],FILE_MODE_READ,&viewer);
73: MatCreate(PETSC_COMM_WORLD,&B);
74: MatSetFromOptions(B);
75: MatLoad(B,viewer);
76: PetscViewerDestroy(&viewer);
78: MatGetType(B,&mattype);
80: MatGetSize(B,&M,&N);
81: nzp = PetscMax((PetscInt)(0.1*M),5);
82: PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
83: a = (PetscScalar*)(idxn + nzp);
85: /* Create vectors v1 and v2 that are compatible with A_save */
86: VecCreate(PETSC_COMM_WORLD,&v1);
87: MatGetLocalSize(A_save,&m,NULL);
88: VecSetSizes(v1,m,PETSC_DECIDE);
89: VecSetFromOptions(v1);
90: VecDuplicate(v1,&v2);
92: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
93: PetscRandomSetFromOptions(rdm);
94: PetscOptionsGetReal(NULL,NULL,"-fill",&fill,NULL);
96: /* Test MatAXPY() */
97: /*-------------------*/
98: PetscOptionsHasName(NULL,NULL,"-test_MatAXPY",&Test_MatAXPY);
99: if (Test_MatAXPY) {
100: Mat Btmp;
101: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
102: MatDuplicate(B,MAT_COPY_VALUES,&Btmp);
103: MatAXPY(A,-1.0,B,DIFFERENT_NONZERO_PATTERN); /* A = -B + A_save */
105: MatScale(A,-1.0); /* A = -A = B - A_save */
106: MatAXPY(Btmp,-1.0,A,DIFFERENT_NONZERO_PATTERN); /* Btmp = -A + B = A_save */
107: MatMultEqual(A_save,Btmp,10,&flg);
108: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatAXPY() is incorrect\n");
109: MatDestroy(&A);
110: MatDestroy(&Btmp);
112: Test_MatMatMult = PETSC_FALSE;
113: Test_MatMatTr = PETSC_FALSE;
114: Test_MatPtAP = PETSC_FALSE;
115: Test_MatRARt = PETSC_FALSE;
116: Test_MatMatMatMult = PETSC_FALSE;
117: }
119: /* 1) Test MatMatMult() */
120: /* ---------------------*/
121: if (Test_MatMatMult) {
122: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
123: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
124: MatSetOptionsPrefix(C,"matmatmult_"); /* enable option '-matmatmult_' for matrix C */
125: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
127: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
128: alpha=1.0;
129: for (i=0; i<2; i++) {
130: alpha -=0.1;
131: MatScale(A,alpha);
132: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
133: }
134: MatMatMultEqual(A,B,C,10,&flg);
135: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()\n");
136: MatDestroy(&A);
138: /* Test MatDuplicate() of C=A*B */
139: MatDuplicate(C,MAT_COPY_VALUES,&C1);
140: MatDestroy(&C1);
141: MatDestroy(&C);
142: } /* if (Test_MatMatMult) */
144: /* 2) Test MatTransposeMatMult() and MatMatTransposeMult() */
145: /* ------------------------------------------------------- */
146: if (Test_MatMatTr) {
147: /* Create P */
148: PetscInt PN,rstart,rend;
149: PN = M/2;
150: nzp = 5; /* num of nonzeros in each row of P */
151: MatCreate(PETSC_COMM_WORLD,&P);
152: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
153: MatSetType(P,mattype);
154: MatSeqAIJSetPreallocation(P,nzp,NULL);
155: MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
156: MatGetOwnershipRange(P,&rstart,&rend);
157: for (i=0; i<nzp; i++) {
158: PetscRandomGetValue(rdm,&a[i]);
159: }
160: for (i=rstart; i<rend; i++) {
161: for (j=0; j<nzp; j++) {
162: PetscRandomGetValue(rdm,&rval);
163: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
164: }
165: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
166: }
167: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
168: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
170: /* Create R = P^T */
171: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
173: { /* Test R = P^T, C1 = R*B */
174: MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
175: MatTranspose(P,MAT_REUSE_MATRIX,&R);
176: MatMatMult(R,B,MAT_REUSE_MATRIX,fill,&C1);
177: MatDestroy(&C1);
178: }
180: /* C = P^T*B */
181: MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,fill,&C);
182: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
184: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
185: MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,fill,&C);
187: /* Compare P^T*B and R*B */
188: MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
189: MatNormDifference(C,C1,&norm);
190: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatTransposeMatMult(): %g\n",(double)norm);
191: MatDestroy(&C1);
193: /* Test MatDuplicate() of C=P^T*B */
194: MatDuplicate(C,MAT_COPY_VALUES,&C1);
195: MatDestroy(&C1);
196: MatDestroy(&C);
198: /* C = B*R^T */
199: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);
200: if (size == 1 && seqaij) {
201: MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,fill,&C);
202: MatSetOptionsPrefix(C,"matmatmulttr_"); /* enable '-matmatmulttr_' for matrix C */
203: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
205: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
206: MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,fill,&C);
208: /* Check */
209: MatMatMult(B,P,MAT_INITIAL_MATRIX,fill,&C1);
210: MatNormDifference(C,C1,&norm);
211: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatTransposeMult() %g\n",(double)norm);
212: MatDestroy(&C1);
213: MatDestroy(&C);
214: }
215: MatDestroy(&P);
216: MatDestroy(&R);
217: }
219: /* 3) Test MatPtAP() */
220: /*-------------------*/
221: if (Test_MatPtAP) {
222: PetscInt PN;
223: Mat Cdup;
225: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
226: MatGetSize(A,&M,&N);
227: MatGetLocalSize(A,&m,&n);
229: PN = M/2;
230: nzp = (PetscInt)(0.1*PN+1); /* num of nozeros in each row of P */
231: MatCreate(PETSC_COMM_WORLD,&P);
232: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
233: MatSetType(P,mattype);
234: MatSeqAIJSetPreallocation(P,nzp,NULL);
235: MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
236: for (i=0; i<nzp; i++) {
237: PetscRandomGetValue(rdm,&a[i]);
238: }
239: MatGetOwnershipRange(P,&rstart,&rend);
240: for (i=rstart; i<rend; i++) {
241: for (j=0; j<nzp; j++) {
242: PetscRandomGetValue(rdm,&rval);
243: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
244: }
245: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
246: }
247: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
248: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
250: /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
251: MatGetSize(P,&pM,&pN);
252: MatGetLocalSize(P,&pm,&pn);
253: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
255: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
256: alpha=1.0;
257: for (i=0; i<2; i++) {
258: alpha -=0.1;
259: MatScale(A,alpha);
260: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
261: }
263: /* Test PtAP ops with P SeqDense and A either SeqAIJ or SeqDense (it assumes MatPtAP_SeqAIJ_SeqAIJ is fine) */
264: if (size == 1) {
265: Mat Cdensetest,Pdense,Cdense,Adense;
266: PetscReal norm;
268: MatConvert(C,MATSEQDENSE,MAT_INITIAL_MATRIX,&Cdensetest);
269: MatConvert(P,MATSEQDENSE,MAT_INITIAL_MATRIX,&Pdense);
271: /* test with A SeqAIJ */
272: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&seqaij);
273: if (seqaij) {
274: MatPtAP(A,Pdense,MAT_INITIAL_MATRIX,fill,&Cdense);
275: MatAXPY(Cdense,-1.0,Cdensetest,SAME_NONZERO_PATTERN);
276: MatNorm(Cdense,NORM_FROBENIUS,&norm);
277: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatPtAP with A SeqAIJ and P SeqDense: %g\n",(double)norm);
278: MatScale(Cdense,-1.);
279: MatPtAP(A,Pdense,MAT_REUSE_MATRIX,fill,&Cdense);
280: MatAXPY(Cdense,-1.0,Cdensetest,SAME_NONZERO_PATTERN);
281: MatNorm(Cdense,NORM_FROBENIUS,&norm);
282: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatPtAP with A SeqAIJ and P SeqDense and MAT_REUSE_MATRIX: %g\n",(double)norm);
283: MatDestroy(&Cdense);
284: }
286: /* test with A SeqDense */
287: MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&Adense);
288: MatPtAP(Adense,Pdense,MAT_INITIAL_MATRIX,fill,&Cdense);
289: MatAXPY(Cdense,-1.0,Cdensetest,SAME_NONZERO_PATTERN);
290: MatNorm(Cdense,NORM_FROBENIUS,&norm);
291: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatPtAP with A SeqDense and P SeqDense: %g\n",(double)norm);
292: MatScale(Cdense,-1.);
293: MatPtAP(Adense,Pdense,MAT_REUSE_MATRIX,fill,&Cdense);
294: MatAXPY(Cdense,-1.0,Cdensetest,SAME_NONZERO_PATTERN);
295: MatNorm(Cdense,NORM_FROBENIUS,&norm);
296: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatPtAP with A SeqDense and P SeqDense and MAT_REUSE_MATRIX: %g\n",(double)norm);
297: MatDestroy(&Cdense);
298: MatDestroy(&Cdensetest);
299: MatDestroy(&Pdense);
300: MatDestroy(&Adense);
301: }
303: /* Test MatDuplicate() of C=PtAP */
304: MatDuplicate(C,MAT_COPY_VALUES,&Cdup);
305: MatDestroy(&Cdup);
307: if (size>1 || !seqaij) Test_MatRARt = PETSC_FALSE;
308: /* 4) Test MatRARt() */
309: /* ----------------- */
310: if (Test_MatRARt) {
311: Mat R, RARt;
312: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
313: MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);
314: MatNormDifference(C,RARt,&norm);
315: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"|PtAP - RARt| = %g",(double)norm);
316: MatDestroy(&R);
317: MatDestroy(&RARt);
318: }
320: if (Test_MatMatMatMult && size == 1) {
321: Mat R, RAP;
322: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
323: MatMatMatMult(R,A,P,MAT_INITIAL_MATRIX,2.0,&RAP);
324: MatNormDifference(C,RAP,&norm);
325: if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PtAP != RAP %g",(double)norm);
326: MatDestroy(&R);
327: MatDestroy(&RAP);
328: }
330: /* Create vector x that is compatible with P */
331: VecCreate(PETSC_COMM_WORLD,&x);
332: MatGetLocalSize(P,&m,&n);
333: VecSetSizes(x,n,PETSC_DECIDE);
334: VecSetFromOptions(x);
336: VecCreate(PETSC_COMM_WORLD,&v3);
337: VecSetSizes(v3,n,PETSC_DECIDE);
338: VecSetFromOptions(v3);
339: VecDuplicate(v3,&v4);
341: norm = 0.0;
342: for (i=0; i<10; i++) {
343: VecSetRandom(x,rdm);
344: MatMult(P,x,v1);
345: MatMult(A,v1,v2); /* v2 = A*P*x */
347: MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
348: MatMult(C,x,v4); /* v3 = C*x */
349: VecNorm(v4,NORM_2,&norm_abs);
350: VecAXPY(v4,none,v3);
351: VecNorm(v4,NORM_2,&norm_tmp);
353: norm_tmp /= norm_abs;
354: if (norm_tmp > norm) norm = norm_tmp;
355: }
356: if (norm >= PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatPtAP(), |v1 - v2|: %g\n",(double)norm);
358: MatDestroy(&A);
359: MatDestroy(&P);
360: MatDestroy(&C);
361: VecDestroy(&v3);
362: VecDestroy(&v4);
363: VecDestroy(&x);
364: }
366: /* Destroy objects */
367: VecDestroy(&v1);
368: VecDestroy(&v2);
369: PetscRandomDestroy(&rdm);
370: PetscFree(idxn);
372: MatDestroy(&A_save);
373: MatDestroy(&B);
375: PetscPreLoadEnd();
376: PetscFinalize();
377: return ierr;
378: }
382: /*TEST
384: test:
385: suffix: 2_mattransposematmult_matmatmult
386: nsize: 3
387: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
388: args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -mattransposematmult_via matmatmult> ex94_2.tmp 2>&1
390: test:
391: suffix: 2_mattransposematmult_scalable
392: nsize: 3
393: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
394: args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -mattransposematmult_via scalable> ex94_2.tmp 2>&1
395: output_file: output/ex94_1.out
397: test:
398: suffix: axpy_mpiaij
399: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
400: nsize: 8
401: args: -f0 ${DATAFILESPATH}/matrices/EigenProblems/jifengzhao/petsc_stiff20.dat -f1 ${DATAFILESPATH}/matrices/EigenProblems/jifengzhao/petsc_mass20.dat -test_MatAXPY
402: output_file: output/ex94_1.out
404: test:
405: suffix: axpy_mpibaij
406: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
407: nsize: 8
408: args: -f0 ${DATAFILESPATH}/matrices/EigenProblems/jifengzhao/petsc_stiff20.dat -f1 ${DATAFILESPATH}/matrices/EigenProblems/jifengzhao/petsc_mass20.dat -test_MatAXPY -mat_type baij
409: output_file: output/ex94_1.out
411: test:
412: suffix: axpy_mpisbaij
413: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
414: nsize: 8
415: args: -f0 ${DATAFILESPATH}/matrices/EigenProblems/jifengzhao/petsc_stiff20.dat -f1 ${DATAFILESPATH}/matrices/EigenProblems/jifengzhao/petsc_mass20.dat -test_MatAXPY -mat_type sbaij
416: output_file: output/ex94_1.out
418: test:
419: suffix: matmatmult
420: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
421: args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info
422: output_file: output/ex94_1.out
424: test:
425: suffix: matmatmult_2
426: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
427: args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -mat_type mpiaij -viewer_binary_skip_info
428: output_file: output/ex94_1.out
430: test:
431: suffix: matmatmult_scalable
432: nsize: 4
433: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
434: args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -matmatmult_via scalable
435: output_file: output/ex94_1.out
437: test:
438: suffix: ptap
439: nsize: 3
440: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
441: args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -matptap_via scalable
442: output_file: output/ex94_1.out
444: test:
445: suffix: rap
446: nsize: 3
447: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
448: args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium
449: output_file: output/ex94_1.out
451: test:
452: suffix: scalable0
453: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
454: args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info
455: output_file: output/ex94_1.out
457: test:
458: suffix: scalable1
459: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
460: args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info -matptap_via scalable
461: output_file: output/ex94_1.out
463: TEST*/