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: PetscMalloc2(m,PetscInt,&im,M,PetscInt,&in);
68: k = 0;
69: for (j=0; j<M; j++){ /* column oriented! */
70: in[j] = j;
71: for (i=0; i<m; i++){
72: im[i] = i+Istart;
73: PetscRandomGetValue(rand,&rval);
74: array[k++] = rval;
75: }
76: }
77: MatSetValues(Cpetsc,m,im,M,in,array,ADD_VALUES);
78: MatSetValues(C,m,im,M,in,array,ADD_VALUES);
79: PetscFree(array);
80: PetscFree2(im,in);
82: MatAssemblyBegin(Cpetsc,MAT_FINAL_ASSEMBLY);
83: MatAssemblyEnd(Cpetsc,MAT_FINAL_ASSEMBLY);
84: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
85: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
86: /*
87: if (!rank) {printf("main, Cpetsc: \n");}
88: MatView(Cpetsc,PETSC_VIEWER_STDOUT_WORLD);
89: */
90: MatGetOrdering(C,MATORDERING_NATURAL,&perm,&iperm);
92: /* Test nonsymmetric MatMult() */
93: VecGetArray(x,&array);
94: for (i=0; i<m; i++){
95: PetscRandomGetValue(rand,&rval);
96: array[i] = rval;
97: }
98: VecRestoreArray(x,&array);
99:
100: MatMult(Cpetsc,x,b);
101: MatMult(C,x,bpla);
102: VecAXPY(bpla,-1.0,b);
103: VecNorm(bpla,NORM_2,&norm);
104: if (norm > 1.e-12 && !rank){
105: PetscPrintf(PETSC_COMM_SELF,"Nonsymmetric MatMult_Plapack error: |b_pla - b|= %g\n",norm);
106: }
108: /* Test LU Factorization */
109: if (nproc == 1){
110: MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&F);
111: } else {
112: MatGetFactor(C,MAT_SOLVER_PLAPACK,MAT_FACTOR_LU,&F);
113: }
114: MatLUFactorSymbolic(F,C,perm,iperm,&info);
115: for (nfact = 0; nfact < 2; nfact++){
116: if (!rank) printf(" LU nfact %d\n",nfact);
117: if (nfact>0){ /* change matrix value for testing repeated MatLUFactorNumeric() */
118: if (!rank){
119: i = j = 0;
120: rval = nfact;
121: MatSetValues(Cpetsc,1,&i,1,&j,&rval,ADD_VALUES);
122: MatSetValues(C,1,&i,1,&j,&rval,ADD_VALUES);
123: } else { /* PLAPACK seems requiring all processors call MatSetValues(), so we add 0.0 on processesses with rank>0! */
124: i = j = 0;
125: rval = 0.0;
126: MatSetValues(Cpetsc,1,&i,1,&j,&rval,ADD_VALUES);
127: MatSetValues(C,1,&i,1,&j,&rval,ADD_VALUES);
128: }
129: MatAssemblyBegin(Cpetsc,MAT_FINAL_ASSEMBLY);
130: MatAssemblyEnd(Cpetsc,MAT_FINAL_ASSEMBLY);
131: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
132: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
133: }
134: MatLUFactorNumeric(F,C,&info);
136: /* Test MatSolve() */
137: for (nsolve = 0; nsolve < 2; nsolve++){
138: VecGetArray(x,&array);
139: for (i=0; i<m; i++){
140: PetscRandomGetValue(rand,&rval);
141: array[i] = rval;
142: /* array[i] = rank + 1; */
143: }
144: VecRestoreArray(x,&array);
145: VecCopy(x,u);
146: MatMult(C,x,b);
147: MatSolve(F,b,x);
149: /* Check the error */
150: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
151: VecNorm(u,NORM_2,&norm);
152: if (norm > tol){
153: if (!rank){
154: PetscPrintf(PETSC_COMM_SELF,"Error: Norm of error %g, LU nfact %d\n",norm,nfact);
155: }
156: }
157: }
158: }
159: MatDestroy(F);
160:
161: /* Test non-symmetric operations */
162: /*-------------------------------*/
163: /* Create a symmetric Plapack dense matrix Csymm */
164: MatCreate(PETSC_COMM_WORLD,&Csymm);
165: MatSetSizes(Csymm,PETSC_DECIDE,PETSC_DECIDE,M,M);
166: MatSetType(Csymm,MATDENSE);
167: MatSetFromOptions(Csymm);
168: MatSetOption(Csymm,MAT_ROW_ORIENTED,PETSC_FALSE);
169: MatSetOption(Csymm,MAT_SYMMETRIC,PETSC_TRUE);
170: MatSetOption(Csymm,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
172: MatZeroEntries(Csymm);
173: MatZeroEntries(Cpetsc);
174: for (i=Istart; i<Iend; i++){
175: for (j=0; j<=i; j++){
176: PetscRandomGetValue(rand,&rval);
177: MatSetValues(Cpetsc,1,&i,1,&j,&rval,ADD_VALUES);
178: MatSetValues(Csymm,1,&i,1,&j,&rval,ADD_VALUES);
179: if (j<i){
180: /* Although PLAPACK only requires lower triangular entries, we must add all the entries.
181: MatSetValues_Plapack() will ignore the upper triangular entries AFTER an index map! */
182: MatSetValues(Cpetsc,1,&j,1,&i,&rval,ADD_VALUES);
183: MatSetValues(Csymm,1,&j,1,&i,&rval,ADD_VALUES);
184: }
185: }
186: }
187: MatAssemblyBegin(Cpetsc,MAT_FINAL_ASSEMBLY);
188: MatAssemblyEnd(Cpetsc,MAT_FINAL_ASSEMBLY);
189: MatAssemblyBegin(Csymm,MAT_FINAL_ASSEMBLY);
190: MatAssemblyEnd(Csymm,MAT_FINAL_ASSEMBLY);
192: /* Test symmetric MatMult() */
193: VecGetArray(x,&array);
194: for (i=0; i<m; i++){
195: PetscRandomGetValue(rand,&rval);
196: array[i] = rval;
197: }
198: VecRestoreArray(x,&array);
199:
200: MatMult(Cpetsc,x,b);
201: MatMult(Csymm,x,bpla);
202: VecAXPY(bpla,-1.0,b);
203: VecNorm(bpla,NORM_2,&norm);
204: if (norm > 1.e-12 && !rank){
205: PetscPrintf(PETSC_COMM_SELF,"Symmetric MatMult_Plapack error: |b_pla - b|= %g\n",norm);
206: }
208: /* Test Cholesky Factorization */
209: MatShift(Csymm,M); /* make Csymm positive definite */
210: if (nproc == 1){
211: MatGetFactor(Csymm,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&F);
212: } else {
213: MatGetFactor(Csymm,MAT_SOLVER_PLAPACK,MAT_FACTOR_CHOLESKY,&F);
214: }
215: MatCholeskyFactorSymbolic(F,Csymm,perm,&info);
216: for (nfact = 0; nfact < 2; nfact++){
217: if (!rank) printf(" Cholesky nfact %d\n",nfact);
218: MatCholeskyFactorNumeric(F,Csymm,&info);
220: /* Test MatSolve() */
221: for (nsolve = 0; nsolve < 2; nsolve++){
222: VecGetArray(x,&array);
223: for (i=0; i<m; i++){
224: PetscRandomGetValue(rand,&rval);
225: array[i] = rval;
226: }
227: VecRestoreArray(x,&array);
228: VecCopy(x,u);
229: MatMult(Csymm,x,b);
230: MatSolve(F,b,x);
232: /* Check the error */
233: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
234: VecNorm(u,NORM_2,&norm);
235: if (norm > tol){
236: if (!rank){
237: PetscPrintf(PETSC_COMM_SELF,"Error: Norm of error %g, Cholesky nfact %d\n",norm,nfact);
238: }
239: }
240: }
241: }
242: MatDestroy(F);
244: /* Free data structures */
245: ISDestroy(perm);
246: ISDestroy(iperm);
247:
248: PetscRandomDestroy(rand);
249: VecDestroy(x);
250: VecDestroy(b);
251: VecDestroy(bpla);
252: VecDestroy(u);
253:
254: MatDestroy(Cpetsc);
255: MatDestroy(C);
256: MatDestroy(Csymm);
258: PetscFinalize();
259: return 0;
260: }