Actual source code: ex94.c
petsc-3.6.4 2016-04-12
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_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>
14: int main(int argc,char **args)
15: {
16: Mat A,A_save,B,P,R,C,C1;
17: Vec x,v1,v2,v3,v4;
18: PetscViewer viewer;
20: PetscMPIInt size,rank;
21: PetscInt i,m,n,j,*idxn,M,N,nzp,rstart,rend;
22: PetscReal norm,norm_abs,norm_tmp,tol=1.e-8,fill=4.0;
23: PetscRandom rdm;
24: char file[4][128];
25: PetscBool flg,preload = PETSC_TRUE;
26: PetscScalar *a,rval,alpha,none = -1.0;
27: PetscBool Test_MatMatMult=PETSC_TRUE,Test_MatMatTr=PETSC_TRUE,Test_MatPtAP=PETSC_TRUE,Test_MatRARt=PETSC_TRUE,Test_MatMatMatMult=PETSC_TRUE;
28: PetscBool Test_MatAXPY=PETSC_FALSE;
29: PetscInt pm,pn,pM,pN;
30: MatInfo info;
32: PetscInitialize(&argc,&args,(char*)0,help);
33: MPI_Comm_size(PETSC_COMM_WORLD,&size);
34: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
36: PetscOptionsGetReal(NULL,"-fill",&fill,NULL);
38: /* Load the matrices A_save and B */
39: PetscOptionsGetString(NULL,"-f0",file[0],128,&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,"-f1",file[1],128,&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,"-f2",file[2],128,&flg);
44: if (!flg) {
45: preload = PETSC_FALSE;
46: } else {
47: PetscOptionsGetString(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: MatGetSize(B,&M,&N);
65: nzp = PetscMax((PetscInt)(0.1*M),5);
66: PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
67: a = (PetscScalar*)(idxn + nzp);
69: /* Create vectors v1 and v2 that are compatible with A_save */
70: VecCreate(PETSC_COMM_WORLD,&v1);
71: MatGetLocalSize(A_save,&m,NULL);
72: VecSetSizes(v1,m,PETSC_DECIDE);
73: VecSetFromOptions(v1);
74: VecDuplicate(v1,&v2);
76: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
77: PetscRandomSetFromOptions(rdm);
78: PetscOptionsGetReal(NULL,"-fill",&fill,NULL);
80: /* Test MatAXPY() */
81: /*-------------------*/
82: PetscOptionsHasName(NULL,"-test_MatAXPY",&Test_MatAXPY);
83: if (Test_MatAXPY) {
84: Mat Btmp;
85: /* if (!rank) printf(" Loading matrices is done...\n"); */
86: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
87: MatDuplicate(B,MAT_COPY_VALUES,&Btmp);
88: MatAXPY(A,-1.0,B,DIFFERENT_NONZERO_PATTERN); /* A = -B + A_save */
90: /* if (!rank) printf(" Test_MatAXPY is done, now checking accuracy ...\n"); */
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: /* 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);
112: /* PetscPrintf(PETSC_COMM_WORLD,"MatMatMult: nz_allocated = %g; nz_used = %g; nz_unneeded = %g\n",info.nz_allocated,info.nz_used, info.nz_unneeded); */
114: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
115: alpha=1.0;
116: for (i=0; i<2; i++) {
117: alpha -=0.1;
118: MatScale(A,alpha);
119: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
120: }
122: /* Create vector x that is compatible with B */
123: VecCreate(PETSC_COMM_WORLD,&x);
124: MatGetLocalSize(B,NULL,&n);
125: VecSetSizes(x,n,PETSC_DECIDE);
126: VecSetFromOptions(x);
128: norm = 0.0;
129: for (i=0; i<10; i++) {
130: VecSetRandom(x,rdm);
131: MatMult(B,x,v1);
132: MatMult(A,v1,v2); /* v2 = A*B*x */
133: MatMult(C,x,v1); /* v1 = C*x */
134: VecNorm(v1,NORM_2,&norm_abs);
135: VecAXPY(v1,none,v2);
136: VecNorm(v1,NORM_2,&norm_tmp);
138: norm_tmp /= norm_abs;
139: if (norm_tmp > norm) norm = norm_tmp;
140: }
141: if (norm >= tol) {
142: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|: %g\n",(double)norm);
143: }
145: VecDestroy(&x);
146: MatDestroy(&A);
148: /* Test MatDuplicate() of C */
149: MatDuplicate(C,MAT_COPY_VALUES,&C1);
150: MatDestroy(&C1);
151: MatDestroy(&C);
152: } /* if (Test_MatMatMult) */
154: /* Test MatTransposeMatMult() and MatMatTransposeMult() */
155: /*------------------------------------------------------*/
156: if (Test_MatMatTr) {
157: /* Create P */
158: PetscInt PN,rstart,rend;
159: PN = M/2;
160: nzp = 5; /* num of nonzeros in each row of P */
161: MatCreate(PETSC_COMM_WORLD,&P);
162: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
163: MatSetType(P,MATAIJ);
164: MatSeqAIJSetPreallocation(P,nzp,NULL);
165: MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
166: MatGetOwnershipRange(P,&rstart,&rend);
167: for (i=0; i<nzp; i++) {
168: PetscRandomGetValue(rdm,&a[i]);
169: }
170: for (i=rstart; i<rend; i++) {
171: for (j=0; j<nzp; j++) {
172: PetscRandomGetValue(rdm,&rval);
173: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
174: }
175: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
176: }
177: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
178: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
180: /* Create R = P^T */
181: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
183: { /* Test R = P^T, C1 = R*B */
184: MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
185: MatTranspose(P,MAT_REUSE_MATRIX,&R);
186: MatMatMult(R,B,MAT_REUSE_MATRIX,fill,&C1);
187: MatDestroy(&C1);
188: }
190: /* C = P^T*B */
191: MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,fill,&C);
192: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
194: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
195: MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,fill,&C);
197: /* Compare P^T*B and R*B */
198: MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
199: MatEqual(C,C1,&flg);
200: if (!flg) {
201: /* Check norm of C1 = (-1.0)*C + C1 */
202: PetscReal nrm;
203: MatAXPY(C1,-1.0,C,DIFFERENT_NONZERO_PATTERN);
204: MatNorm(C1,NORM_INFINITY,&nrm);
205: if (nrm > 1.e-14) {
206: PetscPrintf(PETSC_COMM_WORLD,"Error in MatTransposeMatMult(): %g\n",nrm);
207: }
208: }
209: MatDestroy(&C1);
210: MatDestroy(&C);
212: /* C = B*R^T */
213: if (size == 1) {
214: MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,fill,&C);
215: MatSetOptionsPrefix(C,"matmatmulttr_"); /* enable '-matmatmulttr_' for matrix C */
216: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
218: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
219: MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,fill,&C);
221: /* Check */
222: MatMatMult(B,P,MAT_INITIAL_MATRIX,fill,&C1);
223: MatEqual(C,C1,&flg);
224: if (!flg) {
225: PetscPrintf(PETSC_COMM_WORLD,"Error in MatMatTransposeMult()\n");
226: }
227: MatDestroy(&C1);
228: MatDestroy(&C);
229: }
230: MatDestroy(&P);
231: MatDestroy(&R);
232: }
234: /* Test MatPtAP() */
235: /*----------------------*/
236: if (Test_MatPtAP) {
237: PetscInt PN;
238: Mat Cdup;
240: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
241: MatGetSize(A,&M,&N);
242: MatGetLocalSize(A,&m,&n);
243: /* PetscPrintf(PETSC_COMM_SELF,"[%d] A: %d,%d, %d,%d\n",rank,m,n,M,N); */
245: PN = M/2;
246: nzp = (PetscInt)(0.1*PN+1); /* num of nozeros in each row of P */
247: MatCreate(PETSC_COMM_WORLD,&P);
248: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
249: MatSetType(P,MATAIJ);
250: MatSeqAIJSetPreallocation(P,nzp,NULL);
251: MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
252: for (i=0; i<nzp; i++) {
253: PetscRandomGetValue(rdm,&a[i]);
254: }
255: MatGetOwnershipRange(P,&rstart,&rend);
256: for (i=rstart; i<rend; i++) {
257: for (j=0; j<nzp; j++) {
258: PetscRandomGetValue(rdm,&rval);
259: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
260: }
261: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
262: }
263: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
264: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
266: /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
267: MatGetSize(P,&pM,&pN);
268: MatGetLocalSize(P,&pm,&pn);
269: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
270: /* if (!rank) {PetscPrintf(PETSC_COMM_SELF," MatPtAP() is done, P, %d, %d, %d,%d\n",pm,pn,pM,pN);} */
272: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
273: alpha=1.0;
274: for (i=0; i<2; i++) {
275: alpha -=0.1;
276: MatScale(A,alpha);
277: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
278: }
280: /* Test MatDuplicate() */
281: MatDuplicate(C,MAT_COPY_VALUES,&Cdup);
282: MatDestroy(&Cdup);
284: if (size>1) Test_MatRARt = PETSC_FALSE;
285: /* Test MatRARt() */
286: if (Test_MatRARt) {
287: Mat R, RARt;
288: PetscBool equal;
289: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
290: MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);
291: MatEqual(C,RARt,&equal);
292: if (!equal) {
293: PetscReal norm;
294: MatAXPY(RARt,-1.0,C,DIFFERENT_NONZERO_PATTERN); /* RARt = -RARt + C */
295: MatNorm(RARt,NORM_FROBENIUS,&norm);
296: if (norm > 1.e-12) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"|PtAP - RARt| = %g",norm);
297: }
298: MatDestroy(&R);
299: MatDestroy(&RARt);
300: }
302: if (Test_MatMatMatMult && size == 1) {
303: Mat R, RAP;
304: PetscBool equal;
305: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
306: MatMatMatMult(R,A,P,MAT_INITIAL_MATRIX,2.0,&RAP);
307: MatEqual(C,RAP,&equal);
308: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PtAP != RAP");
309: MatDestroy(&R);
310: MatDestroy(&RAP);
311: }
313: /* Create vector x that is compatible with P */
314: VecCreate(PETSC_COMM_WORLD,&x);
315: MatGetLocalSize(P,&m,&n);
316: VecSetSizes(x,n,PETSC_DECIDE);
317: VecSetFromOptions(x);
319: VecCreate(PETSC_COMM_WORLD,&v3);
320: VecSetSizes(v3,n,PETSC_DECIDE);
321: VecSetFromOptions(v3);
322: VecDuplicate(v3,&v4);
324: norm = 0.0;
325: for (i=0; i<10; i++) {
326: VecSetRandom(x,rdm);
327: MatMult(P,x,v1);
328: MatMult(A,v1,v2); /* v2 = A*P*x */
330: MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
331: MatMult(C,x,v4); /* v3 = C*x */
332: VecNorm(v4,NORM_2,&norm_abs);
333: VecAXPY(v4,none,v3);
334: VecNorm(v4,NORM_2,&norm_tmp);
336: norm_tmp /= norm_abs;
337: if (norm_tmp > norm) norm = norm_tmp;
338: }
339: if (norm >= tol) {
340: PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v1 - v2|: %g\n",(double)norm);
341: }
343: MatDestroy(&A);
344: MatDestroy(&P);
345: MatDestroy(&C);
346: VecDestroy(&v3);
347: VecDestroy(&v4);
348: VecDestroy(&x);
349: }
351: /* Destroy objects */
352: VecDestroy(&v1);
353: VecDestroy(&v2);
354: PetscRandomDestroy(&rdm);
355: PetscFree(idxn);
357: MatDestroy(&A_save);
358: MatDestroy(&B);
360: PetscPreLoadEnd();
361: PetscFinalize();
362: return 0;
363: }