Actual source code: ex4.c


  2: static char help[] = "Solves a linear system in parallel with KSP and HMG.\n\
  3: Input parameters include:\n\
  4:   -view_exact_sol    : write exact solution vector to stdout\n\
  5:   -m  <mesh_x>       : number of mesh points in x-direction\n\
  6:   -n  <mesh_y>       : number of mesh points in y-direction\n\
  7:   -bs                : number of variables on each mesh vertex \n\n";

  9: /*
 10:   Simple example is used to test PCHMG
 11: */
 12: #include <petscksp.h>

 14: int main(int argc,char **args)
 15: {
 16:   Vec            x,b,u;    /* approx solution, RHS, exact solution */
 17:   Mat            A;        /* linear system matrix */
 18:   KSP            ksp;      /* linear solver context */
 19:   PetscReal      norm;     /* norm of solution error */
 20:   PetscInt       i,j,Ii,J,Istart,Iend,m = 8,n = 7,its,bs=1,II,JJ,jj;
 21:   PetscBool      flg,test=PETSC_FALSE,reuse=PETSC_FALSE,viewexpl=PETSC_FALSE;
 22:   PetscScalar    v;
 23:   PC             pc;

 25:   PetscInitialize(&argc,&args,(char*)0,help);
 26:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 27:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 28:   PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
 29:   PetscOptionsGetBool(NULL,NULL,"-test_hmg_interface",&test,NULL);
 30:   PetscOptionsGetBool(NULL,NULL,"-test_reuse_interpolation",&reuse,NULL);
 31:   PetscOptionsGetBool(NULL,NULL,"-view_explicit_mat",&viewexpl,NULL);

 33:   MatCreate(PETSC_COMM_WORLD,&A);
 34:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n*bs,m*n*bs);
 35:   MatSetBlockSize(A,bs);
 36:   MatSetFromOptions(A);
 37:   MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
 38:   MatSeqAIJSetPreallocation(A,5,NULL);
 39: #if defined(PETSC_HAVE_HYPRE)
 40:   MatHYPRESetPreallocation(A,5,NULL,5,NULL);
 41: #endif

 43:   MatGetOwnershipRange(A,&Istart,&Iend);

 45:   for (Ii=Istart/bs; Ii<Iend/bs; Ii++) {
 46:     v = -1.0; i = Ii/n; j = Ii - i*n;
 47:     if (i>0) {
 48:       J = Ii - n;
 49:       for (jj=0; jj<bs; jj++) {
 50:         II = Ii*bs + jj;
 51:         JJ = J*bs + jj;
 52:         MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
 53:       }
 54:     }
 55:     if (i<m-1) {
 56:       J = Ii + n;
 57:       for (jj=0; jj<bs; jj++) {
 58:         II = Ii*bs + jj;
 59:         JJ = J*bs + jj;
 60:         MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
 61:       }
 62:     }
 63:     if (j>0) {
 64:       J = Ii - 1;
 65:       for (jj=0; jj<bs; jj++) {
 66:         II = Ii*bs + jj;
 67:         JJ = J*bs + jj;
 68:         MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
 69:       }
 70:     }
 71:     if (j<n-1) {
 72:       J = Ii + 1;
 73:       for (jj=0; jj<bs; jj++) {
 74:         II = Ii*bs + jj;
 75:         JJ = J*bs + jj;
 76:         MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
 77:       }
 78:     }
 79:     v = 4.0;
 80:     for (jj=0; jj<bs; jj++) {
 81:       II = Ii*bs + jj;
 82:       MatSetValues(A,1,&II,1,&II,&v,ADD_VALUES);
 83:     }
 84:   }

 86:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 87:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 88:   if (viewexpl) {
 89:     Mat E;
 90:     MatComputeOperator(A,MATAIJ,&E);
 91:     MatView(E,NULL);
 92:     MatDestroy(&E);
 93:   }

 95:   MatCreateVecs(A,&u,NULL);
 96:   VecSetFromOptions(u);
 97:   VecDuplicate(u,&b);
 98:   VecDuplicate(b,&x);

100:   VecSet(u,1.0);
101:   MatMult(A,u,b);

103:   flg  = PETSC_FALSE;
104:   PetscOptionsGetBool(NULL,NULL,"-view_exact_sol",&flg,NULL);
105:   if (flg) VecView(u,PETSC_VIEWER_STDOUT_WORLD);

107:   KSPCreate(PETSC_COMM_WORLD,&ksp);
108:   KSPSetOperators(ksp,A,A);
109:   KSPSetTolerances(ksp,1.e-2/((m+1)*(n+1)),1.e-50,PETSC_DEFAULT,PETSC_DEFAULT);

