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