Actual source code: ex94.c
2: static char help[] = "Tests sequential and parallel MatMatMult() and MatPtAP(), sequential MatMatMultTranspose()\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_info -matmatmulttr_mat_view_info
7: mpiexec -n 3 ./ex94 -f0 medium -f1 medium -f2 arco1 -f3 arco1 -matmatmult_mat_view_info
8: */
10: #include <petscmat.h>
14: int main(int argc,char **args)
15: {
16: Mat A,A_save,B,P,C,C1;
17: Vec x,v1,v2;
18: PetscViewer viewer;
20: PetscMPIInt size,rank;
21: PetscInt i,m,n,j,*idxn,M,N,nzp;
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_MatMatMultTr=PETSC_TRUE;
28: Vec v3,v4,v5;
29: PetscInt pm,pn,pM,pN;
30: PetscBool Test_MatPtAP=PETSC_TRUE;
31: MatInfo info;
32:
33: PetscInitialize(&argc,&args,(char *)0,help);
34: MPI_Comm_size(PETSC_COMM_WORLD,&size);
35: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
37: PetscOptionsGetReal(PETSC_NULL,"-fill",&fill,PETSC_NULL);
39: /* Load the matrices A and B */
40: PetscOptionsGetString(PETSC_NULL,"-f0",file[0],128,&flg);
41: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix A with the -f0 option.");
42: PetscOptionsGetString(PETSC_NULL,"-f1",file[1],128,&flg);
43: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix B with the -f1 option.");
44: PetscOptionsGetString(PETSC_NULL,"-f2",file[2],128,&flg);
45: if (!flg) {
46: preload = PETSC_FALSE;
47: } else {
48: PetscOptionsGetString(PETSC_NULL,"-f3",file[3],128,&flg);
49: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for test matrix B with the -f3 option.");
50: }
52: PetscPreLoadBegin(preload,"Load system");
53: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PetscPreLoadIt],FILE_MODE_READ,&viewer);
54: MatCreate(PETSC_COMM_WORLD,&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: MatLoad(B,viewer);
61: PetscViewerDestroy(&viewer);
63: MatGetSize(B,&M,&N);
64: nzp = (PetscInt)(0.1*M);
65: PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
66: a = (PetscScalar*)(idxn + nzp);
67:
68: /* Create vectors v1 and v2 that are compatible with A_save */
69: VecCreate(PETSC_COMM_WORLD,&v1);
70: MatGetLocalSize(A_save,&m,PETSC_NULL);
71: VecSetSizes(v1,m,PETSC_DECIDE);
72: VecSetFromOptions(v1);
73: VecDuplicate(v1,&v2);
75: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
76: PetscRandomSetFromOptions(rdm);
77: PetscOptionsGetReal(PETSC_NULL,"-fill",&fill,PETSC_NULL);
79: /* Test MatMatMult() */
80: /*-------------------*/
81: if (Test_MatMatMult){
82: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
83: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
84: MatSetOptionsPrefix(C,"matmatmult_"); /* enable option '-matmatmult_' for matrix C */
85: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
86: PetscPrintf(PETSC_COMM_WORLD,"MatMatMult: nz_allocated = %g; nz_used = %g; nz_unneeded = %g\n",info.nz_allocated,info.nz_used, info.nz_unneeded);
88: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
89: alpha=1.0;
90: for (i=0; i<2; i++){
91: alpha -=0.1;
92: MatScale(A,alpha);
93: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
94: }
96: /* Create vector x that is compatible with B */
97: VecCreate(PETSC_COMM_WORLD,&x);
98: MatGetLocalSize(B,PETSC_NULL,&n);
99: VecSetSizes(x,n,PETSC_DECIDE);
100: VecSetFromOptions(x);
102: norm = 0.0;
103: for (i=0; i<10; i++) {
104: VecSetRandom(x,rdm);
105: MatMult(B,x,v1);
106: MatMult(A,v1,v2); /* v2 = A*B*x */
107: MatMult(C,x,v1); /* v1 = C*x */
108: VecNorm(v1,NORM_2,&norm_abs);
109: VecAXPY(v1,none,v2);
110: VecNorm(v1,NORM_2,&norm_tmp);
111: norm_tmp /= norm_abs;
112: if (norm_tmp > norm) norm = norm_tmp;
113: }
114: if (norm >= tol) {
115: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|: %G\n",norm);
116: }
118: VecDestroy(&x);
119: MatDestroy(&A);
121: /* Test MatDuplicate() of C */
122: MatDuplicate(C,MAT_COPY_VALUES,&C1);
123: MatDestroy(&C1);
124: MatDestroy(&C);
125: } /* if (Test_MatMatMult) */
127: /* Test MatMatMultTranspose() */
128: /*----------------------------*/
129: if (size>1) Test_MatMatMultTr = PETSC_FALSE;
130: if (Test_MatMatMultTr){
131: PetscInt PN;
132: PN = M/2;
133: nzp = 5; /* num of nonzeros in each row of P */
134: MatCreate(PETSC_COMM_WORLD,&P);
135: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
136: MatSetType(P,MATAIJ);
137: MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
138: MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
139: for (i=0; i<nzp; i++){
140: PetscRandomGetValue(rdm,&a[i]);
141: }
142: for (i=0; i<M; 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);
151:
152: MatMatMultTranspose(P,B,MAT_INITIAL_MATRIX,fill,&C);
153: MatSetOptionsPrefix(C,"matmatmulttr_"); /* enable '-matmatmulttr_' for matrix C */
154: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
155: PetscPrintf(PETSC_COMM_WORLD,"MatMatMultTranspose: nz_allocated = %g; nz_used = %g; nz_unneeded = %g\n",info.nz_allocated,info.nz_used, info.nz_unneeded);
157: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
158: alpha=1.0;
159: for (i=0; i<2; i++){
160: alpha -=0.1;
161: MatScale(P,alpha);
162: MatMatMultTranspose(P,B,MAT_REUSE_MATRIX,fill,&C);
163: }
165: /* Create vector x, v5 that are compatible with B */
166: VecCreate(PETSC_COMM_WORLD,&x);
167: MatGetLocalSize(B,&m,&n);
168: VecSetSizes(x,n,PETSC_DECIDE);
169: VecSetFromOptions(x);
171: VecCreate(PETSC_COMM_WORLD,&v5);
172: VecSetSizes(v5,m,PETSC_DECIDE);
173: VecSetFromOptions(v5);
174:
175: MatGetLocalSize(P,PETSC_NULL,&n);
176: VecCreate(PETSC_COMM_WORLD,&v3);
177: VecSetSizes(v3,n,PETSC_DECIDE);
178: VecSetFromOptions(v3);
179: VecDuplicate(v3,&v4);
181: norm = 0.0;
182: for (i=0; i<10; i++) {
183: VecSetRandom(x,rdm);
184: MatMult(B,x,v5); /* v5 = B*x */
185: MatMultTranspose(P,v5,v3); /* v3 = Pt*B*x */
186: MatMult(C,x,v4); /* v4 = C*x */
187: VecNorm(v4,NORM_2,&norm_abs);
188: VecAXPY(v4,none,v3);
189: VecNorm(v4,NORM_2,&norm_tmp);
190: norm_tmp /= norm_abs;
191: if (norm_tmp > norm) norm = norm_tmp;
192: }
193: if (norm >= tol) {
194: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMultTr(), |v3 - v4|: %G\n",norm);
195: }
196: MatDestroy(&P);
197: MatDestroy(&C);
198: VecDestroy(&v3);
199: VecDestroy(&v4);
200: VecDestroy(&v5);
201: VecDestroy(&x);
202: }
204: /* Test MatPtAP() */
205: /*----------------------*/
206: if (Test_MatPtAP){
207: PetscInt PN;
208: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
209: MatGetSize(A,&M,&N);
210: MatGetLocalSize(A,&m,&n);
211: /* PetscPrintf(PETSC_COMM_SELF,"[%d] A: %d,%d, %d,%d\n",rank,m,n,M,N); */
213: PN = M/2;
214: nzp = (PetscInt)(0.1*PN); /* num of nozeros in each row of P */
215: MatCreate(PETSC_COMM_WORLD,&P);
216: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
217: MatSetType(P,MATAIJ);
218: MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
219: MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
220: for (i=0; i<nzp; i++){
221: PetscRandomGetValue(rdm,&a[i]);
222: }
223: for (i=0; i<M; i++){
224: for (j=0; j<nzp; j++){
225: PetscRandomGetValue(rdm,&rval);
226: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
227: }
228: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
229: }
230: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
231: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
233: /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
234: MatGetSize(P,&pM,&pN);
235: MatGetLocalSize(P,&pm,&pn);
236: /* PetscPrintf(PETSC_COMM_SELF," [%d] P, %d, %d, %d,%d\n",rank,pm,pn,pM,pN); */
237: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
239: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
240: alpha=1.0;
241: for (i=0; i<2; i++){
242: alpha -=0.1;
243: MatScale(A,alpha);
244: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
245: }
247: /* Create vector x that is compatible with P */
248: VecCreate(PETSC_COMM_WORLD,&x);
249: MatGetLocalSize(P,&m,&n);
250: VecSetSizes(x,n,PETSC_DECIDE);
251: VecSetFromOptions(x);
252:
253: VecCreate(PETSC_COMM_WORLD,&v3);
254: VecSetSizes(v3,n,PETSC_DECIDE);
255: VecSetFromOptions(v3);
256: VecDuplicate(v3,&v4);
258: norm = 0.0;
259: for (i=0; i<10; i++) {
260: VecSetRandom(x,rdm);
261: MatMult(P,x,v1);
262: MatMult(A,v1,v2); /* v2 = A*P*x */
264: MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
265: MatMult(C,x,v4); /* v3 = C*x */
266: VecNorm(v4,NORM_2,&norm_abs);
267: VecAXPY(v4,none,v3);
268: VecNorm(v4,NORM_2,&norm_tmp);
269: norm_tmp /= norm_abs;
270: if (norm_tmp > norm) norm = norm_tmp;
271: }
272: if (norm >= tol) {
273: PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v1 - v2|: %G\n",norm);
274: }
275:
276: MatDestroy(&A);
277: MatDestroy(&P);
278: MatDestroy(&C);
279: VecDestroy(&v3);
280: VecDestroy(&v4);
281: VecDestroy(&x);
282: } /* if (Test_MatPtAP) */
284: /* Destroy objects */
285: VecDestroy(&v1);
286: VecDestroy(&v2);
287: PetscRandomDestroy(&rdm);
288: PetscFree(idxn);
289:
290: MatDestroy(&A_save);
291: MatDestroy(&B);
293: PetscPreLoadEnd();
294: PetscFinalize();
296: return 0;
297: }