Actual source code: ex6.c

  1: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
  2: Input arguments are:\n\
  3:   -f <input_file> : file to load. For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\n";

  5: #include <petscksp.h>
  6: #include <petsclog.h>

  8: static PetscErrorCode KSPTestResidualMonitor(KSP ksp, PetscInt i, PetscReal r, void *ctx)
  9: {
 10:   Vec      *t, *v;
 11:   PetscReal err;

 14:   KSPCreateVecs(ksp, 2, &t, 2, &v);
 15:   KSPBuildResidualDefault(ksp, t[0], v[0], &v[0]);
 16:   KSPBuildResidual(ksp, t[1], v[1], &v[1]);
 17:   VecAXPY(v[1], -1.0, v[0]);
 18:   VecNorm(v[1], NORM_INFINITY, &err);
 20:   VecDestroyVecs(2, &t);
 21:   VecDestroyVecs(2, &v);
 22:   return 0;
 23: }

 25: int main(int argc, char **args)
 26: {
 27:   PetscInt its;
 28: #if defined(PETSC_USE_LOG)
 29:   PetscLogStage stage1, stage2;
 30: #endif
 31:   PetscReal   norm;
 32:   Vec         x, b, u;
 33:   Mat         A;
 34:   char        file[PETSC_MAX_PATH_LEN];
 35:   PetscViewer fd;
 36:   PetscBool   table = PETSC_FALSE, flg, test_residual = PETSC_FALSE, b_in_f = PETSC_TRUE;
 37:   KSP         ksp;

 40:   PetscInitialize(&argc, &args, (char *)0, help);
 41:   PetscOptionsGetBool(NULL, NULL, "-table", &table, NULL);
 42:   PetscOptionsGetBool(NULL, NULL, "-test_residual", &test_residual, NULL);
 43:   PetscOptionsGetBool(NULL, NULL, "-b_in_f", &b_in_f, NULL);

 45:   /* Read matrix and RHS */
 46:   PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg);
 48:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);
 49:   MatCreate(PETSC_COMM_WORLD, &A);
 50:   MatLoad(A, fd);
 51:   if (b_in_f) {
 52:     VecCreate(PETSC_COMM_WORLD, &b);
 53:     VecLoad(b, fd);
 54:   } else {
 55:     MatCreateVecs(A, NULL, &b);
 56:     VecSetRandom(b, NULL);
 57:   }
 58:   PetscViewerDestroy(&fd);

 60:   /*
 61:    If the load matrix is larger then the vector, due to being padded
 62:    to match the blocksize then create a new padded vector
 63:   */
 64:   {
 65:     PetscInt     m, n, j, mvec, start, end, indx;
 66:     Vec          tmp;
 67:     PetscScalar *bold;

 69:     MatGetLocalSize(A, &m, &n);
 70:     VecCreate(PETSC_COMM_WORLD, &tmp);
 71:     VecSetSizes(tmp, m, PETSC_DECIDE);
 72:     VecSetFromOptions(tmp);
 73:     VecGetOwnershipRange(b, &start, &end);
 74:     VecGetLocalSize(b, &mvec);
 75:     VecGetArray(b, &bold);
 76:     for (j = 0; j < mvec; j++) {
 77:       indx = start + j;
 78:       VecSetValues(tmp, 1, &indx, bold + j, INSERT_VALUES);
 79:     }
 80:     VecRestoreArray(b, &bold);
 81:     VecDestroy(&b);
 82:     VecAssemblyBegin(tmp);
 83:     VecAssemblyEnd(tmp);
 84:     b = tmp;
 85:   }
 86:   VecDuplicate(b, &x);
 87:   VecDuplicate(b, &u);

 89:   VecSet(x, 0.0);
 90:   PetscBarrier((PetscObject)A);

 92:   PetscLogStageRegister("mystage 1", &stage1);
 93:   PetscLogStagePush(stage1);
 94:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 95:   KSPSetOperators(ksp, A, A);
 96:   KSPSetFromOptions(ksp);
 97:   if (test_residual) KSPMonitorSet(ksp, KSPTestResidualMonitor, NULL, NULL);
 98:   KSPSetUp(ksp);
 99:   KSPSetUpOnBlocks(ksp);
100:   PetscLogStagePop();
101:   PetscBarrier((PetscObject)A);

103:   PetscLogStageRegister("mystage 2", &stage2);
104:   PetscLogStagePush(stage2);
105:   KSPSolve(ksp, b, x);
106:   PetscLogStagePop();

108:   /* Show result */
109:   MatMult(A, x, u);
110:   VecAXPY(u, -1.0, b);
111:   VecNorm(u, NORM_2, &norm);
112:   KSPGetIterationNumber(ksp, &its);
113:   /*  matrix PC   KSP   Options       its    residual  */
114:   if (table) {
115:     char       *matrixname, kspinfo[120];
116:     PetscViewer viewer;
117:     PetscViewerStringOpen(PETSC_COMM_WORLD, kspinfo, sizeof(kspinfo), &viewer);
118:     KSPView(ksp, viewer);
119:     PetscStrrchr(file, '/', &matrixname);
120:     PetscPrintf(PETSC_COMM_WORLD, "%-8.8s %3" PetscInt_FMT " %2.0e %s \n", matrixname, its, (double)norm, kspinfo);
121:     PetscViewerDestroy(&viewer);
122:   } else {
123:     PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its);
124:     PetscPrintf(PETSC_COMM_WORLD, "Residual norm = %g\n", (double)norm);
125:   }

