Actual source code: ex58.c
2: static char help[] = "Solves a tridiagonal linear system with KSP.\n\n";
4: /*T
5: Concepts: KSP^solving a system of linear equations
6: Processors: 1
7: T*/
9: /*
10: Modified from ex1.c for testing matrix operations when matrix structure is changed.
11: Contributed by Jose E. Roman, Feb. 2012.
12: */
13: #include <petscksp.h>
15: int main(int argc,char **args)
16: {
17: Vec x, b, u; /* approx solution, RHS, exact solution */
18: Mat A,B,C; /* linear system matrix */
19: KSP ksp; /* linear solver context */
20: PC pc; /* preconditioner context */
21: PetscReal norm; /* norm of solution error */
23: PetscInt i,n = 20,col[3],its;
24: PetscMPIInt size;
25: PetscScalar one = 1.0,value[3];
26: PetscBool nonzeroguess = PETSC_FALSE;
28: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
29: MPI_Comm_size(PETSC_COMM_WORLD,&size);
30: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
31: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
32: PetscOptionsGetBool(NULL,NULL,"-nonzero_guess",&nonzeroguess,NULL);
34: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
35: Compute the matrix and right-hand-side vector that define
36: the linear system, Ax = b.
37: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
39: /*
40: Create vectors. Note that we form 1 vector from scratch and
41: then duplicate as needed.
42: */
43: VecCreate(PETSC_COMM_WORLD,&x);
44: PetscObjectSetName((PetscObject) x, "Solution");
45: VecSetSizes(x,PETSC_DECIDE,n);
46: VecSetFromOptions(x);
47: VecDuplicate(x,&b);
48: VecDuplicate(x,&u);
50: /*
51: Create matrix. When using MatCreate(), the matrix format can
52: be specified at runtime.
54: Performance tuning note: For problems of substantial size,
55: preallocation of matrix memory is crucial for attaining good
56: performance. See the matrix chapter of the users manual for details.
57: */
58: MatCreate(PETSC_COMM_WORLD,&A);
59: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
60: MatSetFromOptions(A);
61: MatSetUp(A);
63: /*
64: Assemble matrix
65: */
66: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
67: for (i=1; i<n-1; i++) {
68: col[0] = i-1; col[1] = i; col[2] = i+1;
69: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
70: }
71: i = n - 1; col[0] = n - 2; col[1] = n - 1;
72: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
73: i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
74: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
75: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
76: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
78: /*
79: jroman: added matrices
80: */
81: MatCreate(PETSC_COMM_WORLD,&B);
82: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n);
83: MatSetFromOptions(B);
84: MatSetUp(B);
85: for (i=0; i<n; i++) {
86: MatSetValue(B,i,i,value[1],INSERT_VALUES);
87: if (n-i+n/3<n) {
88: MatSetValue(B,n-i+n/3,i,value[0],INSERT_VALUES);
89: MatSetValue(B,i,n-i+n/3,value[0],INSERT_VALUES);
90: }
91: }
92: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
93: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
94: MatDuplicate(A,MAT_COPY_VALUES,&C);
95: MatAXPY(C,2.0,B,DIFFERENT_NONZERO_PATTERN);
97: /*
98: Set exact solution; then compute right-hand-side vector.
99: */
100: VecSet(u,one);
101: MatMult(C,u,b);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Create the linear solver and set various options
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: /*
107: Create linear solver context
108: */
109: KSPCreate(PETSC_COMM_WORLD,&ksp);
111: /*
112: Set operators. Here the matrix that defines the linear system
113: also serves as the preconditioning matrix.
114: */
115: KSPSetOperators(ksp,C,C);
117: /*
118: Set linear solver defaults for this problem (optional).
119: - By extracting the KSP and PC contexts from the KSP context,
120: we can then directly call any KSP and PC routines to set
121: various options.
122: - The following four statements are optional; all of these
123: parameters could alternatively be specified at runtime via
124: KSPSetFromOptions();
125: */
126: KSPGetPC(ksp,&pc);
127: PCSetType(pc,PCJACOBI);
128: KSPSetTolerances(ksp,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
130: /*
131: Set runtime options, e.g.,
132: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
133: These options will override those specified above as long as
134: KSPSetFromOptions() is called _after_ any other customization
135: routines.
136: */
137: KSPSetFromOptions(ksp);
139: if (nonzeroguess) {
140: PetscScalar p = .5;
141: VecSet(x,p);
142: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
143: }
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Solve the linear system
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: /*
149: Solve linear system
150: */
151: KSPSolve(ksp,b,x);
153: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: Check solution and clean up
155: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156: /*
157: Check the error
158: */
159: VecAXPY(x,-1.0,u);
160: VecNorm(x,NORM_2,&norm);
161: KSPGetIterationNumber(ksp,&its);
162: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
164: /*
165: Free work space. All PETSc objects should be destroyed when they
166: are no longer needed.
167: */
168: VecDestroy(&x); VecDestroy(&u);
169: VecDestroy(&b); MatDestroy(&A);
170: MatDestroy(&B);
171: MatDestroy(&C);
172: KSPDestroy(&ksp);
174: /*
175: Always call PetscFinalize() before exiting a program. This routine
176: - finalizes the PETSc libraries as well as MPI
177: - provides summary and diagnostic information if certain runtime
178: options are chosen (e.g., -log_view).
179: */
180: PetscFinalize();
181: return ierr;
182: }
184: /*TEST
186: test:
187: args: -mat_type aij
188: output_file: output/ex58.out
190: test:
191: suffix: baij
192: args: -mat_type baij
193: output_file: output/ex58.out
195: test:
196: suffix: sbaij
197: args: -mat_type sbaij
198: output_file: output/ex58.out
200: TEST*/