Actual source code: ex6.c
petsc-3.13.4 2020-08-01
2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
3: Input arguments are:\n\
4: -f <input_file> : file to load. For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\n";
6: #include <petscksp.h>
7: #include <petsclog.h>
9: int main(int argc,char **args)
10: {
12: PetscInt its;
13: PetscLogStage stage1,stage2;
14: PetscReal norm;
15: Vec x,b,u;
16: Mat A;
17: char file[PETSC_MAX_PATH_LEN];
18: PetscViewer fd;
19: PetscBool table = PETSC_FALSE,flg;
20: KSP ksp;
22: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
23: PetscOptionsGetBool(NULL,NULL,"-table",&table,NULL);
25: /* Read matrix and RHS */
26: PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
27: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate binary file with the -f option");
28: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
29: MatCreate(PETSC_COMM_WORLD,&A);
30: MatLoad(A,fd);
31: VecCreate(PETSC_COMM_WORLD,&b);
32: VecLoad(b,fd);
33: PetscViewerDestroy(&fd);
35: /*
36: If the load matrix is larger then the vector, due to being padded
37: to match the blocksize then create a new padded vector
38: */
39: {
40: PetscInt m,n,j,mvec,start,end,indx;
41: Vec tmp;
42: PetscScalar *bold;
44: MatGetLocalSize(A,&m,&n);
45: VecCreate(PETSC_COMM_WORLD,&tmp);
46: VecSetSizes(tmp,m,PETSC_DECIDE);
47: VecSetFromOptions(tmp);
48: VecGetOwnershipRange(b,&start,&end);
49: VecGetLocalSize(b,&mvec);
50: VecGetArray(b,&bold);
51: for (j=0; j<mvec; j++) {
52: indx = start+j;
53: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
54: }
55: VecRestoreArray(b,&bold);
56: VecDestroy(&b);
57: VecAssemblyBegin(tmp);
58: VecAssemblyEnd(tmp);
59: b = tmp;
60: }
61: VecDuplicate(b,&x);
62: VecDuplicate(b,&u);
64: VecSet(x,0.0);
65: PetscBarrier((PetscObject)A);
67: PetscLogStageRegister("mystage 1",&stage1);
68: PetscLogStagePush(stage1);
69: KSPCreate(PETSC_COMM_WORLD,&ksp);
70: KSPSetOperators(ksp,A,A);
71: KSPSetFromOptions(ksp);
72: KSPSetUp(ksp);
73: KSPSetUpOnBlocks(ksp);
74: PetscLogStagePop();
75: PetscBarrier((PetscObject)A);
77: PetscLogStageRegister("mystage 2",&stage2);
78: PetscLogStagePush(stage2);
79: KSPSolve(ksp,b,x);
80: PetscLogStagePop();
82: /* Show result */
83: MatMult(A,x,u);
84: VecAXPY(u,-1.0,b);
85: VecNorm(u,NORM_2,&norm);
86: KSPGetIterationNumber(ksp,&its);
87: /* matrix PC KSP Options its residual */
88: if (table) {
89: char *matrixname,kspinfo[120];
90: PetscViewer viewer;
91: PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,sizeof(kspinfo),&viewer);
92: KSPView(ksp,viewer);
93: PetscStrrchr(file,'/',&matrixname);
94: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %s \n",matrixname,its,norm,kspinfo);
95: PetscViewerDestroy(&viewer);
96: } else {
97: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
98: PetscPrintf(PETSC_COMM_WORLD,"Residual norm = %g\n",(double)norm);
99: }
101: /* Cleanup */
102: KSPDestroy(&ksp);
103: VecDestroy(&x);
104: VecDestroy(&b);
105: VecDestroy(&u);
106: MatDestroy(&A);
107: PetscFinalize();
108: return ierr;
109: }
111: /*TEST
113: test:
114: args: -ksp_type preonly -pc_type lu -options_left no -f ${DATAFILESPATH}/matrices/arco1
115: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
117: test:
118: suffix: 2
119: 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
120: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
122: test:
123: suffix: 7
124: 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
125: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
127: TEST*/