Actual source code: ex107.c
petsc-3.3-p6 2013-02-11
1: static char help[] = "Tests PLAPACK interface.\n\n";
3: /* Usage:
4: mpiexec -n 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,nfact,nsolve,Istart,Iend,*im,*in,start,end;
17: PetscScalar *array,rval;
18: PetscReal norm,tol=1.e-12;
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: /* Test non-symmetric operations */
28: /*-------------------------------*/
29: /* Create a Plapack dense matrix C */
30: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
31: MatCreate(PETSC_COMM_WORLD,&C);
32: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,M,M);
33: MatSetType(C,MATDENSE);
34: MatSetFromOptions(C);
35: MatSetUp(C);
37: /* Create vectors */
38: MatGetOwnershipRange(C,&start,&end);
39: m = end - start;
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: MatSetUp(Cpetsc);
57: MatSetOption(Cpetsc,MAT_ROW_ORIENTED,PETSC_FALSE);
58: MatSetOption(C,MAT_ROW_ORIENTED,PETSC_FALSE);
60: /* Assembly */
61: /* PLAPACK doesn't support INSERT_VALUES mode, zero all entries before calling MatSetValues() */
62: MatZeroEntries(C);
63: MatZeroEntries(Cpetsc);
64: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
65: PetscRandomSetFromOptions(rand);
66: MatGetOwnershipRange(C,&Istart,&Iend);
67: /* printf(" [%d] C m: %d, Istart/end: %d %d\n",rank,m,Istart,Iend); */
68:
69: PetscMalloc((m*M+1)*sizeof(PetscScalar),&array);
70: PetscMalloc2(m,PetscInt,&im,M,PetscInt,&in);
71: k = 0;
72: for (j=0; j<M; j++){ /* column oriented! */
73: in[j] = j;
74: for (i=0; i<m; i++){
75: im[i] = i+Istart;
76: PetscRandomGetValue(rand,&rval);
77: array[k++] = rval;
78: }
79: }
80: MatSetValues(Cpetsc,m,im,M,in,array,ADD_VALUES);
81: MatSetValues(C,m,im,M,in,array,ADD_VALUES);
82: PetscFree(array);
83: PetscFree2(im,in);
85: MatAssemblyBegin(Cpetsc,MAT_FINAL_ASSEMBLY);
86: MatAssemblyEnd(Cpetsc,MAT_FINAL_ASSEMBLY);
87: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
89: /*
90: if (!rank) {printf("main, Cpetsc: \n");}
91: MatView(Cpetsc,PETSC_VIEWER_STDOUT_WORLD);
92: */
93: MatGetOrdering(C,MATORDERINGNATURAL,&perm,&iperm);
95: /* Test nonsymmetric MatMult() */
96: VecGetArray(x,&array);
97: for (i=0; i<m; i++){
98: PetscRandomGetValue(rand,&rval);
99: array[i] = rval;
100: }
101: VecRestoreArray(x,&array);
102:
103: MatMult(Cpetsc,x,b);
104: MatMult(C,x,bpla);
105: VecAXPY(bpla,-1.0,b);
106: VecNorm(bpla,NORM_2,&norm);
107: if (norm > 1.e-12 && !rank){
108: PetscPrintf(PETSC_COMM_SELF,"Nonsymmetric MatMult_Plapack error: |b_pla - b|= %g\n",norm);
109: }
111: /* Test LU Factorization */
112: if (nproc == 1){
113: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
114: } else {
115: MatGetFactor(C,MATSOLVERPLAPACK,MAT_FACTOR_LU,&F);
116: }
117: MatLUFactorSymbolic(F,C,perm,iperm,&info);
118: for (nfact = 0; nfact < 2; nfact++){
119: if (!rank) printf(" LU nfact %d\n",nfact);
120: if (nfact>0){ /* change matrix value for testing repeated MatLUFactorNumeric() */
121: if (!rank){
122: i = j = 0;
123: rval = nfact;
124: MatSetValues(Cpetsc,1,&i,1,&j,&rval,ADD_VALUES);
125: MatSetValues(C,1,&i,1,&j,&rval,ADD_VALUES);
126: } else { /* PLAPACK seems requiring all processors call MatSetValues(), so we add 0.0 on processesses with rank>0! */
127: i = j = 0;
128: rval = 0.0;
129: MatSetValues(Cpetsc,1,&i,1,&j,&rval,ADD_VALUES);
130: MatSetValues(C,1,&i,1,&j,&rval,ADD_VALUES);
131: }
132: MatAssemblyBegin(Cpetsc,MAT_FINAL_ASSEMBLY);
133: MatAssemblyEnd(Cpetsc,MAT_FINAL_ASSEMBLY);
134: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
135: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
136: }
137: MatLUFactorNumeric(F,C,&info);
139: /* Test MatSolve() */
140: for (nsolve = 0; nsolve < 2; nsolve++){
141: VecGetArray(x,&array);
142: for (i=0; i<m; i++){
143: PetscRandomGetValue(rand,&rval);
144: array[i] = rval;
145: /* array[i] = rank + 1; */
146: }
147: VecRestoreArray(x,&array);
148: VecCopy(x,u);
149: MatMult(C,x,b);
150: MatSolve(F,b,x);
152: /* Check the error */
153: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
154: VecNorm(u,NORM_2,&norm);
155: if (norm > tol){
156: if (!rank){
157: PetscPrintf(PETSC_COMM_SELF,"Error: Norm of error %g, LU nfact %d\n",norm,nfact);
158: }
159: }
160: }
161: }
162: MatDestroy(&F);
163:
164: /* Test non-symmetric operations */
165: /*-------------------------------*/
166: /* Create a symmetric Plapack dense matrix Csymm */
167: MatCreate(PETSC_COMM_WORLD,&Csymm);
168: MatSetSizes(Csymm,PETSC_DECIDE,PETSC_DECIDE,M,M);
169: MatSetType(Csymm,MATDENSE);
170: MatSetFromOptions(Csymm);
171: MatSetUp(Csymm);
173: MatSetOption(Csymm,MAT_ROW_ORIENTED,PETSC_FALSE);
174: MatSetOption(Csymm,MAT_SYMMETRIC,PETSC_TRUE);
175: MatSetOption(Csymm,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
177: MatZeroEntries(Csymm);
178: MatZeroEntries(Cpetsc);
179: for (i=Istart; i<Iend; i++){
180: for (j=0; j<=i; j++){
181: PetscRandomGetValue(rand,&rval);
182: MatSetValues(Cpetsc,1,&i,1,&j,&rval,ADD_VALUES);
183: MatSetValues(Csymm,1,&i,1,&j,&rval,ADD_VALUES);
184: if (j<i){
185: /* Although PLAPACK only requires lower triangular entries, we must add all the entries.
186: MatSetValues_Plapack() will ignore the upper triangular entries AFTER an index map! */
187: MatSetValues(Cpetsc,1,&j,1,&i,&rval,ADD_VALUES);
188: MatSetValues(Csymm,1,&j,1,&i,&rval,ADD_VALUES);
189: }
190: }
191: }
192: MatAssemblyBegin(Cpetsc,MAT_FINAL_ASSEMBLY);
193: MatAssemblyEnd(Cpetsc,MAT_FINAL_ASSEMBLY);
194: MatAssemblyBegin(Csymm,MAT_FINAL_ASSEMBLY);
195: MatAssemblyEnd(Csymm,MAT_FINAL_ASSEMBLY);
197: /* Test symmetric MatMult() */
198: VecGetArray(x,&array);
199: for (i=0; i<m; i++){
200: PetscRandomGetValue(rand,&rval);
201: array[i] = rval;
202: }
203: VecRestoreArray(x,&array);
204:
205: MatMult(Cpetsc,x,b);
206: MatMult(Csymm,x,bpla);
207: VecAXPY(bpla,-1.0,b);
208: VecNorm(bpla,NORM_2,&norm);
209: if (norm > 1.e-12 && !rank){
210: PetscPrintf(PETSC_COMM_SELF,"Symmetric MatMult_Plapack error: |b_pla - b|= %g\n",norm);
211: }
213: /* Test Cholesky Factorization */
214: MatShift(Csymm,M); /* make Csymm positive definite */
215: if (nproc == 1){
216: MatGetFactor(Csymm,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F);
217: } else {
218: MatGetFactor(Csymm,MATSOLVERPLAPACK,MAT_FACTOR_CHOLESKY,&F);
219: }
220: MatCholeskyFactorSymbolic(F,Csymm,perm,&info);
221: for (nfact = 0; nfact < 2; nfact++){
222: if (!rank) printf(" Cholesky nfact %d\n",nfact);
223: MatCholeskyFactorNumeric(F,Csymm,&info);
225: /* Test MatSolve() */
226: for (nsolve = 0; nsolve < 2; nsolve++){
227: VecGetArray(x,&array);
228: for (i=0; i<m; i++){
229: PetscRandomGetValue(rand,&rval);
230: array[i] = rval;
231: }
232: VecRestoreArray(x,&array);
233: VecCopy(x,u);
234: MatMult(Csymm,x,b);
235: MatSolve(F,b,x);
237: /* Check the error */
238: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
239: VecNorm(u,NORM_2,&norm);
240: if (norm > tol){
241: if (!rank){
242: PetscPrintf(PETSC_COMM_SELF,"Error: Norm of error %g, Cholesky nfact %d\n",norm,nfact);
243: }
244: }
245: }
246: }
247: MatDestroy(&F);
249: /* Free data structures */
250: ISDestroy(&perm);
251: ISDestroy(&iperm);
252:
253: PetscRandomDestroy(&rand);
254: VecDestroy(&x);
255: VecDestroy(&b);
256: VecDestroy(&bpla);
257: VecDestroy(&u);
258:
259: MatDestroy(&Cpetsc);
260: MatDestroy(&C);
261: MatDestroy(&Csymm);
263: PetscFinalize();
264: return 0;
265: }