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: /* e.g., mpiexec -n 3 ./ex94 -f0 $D/medium -f1 $D/medium -f2 $D/arco1 -f3 $D/arco1 */
7: #include petscmat.h
11: int main(int argc,char **args)
12: {
13: Mat A,A_save,B,P,C,C1;
14: Vec x,v1,v2;
15: PetscViewer viewer;
17: PetscMPIInt size,rank;
18: PetscInt i,m,n,j,*idxn,M,N,nzp;
19: PetscReal norm,norm_tmp,tol=1.e-8,fill=4.0;
20: PetscRandom rdm;
21: char file[4][128];
22: PetscTruth flg,preload = PETSC_TRUE;
23: PetscScalar *a,rval,alpha,none = -1.0;
24: PetscTruth Test_MatMatMult=PETSC_TRUE,Test_MatMatMultTr=PETSC_TRUE;
25: Vec v3,v4,v5;
26: PetscInt pm,pn,pM,pN;
27: PetscTruth Test_MatPtAP=PETSC_TRUE;
28:
29: PetscInitialize(&argc,&args,(char *)0,help);
30: MPI_Comm_size(PETSC_COMM_WORLD,&size);
31: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
33: /* Load the matrices A and B */
34: PetscOptionsGetString(PETSC_NULL,"-f0",file[0],127,&flg);
35: if (!flg) SETERRQ(1,"Must indicate a file name for small matrix A with the -f0 option.");
36: PetscOptionsGetString(PETSC_NULL,"-f1",file[1],127,&flg);
37: if (!flg) SETERRQ(1,"Must indicate a file name for small matrix B with the -f1 option.");
38: PetscOptionsGetString(PETSC_NULL,"-f2",file[2],127,&flg);
39: if (!flg) {
40: preload = PETSC_FALSE;
41: } else {
42: PetscOptionsGetString(PETSC_NULL,"-f3",file[3],127,&flg);
43: if (!flg) SETERRQ(1,"Must indicate a file name for test matrix B with the -f3 option.");
44: }
46: PreLoadBegin(preload,"Load system");
47: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PreLoadIt],FILE_MODE_READ,&viewer);
48: MatLoad(viewer,MATAIJ,&A_save);
49: PetscViewerDestroy(viewer);
51: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PreLoadIt+1],FILE_MODE_READ,&viewer);
52: MatLoad(viewer,MATAIJ,&B);
53: PetscViewerDestroy(viewer);
55: MatGetSize(B,&M,&N);
56: nzp = (PetscInt)(0.1*M);
57: PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
58: a = (PetscScalar*)(idxn + nzp);
59:
60: /* Create vectors v1 and v2 that are compatible with A_save */
61: VecCreate(PETSC_COMM_WORLD,&v1);
62: MatGetLocalSize(A_save,&m,PETSC_NULL);
63: VecSetSizes(v1,m,PETSC_DECIDE);
64: VecSetFromOptions(v1);
65: VecDuplicate(v1,&v2);
67: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
68: PetscRandomSetFromOptions(rdm);
69: PetscOptionsGetReal(PETSC_NULL,"-fill",&fill,PETSC_NULL);
71: /* Test MatMatMult() */
72: /*-------------------*/
73: if (Test_MatMatMult){
74: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
75: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
77: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
78: alpha=1.0;
79: for (i=0; i<2; i++){
80: alpha -=0.1;
81: MatScale(A,alpha);
82: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
83: }
85: /* Create vector x that is compatible with B */
86: VecCreate(PETSC_COMM_WORLD,&x);
87: MatGetLocalSize(B,PETSC_NULL,&n);
88: VecSetSizes(x,n,PETSC_DECIDE);
89: VecSetFromOptions(x);
91: norm = 0.0;
92: for (i=0; i<10; i++) {
93: VecSetRandom(x,rdm);
94: MatMult(B,x,v1);
95: MatMult(A,v1,v2); /* v2 = A*B*x */
96: MatMult(C,x,v1); /* v1 = C*x */
97: VecAXPY(v1,none,v2);
98: VecNorm(v1,NORM_2,&norm_tmp);
99: if (norm_tmp > norm) norm = norm_tmp;
100: }
101: if (norm >= tol) {
102: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|: %G\n",norm);
103: }
105: VecDestroy(x);
106: MatDestroy(A);
108: /* Test MatDuplicate() of C */
109: MatDuplicate(C,MAT_COPY_VALUES,&C1);
110: MatDestroy(C1);
111: MatDestroy(C);
112: } /* if (Test_MatMatMult) */
114: /* Test MatMatMultTranspose() */
115: /*----------------------------*/
116: if (size>1) Test_MatMatMultTr = PETSC_FALSE;
117: if (Test_MatMatMultTr){
118: PetscInt PN;
119: /* MatGetSize(B,&M,&N); */
120: PN = M/2;
121: nzp = 5; /* num of nonzeros in each row of P */
122: MatCreate(PETSC_COMM_WORLD,&P);
123: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
124: MatSetType(P,MATAIJ);
125: MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
126: MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
127: for (i=0; i<nzp; i++){
128: PetscRandomGetValue(rdm,&a[i]);
129: }
130: for (i=0; i<M; i++){
131: for (j=0; j<nzp; j++){
132: PetscRandomGetValue(rdm,&rval);
133: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
134: }
135: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
136: }
137: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
138: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
139:
140: MatMatMultTranspose(P,B,MAT_INITIAL_MATRIX,fill,&C);
142: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
143: alpha=1.0;
144: for (i=0; i<2; i++){
145: alpha -=0.1;
146: MatScale(P,alpha);
147: MatMatMultTranspose(P,B,MAT_REUSE_MATRIX,fill,&C);
148: }
150: /* Create vector x, v5 that are compatible with B */
151: VecCreate(PETSC_COMM_WORLD,&x);
152: MatGetLocalSize(B,&m,&n);
153: VecSetSizes(x,n,PETSC_DECIDE);
154: VecSetFromOptions(x);
156: VecCreate(PETSC_COMM_WORLD,&v5);
157: VecSetSizes(v5,m,PETSC_DECIDE);
158: VecSetFromOptions(v5);
159:
160: MatGetLocalSize(P,PETSC_NULL,&n);
161: VecCreate(PETSC_COMM_WORLD,&v3);
162: VecSetSizes(v3,n,PETSC_DECIDE);
163: VecSetFromOptions(v3);
164: VecDuplicate(v3,&v4);
166: norm = 0.0;
167: for (i=0; i<10; i++) {
168: VecSetRandom(x,rdm);
169: MatMult(B,x,v5); /* v5 = B*x */
170: MatMultTranspose(P,v5,v3); /* v3 = Pt*B*x */
171: MatMult(C,x,v4); /* v4 = C*x */
172: VecAXPY(v4,none,v3);
173: VecNorm(v4,NORM_2,&norm_tmp);
174: if (norm_tmp > norm) norm = norm_tmp;
175: }
176: if (norm >= tol) {
177: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMultTr(), |v3 - v4|: %G\n",norm);
178: }
179: MatDestroy(P);
180: MatDestroy(C);
181: VecDestroy(v3);
182: VecDestroy(v4);
183: VecDestroy(v5);
184: VecDestroy(x);
185: }
187: /* Test MatPtAP() */
188: /*----------------------*/
189: if (Test_MatPtAP){
190: PetscInt PN;
191: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
192: MatGetSize(A,&M,&N);
193: MatGetLocalSize(A,&m,&n);
194: /* PetscPrintf(PETSC_COMM_SELF,"[%d] A: %d,%d, %d,%d\n",rank,m,n,M,N); */
196: PN = M/2;
197: nzp = (PetscInt)(0.1*PN); /* num of nozeros in each row of P */
198: MatCreate(PETSC_COMM_WORLD,&P);
199: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
200: MatSetType(P,MATAIJ);
201: MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
202: MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
203: for (i=0; i<nzp; i++){
204: PetscRandomGetValue(rdm,&a[i]);
205: }
206: for (i=0; i<M; i++){
207: for (j=0; j<nzp; j++){
208: PetscRandomGetValue(rdm,&rval);
209: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
210: }
211: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
212: }
213: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
214: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
216: /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
217: MatGetSize(P,&pM,&pN);
218: MatGetLocalSize(P,&pm,&pn);
219: /* PetscPrintf(PETSC_COMM_SELF," [%d] P, %d, %d, %d,%d\n",rank,pm,pn,pM,pN); */
220: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
222: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
223: alpha=1.0;
224: for (i=0; i<2; i++){
225: alpha -=0.1;
226: MatScale(A,alpha);
227: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
228: }
230: /* Create vector x that is compatible with P */
231: VecCreate(PETSC_COMM_WORLD,&x);
232: MatGetLocalSize(P,&m,&n);
233: VecSetSizes(x,n,PETSC_DECIDE);
234: VecSetFromOptions(x);
235:
236: VecCreate(PETSC_COMM_WORLD,&v3);
237: VecSetSizes(v3,n,PETSC_DECIDE);
238: VecSetFromOptions(v3);
239: VecDuplicate(v3,&v4);
241: norm = 0.0;
242: for (i=0; i<10; i++) {
243: VecSetRandom(x,rdm);
244: MatMult(P,x,v1);
245: MatMult(A,v1,v2); /* v2 = A*P*x */
247: MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
248: MatMult(C,x,v4); /* v3 = C*x */
249: VecAXPY(v4,none,v3);
250: VecNorm(v4,NORM_2,&norm_tmp);
251: if (norm_tmp > norm) norm = norm_tmp;
252: }
253: if (norm >= tol) {
254: PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v1 - v2|: %G\n",norm);
255: }
256:
257: MatDestroy(A);
258: MatDestroy(P);
259: MatDestroy(C);
260: VecDestroy(v3);
261: VecDestroy(v4);
262: VecDestroy(x);
263: } /* if (Test_MatPtAP) */
265: /* Destroy objects */
266: VecDestroy(v1);
267: VecDestroy(v2);
268: PetscRandomDestroy(rdm);
269: PetscFree(idxn);
270:
271: MatDestroy(A_save);
272: MatDestroy(B);
274: PreLoadEnd();
275: PetscFinalize();
277: return 0;
278: }