Actual source code: ex94.c
petsc-3.3-p6 2013-02-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_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,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;
28: PetscInt pm,pn,pM,pN;
29: MatInfo info;
30:
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(PETSC_NULL,"-fill",&fill,PETSC_NULL);
37: /* Load the matrices A_save and B */
38: PetscOptionsGetString(PETSC_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(PETSC_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(PETSC_NULL,"-f2",file[2],128,&flg);
43: if (!flg) {
44: preload = PETSC_FALSE;
45: } else {
46: PetscOptionsGetString(PETSC_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 = (PetscInt)(0.1*M);
63: PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
64: a = (PetscScalar*)(idxn + nzp);
65:
66: /* Create vectors v1 and v2 that are compatible with A_save */
67: VecCreate(PETSC_COMM_WORLD,&v1);
68: MatGetLocalSize(A_save,&m,PETSC_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(PETSC_NULL,"-fill",&fill,PETSC_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,PETSC_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);
109: norm_tmp /= norm_abs;
110: if (norm_tmp > norm) norm = norm_tmp;
111: }
112: if (norm >= tol) {
113: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|: %G\n",norm);
114: }
116: VecDestroy(&x);
117: MatDestroy(&A);
119: /* Test MatDuplicate() of C */
120: MatDuplicate(C,MAT_COPY_VALUES,&C1);
121: MatDestroy(&C1);
122: MatDestroy(&C);
123: } /* if (Test_MatMatMult) */
125: /* Test MatTransposeMatMult() and MatMatTransposeMult() */
126: /*------------------------------------------------------*/
127: if (Test_MatMatTr){
128: /* Create P */
129: PetscInt PN,rstart,rend;
130: PN = M/2;
131: nzp = 5; /* num of nonzeros in each row of P */
132: MatCreate(PETSC_COMM_WORLD,&P);
133: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
134: MatSetType(P,MATAIJ);
135: MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
136: MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
137: MatGetOwnershipRange(P,&rstart,&rend);
138: for (i=0; i<nzp; i++){
139: PetscRandomGetValue(rdm,&a[i]);
140: }
141: for (i=rstart; i<rend; i++){
142: for (j=0; j<nzp; j++){
143: PetscRandomGetValue(rdm,&rval);
144: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
145: }
146: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
147: }
148: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
149: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
151: /* Create R = P^T */
152: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
153:
154: /* C = P^T*B */
155: MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,fill,&C);
156: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
157:
158: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
159: MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,fill,&C);
160:
161: /* Compare P^T*B and R*B */
162: MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
163: MatEqual(C,C1,&flg);
164: if (!flg){
165: /* Check norm of C1 = (-1.0)*C + C1 */
166: PetscReal nrm;
167: MatAXPY(C1,-1.0,C,DIFFERENT_NONZERO_PATTERN);
168: MatNorm(C1,NORM_INFINITY,&nrm);
169: if (nrm > 1.e-14){
170: PetscPrintf(PETSC_COMM_WORLD,"Error in MatTransposeMatMult(): %g\n",nrm);
171: }
172: }
173: MatDestroy(&C1);
174: MatDestroy(&C);
176: /* C = B*R^T */
177: if (size == 1){
178: MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,fill,&C);
179: MatSetOptionsPrefix(C,"matmatmulttr_"); /* enable '-matmatmulttr_' for matrix C */
180: MatGetInfo(C,MAT_GLOBAL_SUM,&info);
181:
182: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
183: MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,fill,&C);
184:
185: /* Check */
186: MatMatMult(B,P,MAT_INITIAL_MATRIX,fill,&C1);
187: MatEqual(C,C1,&flg);
188: if (!flg){
189: PetscPrintf(PETSC_COMM_WORLD,"Error in MatMatTransposeMult()\n");
190: }
191: MatDestroy(&C1);
192: MatDestroy(&C);
193: }
194: MatDestroy(&P);
195: MatDestroy(&R);
196: }
198: /* Test MatPtAP() */
199: /*----------------------*/
200: if (Test_MatPtAP){
201: PetscInt PN;
202: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
203: MatGetSize(A,&M,&N);
204: MatGetLocalSize(A,&m,&n);
205: /* PetscPrintf(PETSC_COMM_SELF,"[%d] A: %d,%d, %d,%d\n",rank,m,n,M,N); */
207: PN = M/2;
208: nzp = (PetscInt)(0.1*PN+1); /* num of nozeros in each row of P */
209: MatCreate(PETSC_COMM_WORLD,&P);
210: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
211: MatSetType(P,MATAIJ);
212: MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
213: MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
214: for (i=0; i<nzp; i++){
215: PetscRandomGetValue(rdm,&a[i]);
216: }
217: MatGetOwnershipRange(P,&rstart,&rend);
218: for (i=rstart; i<rend; i++){
219: for (j=0; j<nzp; j++){
220: PetscRandomGetValue(rdm,&rval);
221: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
222: }
223: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
224: }
225: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
226: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
228: /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
229: MatGetSize(P,&pM,&pN);
230: MatGetLocalSize(P,&pm,&pn);
231: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
232: /* if (!rank){PetscPrintf(PETSC_COMM_SELF," MatPtAP() is done, P, %d, %d, %d,%d\n",pm,pn,pM,pN);} */
234: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
235: alpha=1.0;
236: for (i=0; i<2; i++){
237: alpha -=0.1;
238: MatScale(A,alpha);
239: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
240: }
242: /* Test MatDuplicate() */
243: Mat Cdup;
244: MatDuplicate(C,MAT_COPY_VALUES,&Cdup);
245: MatDestroy(&Cdup);
247: if (size>1) Test_MatRARt = PETSC_FALSE;
248: /* Test MatRARt() */
249: if (Test_MatRARt){
250: Mat R, RARt;
251: PetscBool equal;
252: MatTranspose(P,MAT_INITIAL_MATRIX,&R);
253: MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);
254: MatEqual(C,RARt,&equal);
255: if (!equal) {
256: PetscReal norm;
257: MatAXPY(RARt,-1.0,C,DIFFERENT_NONZERO_PATTERN); /* RARt = -RARt + C */
258: MatNorm(RARt,NORM_FROBENIUS,&norm);
259: if (norm > 1.e-12) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"|PtAP - RARt| = %g",norm);
260: }
261: MatDestroy(&R);
262: MatDestroy(&RARt);
263: }
265: /* Create vector x that is compatible with P */
266: VecCreate(PETSC_COMM_WORLD,&x);
267: MatGetLocalSize(P,&m,&n);
268: VecSetSizes(x,n,PETSC_DECIDE);
269: VecSetFromOptions(x);
270:
271: VecCreate(PETSC_COMM_WORLD,&v3);
272: VecSetSizes(v3,n,PETSC_DECIDE);
273: VecSetFromOptions(v3);
274: VecDuplicate(v3,&v4);
276: norm = 0.0;
277: for (i=0; i<10; i++) {
278: VecSetRandom(x,rdm);
279: MatMult(P,x,v1);
280: MatMult(A,v1,v2); /* v2 = A*P*x */
282: MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
283: MatMult(C,x,v4); /* v3 = C*x */
284: VecNorm(v4,NORM_2,&norm_abs);
285: VecAXPY(v4,none,v3);
286: VecNorm(v4,NORM_2,&norm_tmp);
287: norm_tmp /= norm_abs;
288: if (norm_tmp > norm) norm = norm_tmp;
289: }
290: if (norm >= tol) {
291: PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v1 - v2|: %G\n",norm);
292: }
293:
294: MatDestroy(&A);
295: MatDestroy(&P);
296: MatDestroy(&C);
297: VecDestroy(&v3);
298: VecDestroy(&v4);
299: VecDestroy(&x);
300: }
302: /* Destroy objects */
303: VecDestroy(&v1);
304: VecDestroy(&v2);
305: PetscRandomDestroy(&rdm);
306: PetscFree(idxn);
307:
308: MatDestroy(&A_save);
309: MatDestroy(&B);
311: PetscPreLoadEnd();
312: PetscFinalize();
313: return 0;
314: }