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 */
19: PetscInitialize(&argc, &args, (char *)0, help);
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, sizeof(filein), 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 = %3" PetscInt_FMT "\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 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*/