Actual source code: ex94.c
petsc-3.4.2 2013-07-02
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: PetscInt pm,pn,pM,pN;
29: MatInfo info;
31: PetscInitialize(&argc,&args,(char*)0,help);
32: MPI_Comm_size(PETSC_COMM_WORLD,&size);
33: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
35: PetscOptionsGetReal(NULL,"-fill",&fill,NULL);
37: /* Load the matrices A_save and B */
38: PetscOptionsGetString(NULL,"-f0",file[0],128,&flg);
39: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix A with the -f0 option.");
40: PetscOptionsGetString(NULL,"-f1",file[1],128,&flg);
41: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix B with the -f1 option.");
42: PetscOptionsGetString(NULL,"-f2",file[2],128,&flg);
43: if (!flg) {
44: preload = PETSC_FALSE;
45: } else {
46: PetscOptionsGetString(NULL,"-f3",file[3],128,&flg);
47: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for test matrix B with the -f3 option.");
48: }
50: PetscPreLoadBegin(preload,"Load system");
51: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PetscPreLoadIt],FILE_MODE_READ,&viewer);
52: MatCreate(PETSC_COMM_WORLD,&A_save);
53: MatLoad(A_save,viewer);
54: PetscViewerDestroy(&viewer);
56: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PetscPreLoadIt+1],FILE_MODE_READ,&viewer);
57: MatCreate(PETSC_COMM_WORLD,&B);
58: MatLoad(B,viewer);
59: PetscViewerDestroy(&viewer);
61: MatGetSize(B,&M,&N);
62: nzp = PetscMax((PetscInt)(0.1*M),5);
63: PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
64: a = (PetscScalar*)(idxn + nzp);
66: /* Create vectors v1 and v2 that are compatible with A_save */
67: VecCreate(PETSC_COMM_WORLD,&v1);
68: MatGetLocalSize(A_save,&m,NULL);
69: VecSetSizes(v1,m,PETSC_DECIDE);
70: VecSetFromOptions(v1);
71: VecDuplicate(v1,&v2);
73: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
74: PetscRandomSetFromOptions(rdm);
75: PetscOptionsGetReal(NULL,"-fill",&fill,NULL);
77: /* Test MatMatMult() */
78: /*-------------------*/
79: if (Test_MatMatMult) {
80: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
81: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
82: MatSetOptionsPrefix(C,"matmatmult_"); /* enable option '-matmatmult_' for matrix C */
83: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
84: /* PetscPrintf(PETSC_COMM_WORLD,"MatMatMult: nz_allocated = %g; nz_used = %g; nz_unneeded = %g\n",info.nz_allocated,info.nz_used, info.nz_unneeded); */
86: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
87: alpha=1.0;
88: for (i=0; i<2; i++) {
89: alpha -=0.1;
90: MatScale(A,alpha);
91: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
92: }
94: /* Create vector x that is compatible with B */
95: VecCreate(PETSC_COMM_WORLD,&x);
96: MatGetLocalSize(B,NULL,&n);
97: VecSetSizes(x,n,PETSC_DECIDE);
98: VecSetFromOptions(x);
100: norm = 0.0;
101: for (i=0; i<10; i++) {
102: VecSetRandom(x,rdm);
103: MatMult(B,x,v1);
104: MatMult(A,v1,v2); /* v2 = A*B*x */
105: MatMult(C,x,v1); /* v1 = C*x */
106: VecNorm(v1,NORM_2,&norm_abs);
107: VecAXPY(v1,none,v2);
108: VecNorm(v1,NORM_2,&norm_tmp);
110: norm_tmp /= norm_abs;
111: if (norm_tmp > norm) norm = norm_tmp;
112: }
113: if (norm >= tol) {
114: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|: %G\n",norm);
115: }
117: VecDestroy(&x);
118: MatDestroy(&A);
120: /* Test MatDuplicate() of C */
121: MatDuplicate(C,MAT_COPY_VALUES,&C1);
122: MatDestroy(&C1);
123: MatDestroy(&C);
124: } /* if (Test_MatMatMult) */
126: /* Test MatTransposeMatMult() and MatMatTransposeMult() */
127: /*------------------------------------------------------*/
128: if (Test_MatMatTr) {
129: /* Create P */
130: PetscInt PN,rstart,rend;
131: PN = M/2;
132: nzp = 5; /* num of nonzeros in each row of P */
133: MatCreate(PETSC_COMM_WORLD,&P);
134: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
135: MatSetType(P,MATAIJ);
136: MatSeqAIJSetPreallocation(P,nzp,NULL);
137: MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
138: MatGetOwnershipRange(P,&rstart,&rend);
139: for (i=0; i<nzp; i++) {
140: PetscRandomGetValue(rdm,&a[i]);
141: }
142: for (i=rstart; i<rend; i++) {
143: for (j=0; j<nzp; j++) {
144: PetscRandomGetValue(rdm,&rval);
145: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
146: }
147: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
148: }
149: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
150: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
152: /* Create R = P^T */
153: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
155: { /* Test R = P^T, C1 = R*B */
156: MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
157: MatTranspose(P,MAT_REUSE_MATRIX,&R);
158: MatMatMult(R,B,MAT_REUSE_MATRIX,fill,&C1);
159: MatDestroy(&C1);
160: }
162: /* C = P^T*B */
163: MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,fill,&C);
164: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
166: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
167: MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,fill,&C);
169: /* Compare P^T*B and R*B */
170: MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
171: MatEqual(C,C1,&flg);
172: if (!flg) {
173: /* Check norm of C1 = (-1.0)*C + C1 */
174: PetscReal nrm;
175: MatAXPY(C1,-1.0,C,DIFFERENT_NONZERO_PATTERN);
176: MatNorm(C1,NORM_INFINITY,&nrm);
177: if (nrm > 1.e-14) {
178: PetscPrintf(PETSC_COMM_WORLD,"Error in MatTransposeMatMult(): %g\n",nrm);
179: }
180: }
181: MatDestroy(&C1);
182: MatDestroy(&C);
184: /* C = B*R^T */
185: if (size == 1) {
186: MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,fill,&C);
187: MatSetOptionsPrefix(C,"matmatmulttr_"); /* enable '-matmatmulttr_' for matrix C */
188: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
190: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
191: MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,fill,&C);
193: /* Check */
194: MatMatMult(B,P,MAT_INITIAL_MATRIX,fill,&C1);
195: MatEqual(C,C1,&flg);
196: if (!flg) {
197: PetscPrintf(PETSC_COMM_WORLD,"Error in MatMatTransposeMult()\n");
198: }
199: MatDestroy(&C1);
200: MatDestroy(&C);
201: }
202: MatDestroy(&P);
203: MatDestroy(&R);
204: }
206: /* Test MatPtAP() */
207: /*----------------------*/
208: if (Test_MatPtAP) {
209: PetscInt PN;
210: Mat Cdup;
212: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
213: MatGetSize(A,&M,&N);
214: MatGetLocalSize(A,&m,&n);
215: /* PetscPrintf(PETSC_COMM_SELF,"[%d] A: %d,%d, %d,%d\n",rank,m,n,M,N); */
217: PN = M/2;
218: nzp = (PetscInt)(0.1*PN+1); /* num of nozeros in each row of P */
219: MatCreate(PETSC_COMM_WORLD,&P);
220: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
221: MatSetType(P,MATAIJ);
222: MatSeqAIJSetPreallocation(P,nzp,NULL);
223: MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
224: for (i=0; i<nzp; i++) {
225: PetscRandomGetValue(rdm,&a[i]);
226: }
227: MatGetOwnershipRange(P,&rstart,&rend);
228: for (i=rstart; i<rend; i++) {
229: for (j=0; j<nzp; j++) {
230: PetscRandomGetValue(rdm,&rval);
231: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
232: }
233: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
234: }
235: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
236: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
238: /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
239: MatGetSize(P,&pM,&pN);
240: MatGetLocalSize(P,&pm,&pn);
241: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
242: /* if (!rank) {PetscPrintf(PETSC_COMM_SELF," MatPtAP() is done, P, %d, %d, %d,%d\n",pm,pn,pM,pN);} */
244: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
245: alpha=1.0;
246: for (i=0; i<2; i++) {
247: alpha -=0.1;
248: MatScale(A,alpha);
249: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
250: }
252: /* Test MatDuplicate() */
253: MatDuplicate(C,MAT_COPY_VALUES,&Cdup);
254: MatDestroy(&Cdup);
256: if (size>1) Test_MatRARt = PETSC_FALSE;
257: /* Test MatRARt() */
258: if (Test_MatRARt) {
259: Mat R, RARt;
260: PetscBool equal;
261: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
262: MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);
263: MatEqual(C,RARt,&equal);
264: if (!equal) {
265: PetscReal norm;
266: MatAXPY(RARt,-1.0,C,DIFFERENT_NONZERO_PATTERN); /* RARt = -RARt + C */
267: MatNorm(RARt,NORM_FROBENIUS,&norm);
268: if (norm > 1.e-12) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"|PtAP - RARt| = %g",norm);
269: }
270: MatDestroy(&R);
271: MatDestroy(&RARt);
272: }
274: if (Test_MatMatMatMult && size == 1) {
275: Mat R, RAP;
276: PetscBool equal;
277: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
278: MatMatMatMult(R,A,P,MAT_INITIAL_MATRIX,2.0,&RAP);
279: MatEqual(C,RAP,&equal);
280: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PtAP != RAP");
281: MatDestroy(&R);
282: MatDestroy(&RAP);
283: }
285: /* Create vector x that is compatible with P */
286: VecCreate(PETSC_COMM_WORLD,&x);
287: MatGetLocalSize(P,&m,&n);
288: VecSetSizes(x,n,PETSC_DECIDE);
289: VecSetFromOptions(x);
291: VecCreate(PETSC_COMM_WORLD,&v3);
292: VecSetSizes(v3,n,PETSC_DECIDE);
293: VecSetFromOptions(v3);
294: VecDuplicate(v3,&v4);
296: norm = 0.0;
297: for (i=0; i<10; i++) {
298: VecSetRandom(x,rdm);
299: MatMult(P,x,v1);
300: MatMult(A,v1,v2); /* v2 = A*P*x */
302: MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
303: MatMult(C,x,v4); /* v3 = C*x */
304: VecNorm(v4,NORM_2,&norm_abs);
305: VecAXPY(v4,none,v3);
306: VecNorm(v4,NORM_2,&norm_tmp);
308: norm_tmp /= norm_abs;
309: if (norm_tmp > norm) norm = norm_tmp;
310: }
311: if (norm >= tol) {
312: PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v1 - v2|: %G\n",norm);
313: }
315: MatDestroy(&A);
316: MatDestroy(&P);
317: MatDestroy(&C);
318: VecDestroy(&v3);
319: VecDestroy(&v4);
320: VecDestroy(&x);
321: }
323: /* Destroy objects */
324: VecDestroy(&v1);
325: VecDestroy(&v2);
326: PetscRandomDestroy(&rdm);
327: PetscFree(idxn);
329: MatDestroy(&A_save);
330: MatDestroy(&B);
332: PetscPreLoadEnd();
333: PetscFinalize();
334: return 0;
335: }