Actual source code: ex107.c
1: static char help[] = "Tests PLAPACK interface.\n\n";
3: /* Usage:
4: mpirun -np 4 ex107 -M 50 -mat_plapack_nprows 2 -mat_plapack_npcols 2 -mat_plapack_nb 1
5: */
7: #include petscmat.h
10: int main(int argc,char **args)
11: {
12: Mat C,F,Cpetsc,Csymm;
13: Vec u,x,b,bpla;
15: PetscMPIInt rank,nproc;
16: PetscInt i,j,k,M = 10,m,n,nfact,nsolve,Istart,Iend,*im,*in;
17: PetscScalar *array,rval;
18: PetscReal norm;
19: IS perm,iperm;
20: MatFactorInfo info;
21: PetscRandom rand;
23: PetscInitialize(&argc,&args,(char *)0,help);
24: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
25: MPI_Comm_size(PETSC_COMM_WORLD, &nproc);
27: #ifdef PETSC_HAVE_PLAPACK
28: /* Test non-symmetric operations */
29: /*-------------------------------*/
30: /* Create a Plapack dense matrix C */
31: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
32: MatCreate(PETSC_COMM_WORLD,&C);
33: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,M,M);
34: MatSetType(C,MATPLAPACK);
35: MatSetFromOptions(C);
37: /* Create vectors */
38: MatGetLocalSize(C,&m,&n);
39: if (m != n) SETERRQ2(PETSC_ERR_ARG_WRONG,"Matrix local size m %d must equal n %d",m,n);
40: /* printf("[%d] C - local size m: %d\n",rank,m); */
41: VecCreate(PETSC_COMM_WORLD,&x);
42: VecSetSizes(x,m,PETSC_DECIDE);
43: VecSetFromOptions(x);
44: VecDuplicate(x,&b);
45: VecDuplicate(x,&bpla);
46: VecDuplicate(x,&u); /* save the true solution */
48: /* Create a petsc dense matrix Cpetsc */
49: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
50: MatCreate(PETSC_COMM_WORLD,&Cpetsc);
51: MatSetSizes(Cpetsc,m,m,M,M);
52: MatSetType(Cpetsc,MATDENSE);
53: MatMPIDenseSetPreallocation(Cpetsc,PETSC_NULL);
54: MatSetFromOptions(Cpetsc);
55: MatSetOption(Cpetsc,MAT_COLUMN_ORIENTED);
56: MatSetOption(C,MAT_COLUMN_ORIENTED);
58: /* Assembly */
59: /* PLAPACK doesn't support INSERT_VALUES mode, zero all entries before calling MatSetValues() */
60: MatZeroEntries(C);
61: MatZeroEntries(Cpetsc);
62: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
63: PetscRandomSetFromOptions(rand);
64: MatGetOwnershipRange(C,&Istart,&Iend);
65: printf(" [%d] C m: %d, Istart/end: %d %d\n",rank,m,Istart,Iend);
66:
67: PetscMalloc((m*M+1)*sizeof(PetscScalar),&array);
68: PetscMalloc((m+M+1)*sizeof(PetscInt),&im);
69: in = im + m;
70: k = 0;
71: for (j=0; j<M; j++){ /* column oriented! */
72: in[j] = j;
73: for (i=0; i<m; i++){
74: im[i] = i+Istart;
75: PetscRandomGetValue(rand,&rval);
76: array[k++] = rval;
77: }
78: }
79: MatSetValues(Cpetsc,m,im,M,in,array,ADD_VALUES);
80: MatSetValues(C,m,im,M,in,array,ADD_VALUES);
81: PetscFree(array);
82: PetscFree(im);
84: MatAssemblyBegin(Cpetsc,MAT_FINAL_ASSEMBLY);
85: MatAssemblyEnd(Cpetsc,MAT_FINAL_ASSEMBLY);
86: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
87: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
88: /*
89: if (!rank) {printf("main, Cpetsc: \n");}
90: MatView(Cpetsc,PETSC_VIEWER_STDOUT_WORLD);
91: */
92: MatGetOrdering(C,MATORDERING_NATURAL,&perm,&iperm);
94: /* Test nonsymmetric MatMult() */
95: VecGetArray(x,&array);
96: for (i=0; i<m; i++){
97: PetscRandomGetValue(rand,&rval);
98: array[i] = rval;
99: }
100: VecRestoreArray(x,&array);
101:
102: MatMult(Cpetsc,x,b);
103: MatMult(C,x,bpla);
104: VecAXPY(bpla,-1.0,b);
105: VecNorm(bpla,NORM_2,&norm);
106: if (norm > 1.e-12 && !rank){
107: PetscPrintf(PETSC_COMM_SELF,"Nonsymmetric MatMult_Plapack error: |b_pla - b|= %g\n",norm);
108: }
110: /* Test LU Factorization */
111: MatLUFactorSymbolic(C,perm,iperm,&info,&F);
112: for (nfact = 0; nfact < 2; nfact++){
113: if (!rank) printf(" LU nfact %d\n",nfact);
114: if (nfact>0){ /* change matrix value for testing repeated MatLUFactorNumeric() */
115: if (!rank){
116: i = j = 0;
117: rval = nfact;
118: MatSetValues(Cpetsc,1,&i,1,&j,&rval,ADD_VALUES);
119: MatSetValues(C,1,&i,1,&j,&rval,ADD_VALUES);
120: } else { /* PLAPACK seems requiring all processors call MatSetValues(), so we add 0.0 on processesses with rank>0! */
121: i = j = 0;
122: rval = 0.0;
123: MatSetValues(Cpetsc,1,&i,1,&j,&rval,ADD_VALUES);
124: MatSetValues(C,1,&i,1,&j,&rval,ADD_VALUES);
125: }
126: MatAssemblyBegin(Cpetsc,MAT_FINAL_ASSEMBLY);
127: MatAssemblyEnd(Cpetsc,MAT_FINAL_ASSEMBLY);
128: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
129: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
130: }
131: MatLUFactorNumeric(C,&info,&F);
133: /* Test MatSolve() */
134: for (nsolve = 0; nsolve < 2; nsolve++){
135: VecGetArray(x,&array);
136: for (i=0; i<m; i++){
137: PetscRandomGetValue(rand,&rval);
138: array[i] = rval;
139: /* array[i] = rank + 1; */
140: }
141: VecRestoreArray(x,&array);
142: VecCopy(x,u);
143: MatMult(C,x,b);
144: MatSolve(F,b,x);
146: /* Check the error */
147: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
148: VecNorm(u,NORM_2,&norm);
149: if (!rank){
150: PetscPrintf(PETSC_COMM_SELF,"Norm of error %A\n",norm);
151: }
152: }
153: }
154: MatDestroy(F);
156: /* Test repeated MatLUFactorSymbolic() -- fail??? */
157: /*
158: MatLUFactorSymbolic(C,perm,iperm,&info,&F);
159: MatLUFactorNumeric(C,&info,&F);
160: MatDestroy(F);
161: */
162:
163: /* Test non-symmetric operations */
164: /*-------------------------------*/
165: /* Create a symmetric Plapack dense matrix Csymm */
166: MatCreate(PETSC_COMM_WORLD,&Csymm);
167: MatSetSizes(Csymm,PETSC_DECIDE,PETSC_DECIDE,M,M);
168: MatSetType(Csymm,MATPLAPACK);
169: MatSetFromOptions(Csymm);
170: MatSetOption(Csymm,MAT_COLUMN_ORIENTED);
171: MatSetOption(Csymm,MAT_SYMMETRIC);
172: MatSetOption(Csymm,MAT_SYMMETRY_ETERNAL);
174: MatZeroEntries(Csymm);
175: MatZeroEntries(Cpetsc);
176: for (i=Istart; i<Iend; i++){
177: for (j=0; j<=i; j++){
178: PetscRandomGetValue(rand,&rval);
179: MatSetValues(Cpetsc,1,&i,1,&j,&rval,ADD_VALUES);
180: MatSetValues(Csymm,1,&i,1,&j,&rval,ADD_VALUES);
181: if (j<i){
182: /* Although PLAPACK only requires lower triangular entries, we must add all the entries.
183: MatSetValues_Plapack() will ignore the upper triangular entries AFTER an index map! */
184: MatSetValues(Cpetsc,1,&j,1,&i,&rval,ADD_VALUES);
185: MatSetValues(Csymm,1,&j,1,&i,&rval,ADD_VALUES);
186: }
187: }
188: }
189: MatAssemblyBegin(Cpetsc,MAT_FINAL_ASSEMBLY);
190: MatAssemblyEnd(Cpetsc,MAT_FINAL_ASSEMBLY);
191: MatAssemblyBegin(Csymm,MAT_FINAL_ASSEMBLY);
192: MatAssemblyEnd(Csymm,MAT_FINAL_ASSEMBLY);
194: /* Test symmetric MatMult() */
195: VecGetArray(x,&array);
196: for (i=0; i<m; i++){
197: PetscRandomGetValue(rand,&rval);
198: array[i] = rval;
199: }
200: VecRestoreArray(x,&array);
201:
202: MatMult(Cpetsc,x,b);
203: MatMult(Csymm,x,bpla);
204: VecAXPY(bpla,-1.0,b);
205: VecNorm(bpla,NORM_2,&norm);
206: if (norm > 1.e-12 && !rank){
207: PetscPrintf(PETSC_COMM_SELF,"Symmetric MatMult_Plapack error: |b_pla - b|= %g\n",norm);
208: }
210: /* Test Cholesky Factorization */
211: MatShift(Csymm,M); /* make Csymm positive definite */
212: MatCholeskyFactorSymbolic(Csymm,perm,&info,&F);
213: for (nfact = 0; nfact < 2; nfact++){
214: if (!rank) printf(" Cholesky nfact %d\n",nfact);
215: MatCholeskyFactorNumeric(Csymm,&info,&F);
217: /* Test MatSolve() */
218: for (nsolve = 0; nsolve < 2; nsolve++){
219: VecGetArray(x,&array);
220: for (i=0; i<m; i++){
221: PetscRandomGetValue(rand,&rval);
222: array[i] = rval;
223: }
224: VecRestoreArray(x,&array);
225: VecCopy(x,u);
226: MatMult(Csymm,x,b);
227: MatSolve(F,b,x);
229: /* Check the error */
230: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
231: VecNorm(u,NORM_2,&norm);
232: if (!rank){
233: PetscPrintf(PETSC_COMM_SELF,"Norm of error %A\n",norm);
234: }
235: }
236: }
237: MatDestroy(F);
239: /* Test repeated MatFactorSymbolic() -- fail??? */
240: /*
241: MatCholeskyFactorSymbolic(Csymm,perm,&info,&F);
242: MatCholeskyFactorNumeric(Csymm,&info,&F);
243: MatDestroy(F);
244: */
246: /* Free data structures */
247: ISDestroy(perm);
248: ISDestroy(iperm);
249:
250: PetscRandomDestroy(rand);
251: VecDestroy(x);
252: VecDestroy(b);
253: VecDestroy(bpla);
254: VecDestroy(u);
255:
256: MatDestroy(Cpetsc);
257: MatDestroy(C);/* works! */
258: MatDestroy(Csymm);
259: /* MatDestroy(C); */
261: #else
262: if (!rank) printf("This example needs PLAPLAPACK\n");
263: #endif
264: PetscFinalize();
265: return 0;
266: }