127:   /* Cleanup */
128:   KSPDestroy(&ksp);
129:   VecDestroy(&x);
130:   VecDestroy(&b);
131:   VecDestroy(&u);
132:   MatDestroy(&A);
133:   PetscFinalize();
134:   return 0;
135: }

137: /*TEST

139:     test:
140:       args: -ksp_type preonly  -pc_type lu -options_left no  -f ${DATAFILESPATH}/matrices/arco1
141:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)

143:     test:
144:       suffix: 2
145:       args: -sub_pc_type ilu -options_left no  -f ${DATAFILESPATH}/matrices/arco1 -ksp_gmres_restart 100 -ksp_gmres_cgs_refinement_type refine_always -sub_ksp_type preonly -pc_type bjacobi -pc_bjacobi_blocks 8 -sub_pc_factor_in_place -ksp_monitor_short
146:       requires: datafilespath double  !complex !defined(PETSC_USE_64BIT_INDICES)

148:     test:
149:       suffix: 7
150:       args: -ksp_gmres_cgs_refinement_type refine_always -pc_type asm -pc_asm_blocks 6 -f ${DATAFILESPATH}/matrices/small -matload_block_size 6  -ksp_monitor_short
151:       requires: datafilespath double  !complex !defined(PETSC_USE_64BIT_INDICES)

153:     test:
154:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
155:       suffix: 3
156:       filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
157:       args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64 -pc_type none -ksp_type {{cg groppcg pipecg pipecgrr pipelcg pipeprcg cgne nash stcg gltr fcg pipefcg gmres pipefgmres fgmres lgmres dgmres pgmres tcqmr bcgs ibcgs qmrcgs fbcgs fbcgsr bcgsl pipebcgs cgs tfqmr cr pipecr lsqr qcg bicg minres symmlq lcd gcr pipegcr cgls}} -ksp_max_it 20 -ksp_error_if_not_converged -ksp_converged_reason -test_residual

159:     test:
160:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
161:       suffix: 3_maxits
162:       output_file: output/ex6_maxits.out
163:       args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64 -pc_type none -ksp_type {{chebyshev cg groppcg pipecg pipecgrr pipelcg pipeprcg cgne nash stcg gltr fcg pipefcg gmres pipefgmres fgmres lgmres dgmres pgmres tcqmr bcgs ibcgs qmrcgs fbcgs fbcgsr bcgsl pipebcgs cgs tfqmr cr pipecr qcg bicg minres symmlq lcd gcr pipegcr cgls richardson}} -ksp_max_it 4 -ksp_error_if_not_converged -ksp_converged_maxits -ksp_converged_reason -test_residual -ksp_norm_type none

165:     testset:
166:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
167:       output_file: output/ex6_skip.out
168:       args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64 -pc_type none -ksp_max_it 8 -ksp_error_if_not_converged -ksp_convergence_test skip -ksp_converged_reason -test_residual
169:       #SYMMLQ converges in 4 iterations and then generate nans
170:       test:
171:         suffix: 3_skip
172:         args: -ksp_type {{chebyshev cg groppcg pipecg pipecgrr pipelcg pipeprcg cgne nash stcg gltr fcg pipefcg gmres fgmres lgmres dgmres pgmres tcqmr bcgs ibcgs qmrcgs fbcgs fbcgsr bcgsl pipebcgs cgs tfqmr cr pipecr qcg bicg minres lcd gcr cgls richardson}}
173:       test:
174:         requires: !pgf90_compiler
175:         suffix: 3_skip_pipefgmres
176:         args: -ksp_type pipefgmres
177:       #PIPEGCR generates nans on linux-knl
178:       test:
179:         requires: !defined(PETSC_USE_AVX512_KERNELS)
180:         suffix: 3_skip_pipegcr
181:         args: -ksp_type pipegcr
182:       test:
183:         requires: hpddm
184:         suffix: 3_skip_hpddm
185:         args: -ksp_type hpddm -ksp_hpddm_type {{cg gmres bgmres bcg bfbcg gcrodr bgcrodr}}

187:     test:
188:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES) hpddm
189:       suffix: 3_hpddm
190:       output_file: output/ex6_3.out
191:       filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
192:       args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64 -pc_type none -ksp_type hpddm -ksp_hpddm_type {{cg gmres bgmres bcg bfbcg gcrodr bgcrodr}} -ksp_max_it 20 -ksp_error_if_not_converged -ksp_converged_reason -test_residual

194:     # test CG shortcut for residual access
195:     test:
196:       suffix: 4
197:       args: -ksp_converged_reason -ksp_max_it 20 -ksp_converged_maxits -ksp_type {{cg pipecg groppcg}} -ksp_norm_type {{preconditioned unpreconditioned natural}separate output} -pc_type {{bjacobi none}separate output} -f ${DATAFILESPATH}/matrices/poisson_2d13p -b_in_f 0 -test_residual
198:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)

200: TEST*/