Actual source code: ex25.c
petsc-3.3-p6 2013-02-11
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,128,PETSC_NULL);
32: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filein,FILE_MODE_READ,&view);
33: MatCreate(PETSC_COMM_WORLD,&C);
34: MatSetType(C,MATMPISBAIJ);
35: MatLoad(C,view);
36: VecCreate(PETSC_COMM_WORLD,&b);
37: VecCreate(PETSC_COMM_WORLD,&u);
38: VecLoad(b,view);
39: VecLoad(u,view);
40: PetscViewerDestroy(&view);
41: /* VecView(b,VIEWER_STDOUT_WORLD); */
42: /* MatView(C,VIEWER_STDOUT_WORLD); */
44: VecDuplicate(u,&u_tmp);
46: /* Check accuracy of the data */
47: /*
48: MatMult(C,u,u_tmp);
49: VecAXPY(u_tmp,none,b);
50: VecNorm(u_tmp,NORM_2,&res_norm);
51: PetscPrintf(PETSC_COMM_WORLD,"Accuracy of the loading data: | b - A*u |_2 : %G \n",res_norm);
52: */
54: /* Setup and solve for system */
55: VecDuplicate(b,&x);
56: for (k=0; k<3; k++){
57: if (k == 0){ /* CG */
58: KSPCreate(PETSC_COMM_WORLD,&ksp);
59: KSPSetType(ksp,KSPCG);
60: KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);
61: PetscPrintf(PETSC_COMM_WORLD,"\n CG: \n");
62: } else if (k == 1){ /* MINRES */
63: KSPCreate(PETSC_COMM_WORLD,&ksp);
64: KSPSetType(ksp,KSPMINRES);
65: KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);
66: PetscPrintf(PETSC_COMM_WORLD,"\n MINRES: \n");
67: } else { /* SYMMLQ */
68: KSPCreate(PETSC_COMM_WORLD,&ksp);
69: KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);
70: KSPSetType(ksp,KSPSYMMLQ);
71: PetscPrintf(PETSC_COMM_WORLD,"\n SYMMLQ: \n");
72: }
74: KSPGetPC(ksp,&pc);
75: PCSetType(pc,PCNONE);
76: /* PCSetType(pc,PCJACOBI); */
77: KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
79: /*
80: Set runtime options, e.g.,
81: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
82: -pc_type jacobi -pc_jacobi_rowmax
83: These options will override those specified above as long as
84: KSPSetFromOptions() is called _after_ any other customization routines.
85: */
86: KSPSetFromOptions(ksp);
88: /* Solve linear system; */
89: KSPSolve(ksp,b,x);
90: KSPGetIterationNumber(ksp,&its);
91:
92: /* Check error */
93: VecCopy(u,u_tmp);
94: VecAXPY(u_tmp,none,x);
95: VecNorm(u_tmp,NORM_2,&err_norm);
96: MatMult(C,x,u_tmp);
97: VecAXPY(u_tmp,none,b);
98: VecNorm(u_tmp,NORM_2,&res_norm);
99:
100: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
101: PetscPrintf(PETSC_COMM_WORLD,"Residual norm: %G;",res_norm);
102: PetscPrintf(PETSC_COMM_WORLD," Error norm: %G.\n",err_norm);
104: KSPDestroy(&ksp);
105: }
106:
107: /*
108: Free work space. All PETSc objects should be destroyed when they
109: are no longer needed.
110: */
111: VecDestroy(&b);
112: VecDestroy(&u);
113: VecDestroy(&x);
114: VecDestroy(&u_tmp);
115: MatDestroy(&C);
117: PetscFinalize();
118: return 0;
119: }