Actual source code: ex24.c
1: /*$Id: ex24.c,v 1.13 2001/08/07 21:30:50 bsmith Exp $*/
3: static char help[] = "Tests CG, MINRES and SYMMLQ on symmetric matrices with SBAIJ format. The preconditioner ICC only works on sequential SBAIJ format. \n\n";
5: #include petscsles.h
10: int main(int argc,char **args)
11: {
12: Mat C;
13: PetscScalar v,none = -1.0;
14: int i,j,I,J,ierr,Istart,Iend,N,m = 4,n = 4,rank,size,its,k;
15: PetscReal err_norm,res_norm;
16: Vec x,b,u,u_tmp;
17: PetscRandom r;
18: SLES sles;
19: PC pc;
20: KSP ksp;
22: PetscInitialize(&argc,&args,(char *)0,help);
23: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
26: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
27: N = m*n;
30: /* Generate matrix */
31: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&C);
32: MatSetFromOptions(C);
33: MatGetOwnershipRange(C,&Istart,&Iend);
34: for (I=Istart; I<Iend; I++) {
35: v = -1.0; i = I/n; j = I - i*n;
36: if (i>0) {J = I - n; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
37: if (i<m-1) {J = I + n; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
38: if (j>0) {J = I - 1; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
39: if (j<n-1) {J = I + 1; MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);}
40: v = 4.0; MatSetValues(C,1,&I,1,&I,&v,ADD_VALUES);
41: }
42: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
43: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
45: /* a shift can make C indefinite. Preconditioners LU, ILU (for BAIJ format) and ICC may fail */
46: /* MatShift(&alpha, C); */
47: /* MatView(C,PETSC_VIEWER_STDOUT_WORLD); */
49: /* Setup and solve for system */
50:
51: /* Create vectors. */
52: VecCreate(PETSC_COMM_WORLD,&x);
53: VecSetSizes(x,PETSC_DECIDE,N);
54: VecSetFromOptions(x);
55: VecDuplicate(x,&b);
56: VecDuplicate(x,&u);
57: VecDuplicate(x,&u_tmp);
59: /* Set exact solution u; then compute right-hand-side vector b. */
60: PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&r);
61: VecSetRandom(r,u);
62: PetscRandomDestroy(r);
63:
64: MatMult(C,u,b);
66: for (k=0; k<3; k++){
67: if (k == 0){ /* CG */
68: SLESCreate(PETSC_COMM_WORLD,&sles);
69: SLESSetOperators(sles,C,C,DIFFERENT_NONZERO_PATTERN);
70: SLESGetKSP(sles,&ksp);
71: PetscPrintf(PETSC_COMM_WORLD,"\n CG: \n");
72: KSPSetType(ksp,KSPCG);
73: } else if (k == 1){ /* MINRES */
74: SLESCreate(PETSC_COMM_WORLD,&sles);
75: SLESSetOperators(sles,C,C,DIFFERENT_NONZERO_PATTERN);
76: SLESGetKSP(sles,&ksp);
77: PetscPrintf(PETSC_COMM_WORLD,"\n MINRES: \n");
78: KSPSetType(ksp,KSPMINRES);
79: } else { /* SYMMLQ */
80: SLESCreate(PETSC_COMM_WORLD,&sles);
81: SLESSetOperators(sles,C,C,DIFFERENT_NONZERO_PATTERN);
82: SLESGetKSP(sles,&ksp);
83: PetscPrintf(PETSC_COMM_WORLD,"\n SYMMLQ: \n");
84: KSPSetType(ksp,KSPSYMMLQ);
85: }
87: SLESGetPC(sles,&pc);
88: /* PCSetType(pc,PCICC); */
89: PCSetType(pc,PCJACOBI);
90: KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
92: /*
93: Set runtime options, e.g.,
94: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
95: These options will override those specified above as long as
96: SLESSetFromOptions() is called _after_ any other customization
97: routines.
98: */
99: SLESSetFromOptions(sles);
101: /* Solve linear system; */
102: SLESSolve(sles,b,x,&its);
103:
104: /* Check error */
105: VecCopy(u,u_tmp);
106: VecAXPY(&none,x,u_tmp);
107: VecNorm(u_tmp,NORM_2,&err_norm);
108: MatMult(C,x,u_tmp);
109: VecAXPY(&none,b,u_tmp);
110: VecNorm(u_tmp,NORM_2,&res_norm);
111:
112: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
113: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %A;",res_norm);
114: PetscPrintf(PETSC_COMM_WORLD," Error norm %A.\n",err_norm);
116: SLESDestroy(sles);
117: }
118:
119: /*
120: Free work space. All PETSc objects should be destroyed when they
121: are no longer needed.
122: */
123: VecDestroy(b);
124: VecDestroy(u);
125: VecDestroy(x);
126: VecDestroy(u_tmp);
127: MatDestroy(C);
129: PetscFinalize();
130: return 0;
131: }