Actual source code: ex25.c
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;
10: PetscInt its, k;
11: PetscReal err_norm, res_norm;
12: Vec x, b, u, u_tmp;
13: PC pc;
14: KSP ksp;
15: PetscViewer view;
16: char filein[128]; /* input file name */
18: PetscFunctionBeginUser;
19: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
20: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
21: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
23: /* Load the binary data file "filein". Set runtime option: -f filein */
24: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n Load dataset ...\n"));
25: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filein, sizeof(filein), NULL));
26: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filein, FILE_MODE_READ, &view));
27: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
28: PetscCall(MatSetType(C, MATMPISBAIJ));
29: PetscCall(MatLoad(C, view));
30: PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
31: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
32: PetscCall(VecLoad(b, view));
33: PetscCall(VecLoad(u, view));
34: PetscCall(PetscViewerDestroy(&view));
35: /* PetscCall(VecView(b,VIEWER_STDOUT_WORLD)); */
36: /* PetscCall(MatView(C,VIEWER_STDOUT_WORLD)); */
38: PetscCall(VecDuplicate(u, &u_tmp));
40: /* Check accuracy of the data */
41: /*
42: PetscCall(MatMult(C,u,u_tmp));
43: PetscCall(VecAXPY(u_tmp,none,b));
44: PetscCall(VecNorm(u_tmp,NORM_2,&res_norm));
45: PetscCall(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: PetscCall(VecDuplicate(b, &x));
50: for (k = 0; k < 3; k++) {
51: if (k == 0) { /* CG */
52: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
53: PetscCall(KSPSetType(ksp, KSPCG));
54: PetscCall(KSPSetOperators(ksp, C, C));
55: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n CG: \n"));
56: } else if (k == 1) { /* MINRES */
57: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
58: PetscCall(KSPSetType(ksp, KSPMINRES));
59: PetscCall(KSPSetOperators(ksp, C, C));
60: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MINRES: \n"));
61: } else { /* SYMMLQ */
62: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
63: PetscCall(KSPSetOperators(ksp, C, C));
64: PetscCall(KSPSetType(ksp, KSPSYMMLQ));
65: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n SYMMLQ: \n"));
66: }
68: PetscCall(KSPGetPC(ksp, &pc));
69: PetscCall(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: PetscCall(KSPSetFromOptions(ksp));
80: /* Solve linear system; */
81: PetscCall(KSPSolve(ksp, b, x));
82: PetscCall(KSPGetIterationNumber(ksp, &its));
84: /* Check error */
85: PetscCall(VecCopy(u, u_tmp));
86: PetscCall(VecAXPY(u_tmp, none, x));
87: PetscCall(VecNorm(u_tmp, NORM_2, &err_norm));
88: PetscCall(MatMult(C, x, u_tmp));
89: PetscCall(VecAXPY(u_tmp, none, b));
90: PetscCall(VecNorm(u_tmp, NORM_2, &res_norm));
92: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its));
93: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm: %g;", (double)res_norm));
94: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm: %g.\n", (double)err_norm));
96: PetscCall(KSPDestroy(&ksp));
97: }
99: /*
100: Free work space. All PETSc objects should be destroyed when they
101: are no longer needed.
102: */
103: PetscCall(VecDestroy(&b));
104: PetscCall(VecDestroy(&u));
105: PetscCall(VecDestroy(&x));
106: PetscCall(VecDestroy(&u_tmp));
107: PetscCall(MatDestroy(&C));
109: PetscCall(PetscFinalize());
110: return 0;
111: }
113: /*TEST
115: test:
116: args: -f ${DATAFILESPATH}/matrices/indefinite/afiro -ksp_rtol 1.e-3
117: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
119: TEST*/