Actual source code: ex23.c
2: static char help[] = "Solves a tridiagonal linear system.\n\n";
4: /*T
5: Concepts: KSP^basic parallel example;
6: Processors: n
7: T*/
9: /*
10: Include "petscksp.h" so that we can use KSP solvers. Note that this file
11: automatically includes:
12: petscsys.h - base PETSc routines petscvec.h - vectors
13: petscmat.h - matrices
14: petscis.h - index sets petscksp.h - Krylov subspace methods
15: petscviewer.h - viewers petscpc.h - preconditioners
17: Note: The corresponding uniprocessor example is ex1.c
18: */
19: #include <petscksp.h>
21: int main(int argc,char **args)
22: {
23: Vec x, b, u; /* approx solution, RHS, exact solution */
24: Mat A; /* linear system matrix */
25: KSP ksp; /* linear solver context */
26: PC pc; /* preconditioner context */
27: PetscReal norm,tol=1000.*PETSC_MACHINE_EPSILON; /* norm of solution error */
28: PetscInt i,n = 10,col[3],its,rstart,rend,nlocal;
29: PetscScalar one = 1.0,value[3];
31: PetscInitialize(&argc,&args,(char*)0,help);
32: PetscOptionsGetInt(NULL,NULL,"-n",&n,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. For this simple case let PETSc decide how
42: many elements of the vector are stored on each processor. The second
43: argument to VecSetSizes() below causes PETSc to decide.
44: */
45: VecCreate(PETSC_COMM_WORLD,&x);
46: VecSetSizes(x,PETSC_DECIDE,n);
47: VecSetFromOptions(x);
48: VecDuplicate(x,&b);
49: VecDuplicate(x,&u);
51: /* Identify the starting and ending mesh points on each
52: processor for the interior part of the mesh. We let PETSc decide
53: above. */
55: VecGetOwnershipRange(x,&rstart,&rend);
56: VecGetLocalSize(x,&nlocal);
58: /*
59: Create matrix. When using MatCreate(), the matrix format can
60: be specified at runtime.
62: Performance tuning note: For problems of substantial size,
63: preallocation of matrix memory is crucial for attaining good
64: performance. See the matrix chapter of the users manual for details.
66: We pass in nlocal as the "local" size of the matrix to force it
67: to have the same parallel layout as the vector created above.
68: */
69: MatCreate(PETSC_COMM_WORLD,&A);
70: MatSetSizes(A,nlocal,nlocal,n,n);
71: MatSetFromOptions(A);
72: MatSetUp(A);
74: /*
75: Assemble matrix.
77: The linear system is distributed across the processors by
78: chunks of contiguous rows, which correspond to contiguous
79: sections of the mesh on which the problem is discretized.
80: For matrix assembly, each processor contributes entries for
81: the part that it owns locally.
82: */
84: if (!rstart) {
85: rstart = 1;
86: i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
87: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
88: }
89: if (rend == n) {
90: rend = n-1;
91: i = n-1; col[0] = n-2; col[1] = n-1; value[0] = -1.0; value[1] = 2.0;
92: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
93: }
95: /* Set entries corresponding to the mesh interior */
96: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
97: for (i=rstart; i<rend; i++) {
98: col[0] = i-1; col[1] = i; col[2] = i+1;
99: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
100: }
102: /* Assemble the matrix */
103: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
104: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
106: /*
107: Set exact solution; then compute right-hand-side vector.
108: */
109: VecSet(u,one);
110: MatMult(A,u,b);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Create the linear solver and set various options
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: /*
116: Create linear solver context
117: */
118: KSPCreate(PETSC_COMM_WORLD,&ksp);
120: /*
121: Set operators. Here the matrix that defines the linear system
122: also serves as the preconditioning matrix.
123: */
124: KSPSetOperators(ksp,A,A);
126: /*
127: Set linear solver defaults for this problem (optional).
128: - By extracting the KSP and PC contexts from the KSP context,
129: we can then directly call any KSP and PC routines to set
130: various options.
131: - The following four statements are optional; all of these
132: parameters could alternatively be specified at runtime via
133: KSPSetFromOptions();
134: */
135: KSPGetPC(ksp,&pc);
136: PCSetType(pc,PCJACOBI);
137: KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
139: /*
140: Set runtime options, e.g.,
141: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
142: These options will override those specified above as long as
143: KSPSetFromOptions() is called _after_ any other customization
144: routines.
145: */
146: KSPSetFromOptions(ksp);
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Solve the linear system
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: /*
152: Solve linear system
153: */
154: KSPSolve(ksp,b,x);
156: /*
157: View solver info; we could instead use the option -ksp_view to
158: print this info to the screen at the conclusion of KSPSolve().
159: */
160: KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
162: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: Check solution and clean up
164: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: /*
166: Check the error
167: */
168: VecAXPY(x,-1.0,u);
169: VecNorm(x,NORM_2,&norm);
170: KSPGetIterationNumber(ksp,&its);
171: if (norm > tol) {
172: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
173: }
175: /*
176: Free work space. All PETSc objects should be destroyed when they
177: are no longer needed.
178: */
179: VecDestroy(&x)); PetscCall(VecDestroy(&u);
180: VecDestroy(&b)); PetscCall(MatDestroy(&A);
181: KSPDestroy(&ksp);
183: /*
184: Always call PetscFinalize() before exiting a program. This routine
185: - finalizes the PETSc libraries as well as MPI
186: - provides summary and diagnostic information if certain runtime
187: options are chosen (e.g., -log_view).
188: */
189: PetscFinalize();
190: return 0;
191: }
193: /*TEST
195: build:
196: requires: !complex !single
198: test:
199: args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
201: test:
202: suffix: 2
203: nsize: 3
204: args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
206: test:
207: suffix: 3
208: nsize: 2
209: args: -ksp_monitor_short -ksp_rtol 1e-6 -ksp_type pipefgmres
211: TEST*/