Actual source code: ex25.c
1: /*$Id: ex25.c,v 1.3 2000/11/15 22:56:05 balay Exp $*/
3: static char help[] =
4: "Tests CG, MINRES and SYMMLQ on the symmetric indefinite matrices: afiro and golan\n\
5: Runtime options: ex25 -fload ~petsc/matrices/indefinite/afiro -pc_type jacobi -pc_jacobi_rowmax\n\
6: See ~petsc/matrices/indefinite/readme \n\n";
8: #include petscksp.h
12: int main(int argc,char **args)
13: {
14: Mat C;
15: PetscScalar v,none = -1.0;
16: int i,j,ierr,Istart,Iend,N,rank,size,its,k;
17: double err_norm,res_norm;
18: Vec x,b,u,u_tmp;
19: PetscRandom r;
20: PC pc;
21: KSP ksp;
22: PetscViewer view;
23: char filein[128]; /* input file name */
25: PetscInitialize(&argc,&args,(char *)0,help);
26: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
27: MPI_Comm_size(PETSC_COMM_WORLD,&size);
28:
29: /* Load the binary data file "filein". Set runtime option: -fload filein */
30: PetscPrintf(PETSC_COMM_WORLD,"\n Load dataset ...\n");
31: PetscOptionsGetString(PETSC_NULL,"-fload",filein,127,PETSC_NULL);
32: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filein,FILE_MODE_READ,&view);
33: MatLoad(view,MATMPISBAIJ,&C);
34: VecLoad(view,VECMPI,&b);
35: VecLoad(view,VECMPI,&u);
36: PetscViewerDestroy(view);
37: /* VecView(b,VIEWER_STDOUT_WORLD); */
38: /* MatView(C,VIEWER_STDOUT_WORLD); */
40: VecDuplicate(u,&u_tmp);
42: /* Check accuracy of the data */
43: /*
44: MatMult(C,u,u_tmp);
45: VecAXPY(u_tmp,none,b);
46: VecNorm(u_tmp,NORM_2,&res_norm);
47: PetscPrintf(PETSC_COMM_WORLD,"Accuracy of the loading data: | b - A*u |_2 : %A \n",res_norm);
48: */
50: /* Setup and solve for system */
51: VecDuplicate(b,&x);
52: for (k=0; k<3; k++){
53: if (k == 0){ /* CG */
54: KSPCreate(PETSC_COMM_WORLD,&ksp);
55: KSPSetType(ksp,KSPCG);
56: KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);
57: PetscPrintf(PETSC_COMM_WORLD,"\n CG: \n");
58: } else if (k == 1){ /* MINRES */
59: KSPCreate(PETSC_COMM_WORLD,&ksp);
60: KSPSetType(ksp,KSPMINRES);
61: KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);
62: PetscPrintf(PETSC_COMM_WORLD,"\n MINRES: \n");
63: } else { /* SYMMLQ */
64: KSPCreate(PETSC_COMM_WORLD,&ksp);
65: KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);
66: KSPSetType(ksp,KSPSYMMLQ);
67: PetscPrintf(PETSC_COMM_WORLD,"\n SYMMLQ: \n");
68: }
70: KSPGetPC(ksp,&pc);
71: PCSetType(pc,PCNONE);
72: /* PCSetType(pc,PCJACOBI); */
73: KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
75: /*
76: Set runtime options, e.g.,
77: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
78: -pc_type jacobi -pc_jacobi_rowmax
79: These options will override those specified above as long as
80: KSPSetFromOptions() is called _after_ any other customization routines.
81: */
82: KSPSetFromOptions(ksp);
84: /* Solve linear system; */
85: KSPSolve(ksp,b,x);
86: KSPGetIterationNumber(ksp,&its);
87:
88: /* Check error */
89: VecCopy(u,u_tmp);
90: VecAXPY(u_tmp,none,x);
91: VecNorm(u_tmp,NORM_2,&err_norm);
92: MatMult(C,x,u_tmp);
93: VecAXPY(u_tmp,none,b);
94: VecNorm(u_tmp,NORM_2,&res_norm);
95:
96: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
97: PetscPrintf(PETSC_COMM_WORLD,"Residual norm: %A;",res_norm);
98: PetscPrintf(PETSC_COMM_WORLD," Error norm: %A.\n",err_norm);
100: KSPDestroy(ksp);
101: }
102:
103: /*
104: Free work space. All PETSc objects should be destroyed when they
105: are no longer needed.
106: */
107: VecDestroy(b);
108: VecDestroy(u);
109: VecDestroy(x);
110: VecDestroy(u_tmp);
111: MatDestroy(C);
113: PetscFinalize();
114: return 0;
115: }