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., ex94 -f0 $D/small -f1 $D/small -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;
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;
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,RANDOM_DEFAULT,&rdm);
68: PetscOptionsGetReal(PETSC_NULL,"-fill",&fill,PETSC_NULL);
70: /* Test MatMatMult() */
71: /*-------------------*/
72: if (Test_MatMatMult){
73: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
74: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
75:
76: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
77: alpha=1.0;
78: for (i=0; i<2; i++){
79: alpha -=0.1;
80: MatScale(A,alpha);
81: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
82: }
83:
84: /* Create vector x that is compatible with B */
85: VecCreate(PETSC_COMM_WORLD,&x);
86: MatGetLocalSize(B,PETSC_NULL,&n);
87: VecSetSizes(x,n,PETSC_DECIDE);
88: VecSetFromOptions(x);
90: norm = 0.0;
91: for (i=0; i<10; i++) {
92: VecSetRandom(x,rdm);
93: MatMult(B,x,v1);
94: MatMult(A,v1,v2); /* v2 = A*B*x */
95: MatMult(C,x,v1); /* v1 = C*x */
96: VecAXPY(v1,none,v2);
97: VecNorm(v1,NORM_2,&norm_tmp);
98: if (norm_tmp > norm) norm = norm_tmp;
99: }
100: if (norm >= tol) {
101: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|: %G\n",norm);
102: }
103: MatDestroy(A);
104: MatDestroy(C);
105: VecDestroy(x);
106: } /* if (Test_MatMatMult) */
108: /* Test MatMatMultTranspose() */
109: /*----------------------------*/
110: if (size>1) Test_MatMatMultTr = PETSC_FALSE;
111: if (Test_MatMatMultTr){
112: PetscInt PN;
113: /* MatGetSize(B,&M,&N); */
114: PN = M/2;
115: nzp = 5; /* num of nonzeros in each row of P */
116: MatCreate(PETSC_COMM_WORLD,&P);
117: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
118: MatSetType(P,MATAIJ);
119: MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
120: MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
121: for (i=0; i<nzp; i++){
122: PetscRandomGetValue(rdm,&a[i]);
123: }
124: for (i=0; i<M; i++){
125: for (j=0; j<nzp; j++){
126: PetscRandomGetValue(rdm,&rval);
127: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
128: }
129: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
130: }
131: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
132: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
133:
134: MatMatMultTranspose(P,B,MAT_INITIAL_MATRIX,fill,&C);
136: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
137: alpha=1.0;
138: for (i=0; i<2; i++){
139: alpha -=0.1;
140: MatScale(P,alpha);
141: MatMatMultTranspose(P,B,MAT_REUSE_MATRIX,fill,&C);
142: }
144: /* Create vector x, v5 that are compatible with B */
145: VecCreate(PETSC_COMM_WORLD,&x);
146: MatGetLocalSize(B,&m,&n);
147: VecSetSizes(x,n,PETSC_DECIDE);
148: VecSetFromOptions(x);
150: VecCreate(PETSC_COMM_WORLD,&v5);
151: VecSetSizes(v5,m,PETSC_DECIDE);
152: VecSetFromOptions(v5);
153:
154: MatGetLocalSize(P,PETSC_NULL,&n);
155: VecCreate(PETSC_COMM_WORLD,&v3);
156: VecSetSizes(v3,n,PETSC_DECIDE);
157: VecSetFromOptions(v3);
158: VecDuplicate(v3,&v4);
160: norm = 0.0;
161: for (i=0; i<10; i++) {
162: VecSetRandom(x,rdm);
163: MatMult(B,x,v5); /* v5 = B*x */
164: MatMultTranspose(P,v5,v3); /* v3 = Pt*B*x */
165: MatMult(C,x,v4); /* v4 = C*x */
166: VecAXPY(v4,none,v3);
167: VecNorm(v4,NORM_2,&norm_tmp);
168: if (norm_tmp > norm) norm = norm_tmp;
169: }
170: if (norm >= tol) {
171: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMultTr(), |v3 - v4|: %G\n",norm);
172: }
173: MatDestroy(P);
174: MatDestroy(C);
175: VecDestroy(v3);
176: VecDestroy(v4);
177: VecDestroy(v5);
178: VecDestroy(x);
179: }
181: /* Test MatPtAP() */
182: /*----------------------*/
183: if (Test_MatPtAP){
184: PetscInt PN;
185: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
186: MatGetSize(A,&M,&N);
187: MatGetLocalSize(A,&m,&n);
188: /* PetscPrintf(PETSC_COMM_SELF,"[%d] A: %d,%d, %d,%d\n",rank,m,n,M,N); */
190: PN = M/2;
191: nzp = (PetscInt)(0.1*PN); /* num of nozeros in each row of P */
192: MatCreate(PETSC_COMM_WORLD,&P);
193: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
194: MatSetType(P,MATAIJ);
195: MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
196: MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
197: for (i=0; i<nzp; i++){
198: PetscRandomGetValue(rdm,&a[i]);
199: }
200: for (i=0; i<M; i++){
201: for (j=0; j<nzp; j++){
202: PetscRandomGetValue(rdm,&rval);
203: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
204: }
205: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
206: }
207: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
208: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
210: /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
211: MatGetSize(P,&pM,&pN);
212: MatGetLocalSize(P,&pm,&pn);
213: /* PetscPrintf(PETSC_COMM_SELF," [%d] P, %d, %d, %d,%d\n",rank,pm,pn,pM,pN); */
214: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
216: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
217: alpha=1.0;
218: for (i=0; i<2; i++){
219: alpha -=0.1;
220: MatScale(A,alpha);
221: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
222: }
224: /* Create vector x that is compatible with P */
225: VecCreate(PETSC_COMM_WORLD,&x);
226: MatGetLocalSize(P,&m,&n);
227: VecSetSizes(x,n,PETSC_DECIDE);
228: VecSetFromOptions(x);
229:
230: VecCreate(PETSC_COMM_WORLD,&v3);
231: VecSetSizes(v3,n,PETSC_DECIDE);
232: VecSetFromOptions(v3);
233: VecDuplicate(v3,&v4);
235: norm = 0.0;
236: for (i=0; i<10; i++) {
237: VecSetRandom(x,rdm);
238: MatMult(P,x,v1);
239: MatMult(A,v1,v2); /* v2 = A*P*x */
241: MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
242: MatMult(C,x,v4); /* v3 = C*x */
243: VecAXPY(v4,none,v3);
244: VecNorm(v4,NORM_2,&norm_tmp);
245: if (norm_tmp > norm) norm = norm_tmp;
246: }
247: if (norm >= tol) {
248: PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v1 - v2|: %G\n",norm);
249: }
250:
251: MatDestroy(A);
252: MatDestroy(P);
253: MatDestroy(C);
254: VecDestroy(v3);
255: VecDestroy(v4);
256: VecDestroy(x);
257: } /* if (Test_MatPtAP) */
259: /* Destroy objects */
260: VecDestroy(v1);
261: VecDestroy(v2);
262: PetscRandomDestroy(rdm);
263: PetscFree(idxn);
264:
265: MatDestroy(A_save);
266: MatDestroy(B);
268: PreLoadEnd();
269: PetscFinalize();
271: return 0;
272: }