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