Actual source code: ex94.c
petsc-3.8.4 2018-03-24
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: */
10: #include <petscmat.h>
12: int main(int argc,char **args)
13: {
14: Mat A,A_save,B,P,R,C,C1;
15: Vec x,v1,v2,v3,v4;
16: PetscViewer viewer;
18: PetscMPIInt size,rank;
19: PetscInt i,m,n,j,*idxn,M,N,nzp,rstart,rend;
20: PetscReal norm,norm_abs,norm_tmp,tol=1.e-8,fill=4.0;
21: PetscRandom rdm;
22: char file[4][128];
23: PetscBool flg,preload = PETSC_TRUE;
24: PetscScalar *a,rval,alpha,none = -1.0;
25: PetscBool Test_MatMatMult=PETSC_TRUE,Test_MatMatTr=PETSC_TRUE,Test_MatPtAP=PETSC_TRUE,Test_MatRARt=PETSC_TRUE,Test_MatMatMatMult=PETSC_TRUE;
26: PetscBool Test_MatAXPY=PETSC_FALSE;
27: PetscInt pm,pn,pM,pN;
28: MatInfo info;
29: PetscBool seqaij;
30: MatType mattype;
32: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
33: MPI_Comm_size(PETSC_COMM_WORLD,&size);
34: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
36: PetscOptionsGetReal(NULL,NULL,"-fill",&fill,NULL);
38: /* Load the matrices A_save and B */
39: PetscOptionsGetString(NULL,NULL,"-f0",file[0],sizeof(file[0]),&flg);
40: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix A with the -f0 option.");
41: PetscOptionsGetString(NULL,NULL,"-f1",file[1],sizeof(file[1]),&flg);
42: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix B with the -f1 option.");
43: PetscOptionsGetString(NULL,NULL,"-f2",file[2],sizeof(file[2]),&flg);
44: if (!flg) {
45: preload = PETSC_FALSE;
46: } else {
47: PetscOptionsGetString(NULL,NULL,"-f3",file[3],128,&flg);
48: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for test matrix B with the -f3 option.");
49: }
51: PetscPreLoadBegin(preload,"Load system");
52: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PetscPreLoadIt],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[2*PetscPreLoadIt+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: /* Create vectors v1 and v2 that are compatible with A_save */
72: VecCreate(PETSC_COMM_WORLD,&v1);
73: MatGetLocalSize(A_save,&m,NULL);
74: VecSetSizes(v1,m,PETSC_DECIDE);
75: VecSetFromOptions(v1);
76: VecDuplicate(v1,&v2);
78: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
79: PetscRandomSetFromOptions(rdm);
80: PetscOptionsGetReal(NULL,NULL,"-fill",&fill,NULL);
82: /* Test MatAXPY() */
83: /*-------------------*/
84: PetscOptionsHasName(NULL,NULL,"-test_MatAXPY",&Test_MatAXPY);
85: if (Test_MatAXPY) {
86: Mat Btmp;
87: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
88: MatDuplicate(B,MAT_COPY_VALUES,&Btmp);
89: MatAXPY(A,-1.0,B,DIFFERENT_NONZERO_PATTERN); /* A = -B + A_save */
91: MatScale(A,-1.0); /* A = -A = B - A_save */
92: MatAXPY(Btmp,-1.0,A,DIFFERENT_NONZERO_PATTERN); /* Btmp = -A + B = A_save */
93: MatMultEqual(A_save,Btmp,10,&flg);
94: if (!flg) SETERRQ(PETSC_COMM_SELF,0,"MatAXPY() is incorrect\n");
95: MatDestroy(&A);
96: MatDestroy(&Btmp);
98: Test_MatMatMult = PETSC_FALSE;
99: Test_MatMatTr = PETSC_FALSE;
100: Test_MatPtAP = PETSC_FALSE;
101: Test_MatRARt = PETSC_FALSE;
102: Test_MatMatMatMult = PETSC_FALSE;
103: }
105: /* 1) Test MatMatMult() */
106: /* ---------------------*/
107: if (Test_MatMatMult) {
108: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
109: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
110: MatSetOptionsPrefix(C,"matmatmult_"); /* enable option '-matmatmult_' for matrix C */
111: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
113: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
114: alpha=1.0;
115: for (i=0; i<2; i++) {
116: alpha -=0.1;
117: MatScale(A,alpha);
118: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
119: }
120: MatMatMultEqual(A,B,C,10,&flg);
121: if (!flg) {
122: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult()\n");
123: }
124: MatDestroy(&A);
126: /* Test MatDuplicate() of C=A*B */
127: MatDuplicate(C,MAT_COPY_VALUES,&C1);
128: MatDestroy(&C1);
129: MatDestroy(&C);
130: } /* if (Test_MatMatMult) */
132: /* 2) Test MatTransposeMatMult() and MatMatTransposeMult() */
133: /* ------------------------------------------------------- */
134: if (Test_MatMatTr) {
135: /* Create P */
136: PetscInt PN,rstart,rend;
137: PN = M/2;
138: nzp = 5; /* num of nonzeros in each row of P */
139: MatCreate(PETSC_COMM_WORLD,&P);
140: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
141: MatSetType(P,mattype);
142: MatSeqAIJSetPreallocation(P,nzp,NULL);
143: MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
144: MatGetOwnershipRange(P,&rstart,&rend);
145: for (i=0; i<nzp; i++) {
146: PetscRandomGetValue(rdm,&a[i]);
147: }
148: for (i=rstart; i<rend; i++) {
149: for (j=0; j<nzp; j++) {
150: PetscRandomGetValue(rdm,&rval);
151: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
152: }
153: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
154: }
155: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
156: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
158: /* Create R = P^T */
159: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
161: { /* Test R = P^T, C1 = R*B */
162: MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
163: MatTranspose(P,MAT_REUSE_MATRIX,&R);
164: MatMatMult(R,B,MAT_REUSE_MATRIX,fill,&C1);
165: MatDestroy(&C1);
166: }
168: /* C = P^T*B */
169: MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,fill,&C);
170: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
172: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
173: MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,fill,&C);
175: /* Compare P^T*B and R*B */
176: MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
177: MatEqual(C,C1,&flg);
178: if (!flg) {
179: /* Check norm of C1 = (-1.0)*C + C1 */
180: PetscReal nrm;
181: MatAXPY(C1,-1.0,C,DIFFERENT_NONZERO_PATTERN);
182: MatNorm(C1,NORM_INFINITY,&nrm);
183: if (nrm > 1.e-14) {
184: PetscPrintf(PETSC_COMM_WORLD,"Error in MatTransposeMatMult(): %g\n",nrm);
185: }
186: }
187: MatDestroy(&C1);
189: /* Test MatDuplicate() of C=P^T*B */
190: MatDuplicate(C,MAT_COPY_VALUES,&C1);
191: MatDestroy(&C1);
192: MatDestroy(&C);
194: /* C = B*R^T */
195: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);
196: if (size == 1 && seqaij) {
197: MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,fill,&C);
198: MatSetOptionsPrefix(C,"matmatmulttr_"); /* enable '-matmatmulttr_' for matrix C */
199: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
201: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
202: MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,fill,&C);
204: /* Check */
205: MatMatMult(B,P,MAT_INITIAL_MATRIX,fill,&C1);
206: MatEqual(C,C1,&flg);
207: if (!flg) {
208: PetscPrintf(PETSC_COMM_WORLD,"Error in MatMatTransposeMult()\n");
209: }
210: MatDestroy(&C1);
211: MatDestroy(&C);
212: }
213: MatDestroy(&P);
214: MatDestroy(&R);
215: }
217: /* 3) Test MatPtAP() */
218: /*-------------------*/
219: if (Test_MatPtAP) {
220: PetscInt PN;
221: Mat Cdup;
223: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
224: MatGetSize(A,&M,&N);
225: MatGetLocalSize(A,&m,&n);
227: PN = M/2;
228: nzp = (PetscInt)(0.1*PN+1); /* num of nozeros in each row of P */
229: MatCreate(PETSC_COMM_WORLD,&P);
230: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
231: MatSetType(P,mattype);
232: MatSeqAIJSetPreallocation(P,nzp,NULL);
233: MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
234: for (i=0; i<nzp; i++) {
235: PetscRandomGetValue(rdm,&a[i]);
236: }
237: MatGetOwnershipRange(P,&rstart,&rend);
238: for (i=rstart; i<rend; i++) {
239: for (j=0; j<nzp; j++) {
240: PetscRandomGetValue(rdm,&rval);
241: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
242: }
243: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
244: }
245: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
246: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
248: /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
249: MatGetSize(P,&pM,&pN);
250: MatGetLocalSize(P,&pm,&pn);
251: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
253: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
254: alpha=1.0;
255: for (i=0; i<2; i++) {
256: alpha -=0.1;
257: MatScale(A,alpha);
258: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
259: }
261: /* Test PtAP ops with P SeqDense and A either SeqAIJ or SeqDense (it assumes MatPtAP_SeqAIJ_SeqAIJ is fine) */
262: if (size == 1) {
263: Mat Cdensetest,Pdense,Cdense,Adense;
264: PetscReal norm;
266: MatConvert(C,MATSEQDENSE,MAT_INITIAL_MATRIX,&Cdensetest);
267: MatConvert(P,MATSEQDENSE,MAT_INITIAL_MATRIX,&Pdense);
269: /* test with A SeqAIJ */
270: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&seqaij);
271: if (seqaij) {
272: MatPtAP(A,Pdense,MAT_INITIAL_MATRIX,fill,&Cdense);
273: MatAXPY(Cdense,-1.0,Cdensetest,SAME_NONZERO_PATTERN);
274: MatNorm(Cdense,NORM_FROBENIUS,&norm);
275: if (norm > PETSC_SMALL) {
276: PetscPrintf(PETSC_COMM_SELF,"Error in MatPtAP with A SeqAIJ and P SeqDense: %g\n",norm);
277: }
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) {
283: PetscPrintf(PETSC_COMM_SELF,"Error in MatPtAP with A SeqAIJ and P SeqDense and MAT_REUSE_MATRIX: %g\n",norm);
284: }
285: MatDestroy(&Cdense);
286: }
288: /* test with A SeqDense */
289: MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&Adense);
290: MatPtAP(Adense,Pdense,MAT_INITIAL_MATRIX,fill,&Cdense);
291: MatAXPY(Cdense,-1.0,Cdensetest,SAME_NONZERO_PATTERN);
292: MatNorm(Cdense,NORM_FROBENIUS,&norm);
293: if (norm > PETSC_SMALL) {
294: PetscPrintf(PETSC_COMM_SELF,"Error in MatPtAP with A SeqDense and P SeqDense: %g\n",norm);
295: }
296: MatScale(Cdense,-1.);
297: MatPtAP(Adense,Pdense,MAT_REUSE_MATRIX,fill,&Cdense);
298: MatAXPY(Cdense,-1.0,Cdensetest,SAME_NONZERO_PATTERN);
299: MatNorm(Cdense,NORM_FROBENIUS,&norm);
300: if (norm > PETSC_SMALL) {
301: PetscPrintf(PETSC_COMM_SELF,"Error in MatPtAP with A SeqDense and P SeqDense and MAT_REUSE_MATRIX: %g\n",norm);
302: }
303: MatDestroy(&Cdense);
304: MatDestroy(&Cdensetest);
305: MatDestroy(&Pdense);
306: MatDestroy(&Adense);
307: }
309: /* Test MatDuplicate() of C=PtAP */
310: MatDuplicate(C,MAT_COPY_VALUES,&Cdup);
311: MatDestroy(&Cdup);
313: if (size>1 || !seqaij) Test_MatRARt = PETSC_FALSE;
314: /* 4) Test MatRARt() */
315: /* ----------------- */
316: if (Test_MatRARt) {
317: Mat R, RARt;
318: PetscBool equal;
319: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
320: MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);
321: MatEqual(C,RARt,&equal);
322: if (!equal) {
323: PetscReal norm;
324: MatAXPY(RARt,-1.0,C,DIFFERENT_NONZERO_PATTERN); /* RARt = -RARt + C */
325: MatNorm(RARt,NORM_FROBENIUS,&norm);
326: if (norm > 1.e-12) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"|PtAP - RARt| = %g",norm);
327: }
328: MatDestroy(&R);
329: MatDestroy(&RARt);
330: }
332: if (Test_MatMatMatMult && size == 1) {
333: Mat R, RAP;
334: PetscBool equal;
335: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
336: MatMatMatMult(R,A,P,MAT_INITIAL_MATRIX,2.0,&RAP);
337: MatEqual(C,RAP,&equal);
338: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PtAP != RAP");
339: MatDestroy(&R);
340: MatDestroy(&RAP);
341: }
343: /* Create vector x that is compatible with P */
344: VecCreate(PETSC_COMM_WORLD,&x);
345: MatGetLocalSize(P,&m,&n);
346: VecSetSizes(x,n,PETSC_DECIDE);
347: VecSetFromOptions(x);
349: VecCreate(PETSC_COMM_WORLD,&v3);
350: VecSetSizes(v3,n,PETSC_DECIDE);
351: VecSetFromOptions(v3);
352: VecDuplicate(v3,&v4);
354: norm = 0.0;
355: for (i=0; i<10; i++) {
356: VecSetRandom(x,rdm);
357: MatMult(P,x,v1);
358: MatMult(A,v1,v2); /* v2 = A*P*x */
360: MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
361: MatMult(C,x,v4); /* v3 = C*x */
362: VecNorm(v4,NORM_2,&norm_abs);
363: VecAXPY(v4,none,v3);
364: VecNorm(v4,NORM_2,&norm_tmp);
366: norm_tmp /= norm_abs;
367: if (norm_tmp > norm) norm = norm_tmp;
368: }
369: if (norm >= tol) {
370: PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v1 - v2|: %g\n",(double)norm);
371: }
373: MatDestroy(&A);
374: MatDestroy(&P);
375: MatDestroy(&C);
376: VecDestroy(&v3);
377: VecDestroy(&v4);
378: VecDestroy(&x);
379: }
381: /* Destroy objects */
382: VecDestroy(&v1);
383: VecDestroy(&v2);
384: PetscRandomDestroy(&rdm);
385: PetscFree(idxn);
387: MatDestroy(&A_save);
388: MatDestroy(&B);
390: PetscPreLoadEnd();
391: PetscFinalize();
392: return ierr;
393: }