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,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;
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);
76:
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: }
84:
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: /* Test MatDuplicate() of C */
106: MatDuplicate(C,MAT_COPY_VALUES,&C1);
107: MatDestroy(C1);
109: /* Test MatConvert() of C to its inherited matrix classes */
110: #if defined(PETSC_HAVE_SUPERLU)
111: if (size == 1){
112: MatConvert(C,MATSUPERLU,MAT_INITIAL_MATRIX,&C1);
113: MatDestroy(C1);
114: }
115: #endif
116: #if defined(PETSC_HAVE_SUPERLU_DIST)
117: MatConvert(C,MATSUPERLU_DIST,MAT_INITIAL_MATRIX,&C1);
118: MatDestroy(C1);
119: #endif
120: #if defined(PETSC_HAVE_MUMPS)
121: MatConvert(C,MATAIJMUMPS,MAT_INITIAL_MATRIX,&C1);
122: MatDestroy(C1);
123: #endif
124: #if defined(PETSC_HAVE_SPOOLES)
125: MatConvert(C,MATAIJSPOOLES,MAT_INITIAL_MATRIX,&C1);
126: MatDestroy(C1);
127: #endif
128: MatDestroy(A);
129: MatDestroy(C);
130: VecDestroy(x);
131: } /* if (Test_MatMatMult) */
133: /* Test MatMatMultTranspose() */
134: /*----------------------------*/
135: if (size>1) Test_MatMatMultTr = PETSC_FALSE;
136: if (Test_MatMatMultTr){
137: PetscInt PN;
138: /* MatGetSize(B,&M,&N); */
139: PN = M/2;
140: nzp = 5; /* num of nonzeros in each row of P */
141: MatCreate(PETSC_COMM_WORLD,&P);
142: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
143: MatSetType(P,MATAIJ);
144: MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
145: MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
146: for (i=0; i<nzp; i++){
147: PetscRandomGetValue(rdm,&a[i]);
148: }
149: for (i=0; i<M; i++){
150: for (j=0; j<nzp; j++){
151: PetscRandomGetValue(rdm,&rval);
152: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
153: }
154: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
155: }
156: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
157: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
158:
159: MatMatMultTranspose(P,B,MAT_INITIAL_MATRIX,fill,&C);
161: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
162: alpha=1.0;
163: for (i=0; i<2; i++){
164: alpha -=0.1;
165: MatScale(P,alpha);
166: MatMatMultTranspose(P,B,MAT_REUSE_MATRIX,fill,&C);
167: }
169: /* Create vector x, v5 that are compatible with B */
170: VecCreate(PETSC_COMM_WORLD,&x);
171: MatGetLocalSize(B,&m,&n);
172: VecSetSizes(x,n,PETSC_DECIDE);
173: VecSetFromOptions(x);
175: VecCreate(PETSC_COMM_WORLD,&v5);
176: VecSetSizes(v5,m,PETSC_DECIDE);
177: VecSetFromOptions(v5);
178:
179: MatGetLocalSize(P,PETSC_NULL,&n);
180: VecCreate(PETSC_COMM_WORLD,&v3);
181: VecSetSizes(v3,n,PETSC_DECIDE);
182: VecSetFromOptions(v3);
183: VecDuplicate(v3,&v4);
185: norm = 0.0;
186: for (i=0; i<10; i++) {
187: VecSetRandom(x,rdm);
188: MatMult(B,x,v5); /* v5 = B*x */
189: MatMultTranspose(P,v5,v3); /* v3 = Pt*B*x */
190: MatMult(C,x,v4); /* v4 = C*x */
191: VecAXPY(v4,none,v3);
192: VecNorm(v4,NORM_2,&norm_tmp);
193: if (norm_tmp > norm) norm = norm_tmp;
194: }
195: if (norm >= tol) {
196: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMultTr(), |v3 - v4|: %G\n",norm);
197: }
198: MatDestroy(P);
199: MatDestroy(C);
200: VecDestroy(v3);
201: VecDestroy(v4);
202: VecDestroy(v5);
203: VecDestroy(x);
204: }
206: /* Test MatPtAP() */
207: /*----------------------*/
208: if (Test_MatPtAP){
209: PetscInt PN;
210: MatDuplicate(A_save,MAT_COPY_VALUES,&A);
211: MatGetSize(A,&M,&N);
212: MatGetLocalSize(A,&m,&n);
213: /* PetscPrintf(PETSC_COMM_SELF,"[%d] A: %d,%d, %d,%d\n",rank,m,n,M,N); */
215: PN = M/2;
216: nzp = (PetscInt)(0.1*PN); /* num of nozeros in each row of P */
217: MatCreate(PETSC_COMM_WORLD,&P);
218: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
219: MatSetType(P,MATAIJ);
220: MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
221: MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
222: for (i=0; i<nzp; i++){
223: PetscRandomGetValue(rdm,&a[i]);
224: }
225: for (i=0; i<M; i++){
226: for (j=0; j<nzp; j++){
227: PetscRandomGetValue(rdm,&rval);
228: idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
229: }
230: MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
231: }
232: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
233: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
235: /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
236: MatGetSize(P,&pM,&pN);
237: MatGetLocalSize(P,&pm,&pn);
238: /* PetscPrintf(PETSC_COMM_SELF," [%d] P, %d, %d, %d,%d\n",rank,pm,pn,pM,pN); */
239: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
241: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
242: alpha=1.0;
243: for (i=0; i<2; i++){
244: alpha -=0.1;
245: MatScale(A,alpha);
246: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
247: }
249: /* Create vector x that is compatible with P */
250: VecCreate(PETSC_COMM_WORLD,&x);
251: MatGetLocalSize(P,&m,&n);
252: VecSetSizes(x,n,PETSC_DECIDE);
253: VecSetFromOptions(x);
254:
255: VecCreate(PETSC_COMM_WORLD,&v3);
256: VecSetSizes(v3,n,PETSC_DECIDE);
257: VecSetFromOptions(v3);
258: VecDuplicate(v3,&v4);
260: norm = 0.0;
261: for (i=0; i<10; i++) {
262: VecSetRandom(x,rdm);
263: MatMult(P,x,v1);
264: MatMult(A,v1,v2); /* v2 = A*P*x */
266: MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
267: MatMult(C,x,v4); /* v3 = C*x */
268: VecAXPY(v4,none,v3);
269: VecNorm(v4,NORM_2,&norm_tmp);
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: PreLoadEnd();
294: PetscFinalize();
296: return 0;
297: }