111:   if (test) {
112:     KSPGetPC(ksp,&pc);
113:     PCSetType(pc,PCHMG);
114:     PCHMGSetInnerPCType(pc,PCGAMG);
115:     PCHMGSetReuseInterpolation(pc,PETSC_TRUE);
116:     PCHMGSetUseSubspaceCoarsening(pc,PETSC_TRUE);
117:     PCHMGUseMatMAIJ(pc,PETSC_FALSE);
118:     PCHMGSetCoarseningComponent(pc,0);
119:   }

121:   KSPSetFromOptions(ksp);
122:   KSPSolve(ksp,b,x);

124:   if (reuse) {
125:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
126:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
127:     KSPSolve(ksp,b,x);
128:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
129:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
130:     KSPSolve(ksp,b,x);
131:     /* Make sparsity pattern different and reuse interpolation */
132:     MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
133:     MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);
134:     MatGetSize(A,&m,NULL);
135:     n = 0;
136:     v = 0;
137:     m--;
138:     /* Connect the last element to the first element */
139:     MatSetValue(A,m,n,v,ADD_VALUES);
140:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
141:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
142:     KSPSolve(ksp,b,x);
143:   }

145:   VecAXPY(x,-1.0,u);
146:   VecNorm(x,NORM_2,&norm);
147:   KSPGetIterationNumber(ksp,&its);

149:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);

151:   KSPDestroy(&ksp);
152:   VecDestroy(&u);
153:   VecDestroy(&x);
154:   VecDestroy(&b);
155:   MatDestroy(&A);

157:   PetscFinalize();
158:   return 0;
159: }

161: /*TEST

163:    build:
164:       requires: !complex !single

166:    test:
167:       suffix: hypre
168:       nsize: 2
169:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
170:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre

172:    test:
173:       suffix: hypre_bs4
174:       nsize: 2
175:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
176:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -pc_hmg_use_subspace_coarsening 1

178:    test:
179:       suffix: hypre_asm
180:       nsize: 2
181:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
182:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -pc_hmg_use_subspace_coarsening 1 -mg_levels_3_pc_type asm

184:    test:
185:       suffix: hypre_fieldsplit
186:       nsize: 2
187:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
188:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -mg_levels_4_pc_type fieldsplit

190:    test:
191:       suffix: gamg
192:       nsize: 2
193:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg

195:    test:
196:       suffix: gamg_bs4
197:       nsize: 2
198:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -pc_hmg_use_subspace_coarsening 1

200:    test:
201:       suffix: gamg_asm
202:       nsize: 2
203:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -pc_hmg_use_subspace_coarsening 1 -mg_levels_1_pc_type asm

205:    test:
206:       suffix: gamg_fieldsplit
207:       nsize: 2
208:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -mg_levels_1_pc_type fieldsplit

210:    test:
211:       suffix: interface
212:       nsize: 2
213:       args: -ksp_monitor -ksp_rtol 1e-6 -test_hmg_interface 1 -bs 4

215:    test:
216:       suffix: reuse
217:       nsize: 2
218:       args: -ksp_monitor -ksp_rtol 1e-6   -pc_type hmg -pc_hmg_reuse_interpolation 1 -test_reuse_interpolation 1 -hmg_inner_pc_type gamg

220:    test:
221:       suffix: component
222:       nsize: 2
223:       args: -ksp_monitor -ksp_rtol 1e-6 -pc_type hmg -pc_hmg_coarsening_component 2  -pc_hmg_use_subspace_coarsening 1 -bs 4 -hmg_inner_pc_type gamg

225:    testset:
226:       output_file: output/ex4_expl.out
227:       nsize: {{1 2}}
228:       filter: grep -v "MPI processes" | grep -v " type:" | grep -v "Mat Object"
229:       args: -ksp_converged_reason -view_explicit_mat -pc_type none -ksp_type {{cg gmres}}
230:       test:
231:         suffix: expl_aij
232:         args: -mat_type aij
233:       test:
234:         suffix: expl_hypre
235:         requires: hypre
236:         args: -mat_type hypre

238:    test:
239:       suffix: hypre_device
240:       nsize: {{1 2}}
241:       requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
242:       args: -mat_type hypre -ksp_converged_reason -pc_type hypre -m 13 -n 17

244:    test:
245:       suffix: hypre_device_cusparse
246:       output_file: output/ex4_hypre_device.out
247:       nsize: {{1 2}}
248:       requires: hypre cuda defined(PETSC_HAVE_HYPRE_DEVICE)
249:       args: -mat_type {{aij aijcusparse}} -vec_type cuda -ksp_converged_reason -pc_type hypre -m 13 -n 17

251: TEST*/