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