Actual source code: ex25.c
petsc-3.12.2 2019-11-22
1: static char help[] = "Tests CG, MINRES and SYMMLQ on the symmetric indefinite matrices: afiro \n\n";
3: #include <petscksp.h>
5: int main(int argc,char **args)
6: {
7: Mat C;
8: PetscScalar none = -1.0;
9: PetscMPIInt rank,size;
11: PetscInt its,k;
12: PetscReal err_norm,res_norm;
13: Vec x,b,u,u_tmp;
14: PC pc;
15: KSP ksp;
16: PetscViewer view;
17: char filein[128]; /* input file name */
19: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
20: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
21: MPI_Comm_size(PETSC_COMM_WORLD,&size);
23: /* Load the binary data file "filein". Set runtime option: -f filein */
24: PetscPrintf(PETSC_COMM_WORLD,"\n Load dataset ...\n");
25: PetscOptionsGetString(NULL,NULL,"-f",filein,128,NULL);
26: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filein,FILE_MODE_READ,&view);
27: MatCreate(PETSC_COMM_WORLD,&C);
28: MatSetType(C,MATMPISBAIJ);
29: MatLoad(C,view);
30: VecCreate(PETSC_COMM_WORLD,&b);
31: VecCreate(PETSC_COMM_WORLD,&u);
32: VecLoad(b,view);
33: VecLoad(u,view);
34: PetscViewerDestroy(&view);
35: /* VecView(b,VIEWER_STDOUT_WORLD); */
36: /* MatView(C,VIEWER_STDOUT_WORLD); */
38: VecDuplicate(u,&u_tmp);
40: /* Check accuracy of the data */
41: /*
42: MatMult(C,u,u_tmp);
43: VecAXPY(u_tmp,none,b);
44: VecNorm(u_tmp,NORM_2,&res_norm);
45: PetscPrintf(PETSC_COMM_WORLD,"Accuracy of the loading data: | b - A*u |_2 : %g \n",(double)res_norm);
46: */
48: /* Setup and solve for system */
49: VecDuplicate(b,&x);
50: for (k=0; k<3; k++) {
51: if (k == 0) { /* CG */
52: KSPCreate(PETSC_COMM_WORLD,&ksp);
53: KSPSetType(ksp,KSPCG);
54: KSPSetOperators(ksp,C,C);
55: PetscPrintf(PETSC_COMM_WORLD,"\n CG: \n");
56: } else if (k == 1) { /* MINRES */
57: KSPCreate(PETSC_COMM_WORLD,&ksp);
58: KSPSetType(ksp,KSPMINRES);
59: KSPSetOperators(ksp,C,C);
60: PetscPrintf(PETSC_COMM_WORLD,"\n MINRES: \n");
61: } else { /* SYMMLQ */
62: KSPCreate(PETSC_COMM_WORLD,&ksp);
63: KSPSetOperators(ksp,C,C);
64: KSPSetType(ksp,KSPSYMMLQ);
65: PetscPrintf(PETSC_COMM_WORLD,"\n SYMMLQ: \n");
66: }
68: KSPGetPC(ksp,&pc);
69: PCSetType(pc,PCNONE);
71: /*
72: Set runtime options, e.g.,
73: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
74: -pc_type jacobi -pc_jacobi_type rowmax
75: These options will override those specified above as long as
76: KSPSetFromOptions() is called _after_ any other customization routines.
77: */
78: KSPSetFromOptions(ksp);
80: /* Solve linear system; */
81: KSPSolve(ksp,b,x);
82: KSPGetIterationNumber(ksp,&its);
84: /* Check error */
85: VecCopy(u,u_tmp);
86: VecAXPY(u_tmp,none,x);
87: VecNorm(u_tmp,NORM_2,&err_norm);
88: MatMult(C,x,u_tmp);
89: VecAXPY(u_tmp,none,b);
90: VecNorm(u_tmp,NORM_2,&res_norm);
92: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
93: PetscPrintf(PETSC_COMM_WORLD,"Residual norm: %g;",(double)res_norm);
94: PetscPrintf(PETSC_COMM_WORLD," Error norm: %g.\n",(double)err_norm);
96: KSPDestroy(&ksp);
97: }
99: /*
100: Free work space. All PETSc objects should be destroyed when they
101: are no longer needed.
102: */
103: VecDestroy(&b);
104: VecDestroy(&u);
105: VecDestroy(&x);
106: VecDestroy(&u_tmp);
107: MatDestroy(&C);
109: PetscFinalize();
110: return ierr;
111: }
113: /*TEST
115: test:
116: args: -f ${DATAFILESPATH}/matrices/indefinite/afiro -ksp_rtol 1.e-3
117: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
119: TEST*/