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: }