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: PetscInitialize(&argc,&args,(char*)0,help);
19: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
20: MPI_Comm_size(PETSC_COMM_WORLD,&size);
22: /* Load the binary data file "filein". Set runtime option: -f filein */
23: PetscPrintf(PETSC_COMM_WORLD,"\n Load dataset ...\n");
24: PetscOptionsGetString(NULL,NULL,"-f",filein,sizeof(filein),NULL);
25: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filein,FILE_MODE_READ,&view);
26: MatCreate(PETSC_COMM_WORLD,&C);
27: MatSetType(C,MATMPISBAIJ);
28: MatLoad(C,view);
29: VecCreate(PETSC_COMM_WORLD,&b);
30: VecCreate(PETSC_COMM_WORLD,&u);
31: VecLoad(b,view);
32: VecLoad(u,view);
33: PetscViewerDestroy(&view);
34: /* VecView(b,VIEWER_STDOUT_WORLD); */
35: /* MatView(C,VIEWER_STDOUT_WORLD); */
37: VecDuplicate(u,&u_tmp);
39: /* Check accuracy of the data */
40: /*
41: MatMult(C,u,u_tmp);
42: VecAXPY(u_tmp,none,b);
43: VecNorm(u_tmp,NORM_2,&res_norm);
44: PetscPrintf(PETSC_COMM_WORLD,"Accuracy of the loading data: | b - A*u |_2 : %g \n",(double)res_norm);
45: */
47: /* Setup and solve for system */
48: VecDuplicate(b,&x);
49: for (k=0; k<3; k++) {
50: if (k == 0) { /* CG */
51: KSPCreate(PETSC_COMM_WORLD,&ksp);
52: KSPSetType(ksp,KSPCG);
53: KSPSetOperators(ksp,C,C);
54: PetscPrintf(PETSC_COMM_WORLD,"\n CG: \n");
55: } else if (k == 1) { /* MINRES */
56: KSPCreate(PETSC_COMM_WORLD,&ksp);
57: KSPSetType(ksp,KSPMINRES);
58: KSPSetOperators(ksp,C,C);
59: PetscPrintf(PETSC_COMM_WORLD,"\n MINRES: \n");
60: } else { /* SYMMLQ */
61: KSPCreate(PETSC_COMM_WORLD,&ksp);
62: KSPSetOperators(ksp,C,C);
63: KSPSetType(ksp,KSPSYMMLQ);
64: PetscPrintf(PETSC_COMM_WORLD,"\n SYMMLQ: \n");
65: }
67: KSPGetPC(ksp,&pc);
68: PCSetType(pc,PCNONE);
70: /*
71: Set runtime options, e.g.,
72: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
73: -pc_type jacobi -pc_jacobi_type rowmax
74: These options will override those specified above as long as
75: KSPSetFromOptions() is called _after_ any other customization routines.
76: */
77: KSPSetFromOptions(ksp);
79: /* Solve linear system; */
80: KSPSolve(ksp,b,x);
81: KSPGetIterationNumber(ksp,&its);
83: /* Check error */
84: VecCopy(u,u_tmp);
85: VecAXPY(u_tmp,none,x);
86: VecNorm(u_tmp,NORM_2,&err_norm);
87: MatMult(C,x,u_tmp);
88: VecAXPY(u_tmp,none,b);
89: VecNorm(u_tmp,NORM_2,&res_norm);
91: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
92: PetscPrintf(PETSC_COMM_WORLD,"Residual norm: %g;",(double)res_norm);
93: PetscPrintf(PETSC_COMM_WORLD," Error norm: %g.\n",(double)err_norm);
95: KSPDestroy(&ksp);
96: }
98: /*
99: Free work space. All PETSc objects should be destroyed when they
100: are no longer needed.
101: */
102: VecDestroy(&b);
103: VecDestroy(&u);
104: VecDestroy(&x);
105: VecDestroy(&u_tmp);
106: MatDestroy(&C);
108: PetscFinalize();
109: return 0;
110: }
112: /*TEST
114: test:
115: args: -f ${DATAFILESPATH}/matrices/indefinite/afiro -ksp_rtol 1.e-3
116: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
118: TEST*/