Actual source code: ex6.c
petsc-3.14.0 2020-09-29
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;
15: KSPCreateVecs(ksp,2,&t,2,&v);
16: KSPBuildResidualDefault(ksp,t[0],v[0],&v[0]);
17: KSPBuildResidual(ksp,t[1],v[1],&v[1]);
18: VecAXPY(v[1],-1.0,v[0]);
19: VecNorm(v[1],NORM_INFINITY,&err);
20: if (err > PETSC_SMALL) SETERRQ3(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Inconsistent residual computed at step %D: %g (KSP %g)\n",i,(double)err,(double)r);
21: VecDestroyVecs(2,&t);
22: VecDestroyVecs(2,&v);
23: return(0);
24: }
26: int main(int argc,char **args)
27: {
29: PetscInt its;
30: #if defined(PETSC_USE_LOG)
31: PetscLogStage stage1,stage2;
32: #endif
33: PetscReal norm;
34: Vec x,b,u;
35: Mat A;
36: char file[PETSC_MAX_PATH_LEN];
37: PetscViewer fd;
38: PetscBool table = PETSC_FALSE,flg,test_residual = PETSC_FALSE,b_in_f = PETSC_TRUE;
39: KSP ksp;
41: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
42: PetscOptionsGetBool(NULL,NULL,"-table",&table,NULL);
43: PetscOptionsGetBool(NULL,NULL,"-test_residual",&test_residual,NULL);
44: PetscOptionsGetBool(NULL,NULL,"-b_in_f",&b_in_f,NULL);
46: /* Read matrix and RHS */
47: PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
48: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate binary file with the -f option");
49: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
50: MatCreate(PETSC_COMM_WORLD,&A);
51: MatLoad(A,fd);
52: if (b_in_f) {
53: VecCreate(PETSC_COMM_WORLD,&b);
54: VecLoad(b,fd);
55: } else {
56: MatCreateVecs(A,NULL,&b);
57: VecSetRandom(b,NULL);
58: }
59: PetscViewerDestroy(&fd);
61: /*
62: If the load matrix is larger then the vector, due to being padded
63: to match the blocksize then create a new padded vector
64: */
65: {
66: PetscInt m,n,j,mvec,start,end,indx;
67: Vec tmp;
68: PetscScalar *bold;
70: MatGetLocalSize(A,&m,&n);
71: VecCreate(PETSC_COMM_WORLD,&tmp);
72: VecSetSizes(tmp,m,PETSC_DECIDE);
73: VecSetFromOptions(tmp);
74: VecGetOwnershipRange(b,&start,&end);
75: VecGetLocalSize(b,&mvec);
76: VecGetArray(b,&bold);
77: for (j=0; j<mvec; j++) {
78: indx = start+j;
79: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
80: }
81: VecRestoreArray(b,&bold);
82: VecDestroy(&b);
83: VecAssemblyBegin(tmp);
84: VecAssemblyEnd(tmp);
85: b = tmp;
86: }
87: VecDuplicate(b,&x);
88: VecDuplicate(b,&u);
90: VecSet(x,0.0);
91: PetscBarrier((PetscObject)A);
93: PetscLogStageRegister("mystage 1",&stage1);
94: PetscLogStagePush(stage1);
95: KSPCreate(PETSC_COMM_WORLD,&ksp);
96: KSPSetOperators(ksp,A,A);
97: KSPSetFromOptions(ksp);
98: if (test_residual) {
99: KSPMonitorSet(ksp,KSPTestResidualMonitor,NULL,NULL);
100: }
101: KSPSetUp(ksp);
102: KSPSetUpOnBlocks(ksp);
103: PetscLogStagePop();
104: PetscBarrier((PetscObject)A);
106: PetscLogStageRegister("mystage 2",&stage2);
107: PetscLogStagePush(stage2);
108: KSPSolve(ksp,b,x);
109: PetscLogStagePop();
111: /* Show result */
112: MatMult(A,x,u);
113: VecAXPY(u,-1.0,b);
114: VecNorm(u,NORM_2,&norm);
115: KSPGetIterationNumber(ksp,&its);
116: /* matrix PC KSP Options its residual */
117: if (table) {
118: char *matrixname,kspinfo[120];
119: PetscViewer viewer;
120: PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,sizeof(kspinfo),&viewer);
121: KSPView(ksp,viewer);
122: PetscStrrchr(file,'/',&matrixname);
123: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %s \n",matrixname,its,norm,kspinfo);
124: PetscViewerDestroy(&viewer);
125: } else {
126: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
127: PetscPrintf(PETSC_COMM_WORLD,"Residual norm = %g\n",(double)norm);
128: }
130: /* Cleanup */
131: KSPDestroy(&ksp);
132: VecDestroy(&x);
133: VecDestroy(&b);
134: VecDestroy(&u);
135: MatDestroy(&A);
136: PetscFinalize();
137: return ierr;
138: }
140: /*TEST
142: test:
143: args: -ksp_type preonly -pc_type lu -options_left no -f ${DATAFILESPATH}/matrices/arco1
144: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
146: test:
147: suffix: 2
148: 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
149: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
151: test:
152: suffix: 7
153: 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
154: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
156: test:
157: requires: double !complex !define(PETSC_USE_64BIT_INDICES)
158: suffix: 3
159: filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
160: 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 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
162: test:
163: requires: double !complex !define(PETSC_USE_64BIT_INDICES)
164: suffix: 3_maxits
165: output_file: output/ex6_maxits.out
166: 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 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
168: testset:
169: requires: double !complex !define(PETSC_USE_64BIT_INDICES)
170: output_file: output/ex6_skip.out
171: 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
172: #SYMMLQ converges in 4 iterations and then generate nans
173: test:
174: suffix: 3_skip
175: args: -ksp_type {{chebyshev cg groppcg pipecg pipecgrr pipelcg pipeprcg cgne nash stcg gltr fcg pipefcg gmres pipefgmres fgmres lgmres dgmres pgmres tcqmr bcgs ibcgs fbcgs fbcgsr bcgsl pipebcgs cgs tfqmr cr pipecr qcg bicg minres lcd gcr cgls richardson}}
176: #PIPEGCR generates nans on linux-knl
177: test:
178: requires: !define(PETSC_USE_AVX512_KERNELS)
179: suffix: 3_skip_pipegcr
180: args: -ksp_type pipegcr
181: test:
182: requires: hpddm
183: suffix: 3_skip_hpddm
184: args: -ksp_type hpddm -ksp_hpddm_type {{cg gmres bgmres bcg bfbcg gcrodr bgcrodr}}
186: test:
187: requires: double !complex !define(PETSC_USE_64BIT_INDICES) hpddm
188: suffix: 3_hpddm
189: output_file: output/ex6_3.out
190: filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
191: 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
193: # test CG shortcut for residual access
194: test:
195: suffix: 4
196: 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
197: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
200: TEST*/