2: /* Usage: mpiexec ex16 [-help] [all PETSc options] */
4: static char help[] = "Solves a sequence of linear systems with different right-hand-side vectors.\n\
5: Input parameters include:\n\
6: -ntimes <ntimes> : number of linear systems to solve\n\
7: -view_exact_sol : write exact solution vector to stdout\n\
8: -m <mesh_x> : number of mesh points in x-direction\n\
9: -n <mesh_n> : number of mesh points in y-direction\n\n";
11: /*T
12: Concepts: KSP^repeatedly solving linear systems;
13: Concepts: KSP^Laplacian, 2d
14: Concepts: Laplacian, 2d
15: Processors: n
16: T*/
18: /*
19: Include "petscksp.h" so that we can use KSP solvers. Note that this file
20: automatically includes:
21: petscsys.h - base PETSc routines petscvec.h - vectors
22: petscmat.h - matrices
23: petscis.h - index sets petscksp.h - Krylov subspace methods
24: petscviewer.h - viewers petscpc.h - preconditioners
25: */
26: #include <petscksp.h>
30: int main(int argc,char **args) 31: {
32: Vec x,b,u; /* approx solution, RHS, exact solution */
33: Mat A; /* linear system matrix */
34: KSP ksp; /* linear solver context */
35: PetscReal norm; /* norm of solution error */
37: PetscInt ntimes,i,j,k,Ii,J,Istart,Iend;
38: PetscInt m = 8,n = 7,its;
39: PetscBool flg = PETSC_FALSE;
40: PetscScalar v,one = 1.0,neg_one = -1.0,rhs;
42: PetscInitialize(&argc,&args,(char*)0,help);
43: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
44: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
46: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: Compute the matrix for use in solving a series of
48: linear systems of the form, A x_i = b_i, for i=1,2,...
49: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50: /*
51: Create parallel matrix, specifying only its global dimensions.
52: When using MatCreate(), the matrix format can be specified at
53: runtime. Also, the parallel partitioning of the matrix is
54: determined by PETSc at runtime.
55: */
56: MatCreate(PETSC_COMM_WORLD,&A);
57: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
58: MatSetFromOptions(A);
59: MatSetUp(A);
61: /*
62: Currently, all PETSc parallel matrix formats are partitioned by
63: contiguous chunks of rows across the processors. Determine which
64: rows of the matrix are locally owned.
65: */
66: MatGetOwnershipRange(A,&Istart,&Iend);
68: /*
69: Set matrix elements for the 2-D, five-point stencil in parallel.
70: - Each processor needs to insert only elements that it owns
71: locally (but any non-local elements will be sent to the
72: appropriate processor during matrix assembly).
73: - Always specify global rows and columns of matrix entries.
74: */
75: for (Ii=Istart; Ii<Iend; Ii++) {
76: v = -1.0; i = Ii/n; j = Ii - i*n;
77: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
78: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
79: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
80: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
81: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
82: }
84: /*
85: Assemble matrix, using the 2-step process:
86: MatAssemblyBegin(), MatAssemblyEnd()
87: Computations can be done while messages are in transition
88: by placing code between these two statements.
89: */
90: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
91: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
93: /*
94: Create parallel vectors.
95: - When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
96: we specify only the vector's global
97: dimension; the parallel partitioning is determined at runtime.
98: - When solving a linear system, the vectors and matrices MUST
99: be partitioned accordingly. PETSc automatically generates
100: appropriately partitioned matrices and vectors when MatCreate()
101: and VecCreate() are used with the same communicator.
102: - Note: We form 1 vector from scratch and then duplicate as needed.
103: */
104: VecCreate(PETSC_COMM_WORLD,&u);
105: VecSetSizes(u,PETSC_DECIDE,m*n);
106: VecSetFromOptions(u);
107: VecDuplicate(u,&b);
108: VecDuplicate(b,&x);
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Create the linear solver and set various options
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: /*
115: Create linear solver context
116: */
117: KSPCreate(PETSC_COMM_WORLD,&ksp);
119: /*
120: Set operators. Here the matrix that defines the linear system
121: also serves as the preconditioning matrix.
122: */
123: KSPSetOperators(ksp,A,A);
125: /*
126: Set runtime options, e.g.,
127: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
128: These options will override those specified above as long as
129: KSPSetFromOptions() is called _after_ any other customization
130: routines.
131: */
132: KSPSetFromOptions(ksp);
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Solve several linear systems of the form A x_i = b_i
136: I.e., we retain the same matrix (A) for all systems, but
137: change the right-hand-side vector (b_i) at each step.
139: In this case, we simply call KSPSolve() multiple times. The
140: preconditioner setup operations (e.g., factorization for ILU)
141: be done during the first call to KSPSolve() only; such operations
142: will NOT be repeated for successive solves.
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: ntimes = 2;
146: PetscOptionsGetInt(NULL,NULL,"-ntimes",&ntimes,NULL);
147: for (k=1; k<ntimes+1; k++) {
149: /*
150: Set exact solution; then compute right-hand-side vector. We use
151: an exact solution of a vector with all elements equal to 1.0*k.
152: */
153: rhs = one * (PetscReal)k;
154: VecSet(u,rhs);
155: MatMult(A,u,b);
157: /*
158: View the exact solution vector if desired
159: */
160: PetscOptionsGetBool(NULL,NULL,"-view_exact_sol",&flg,NULL);
161: if (flg) {VecView(u,PETSC_VIEWER_STDOUT_WORLD);}
163: KSPSolve(ksp,b,x);
165: /*
166: Check the error
167: */
168: VecAXPY(x,neg_one,u);
169: VecNorm(x,NORM_2,&norm);
170: KSPGetIterationNumber(ksp,&its);
171: /*
172: Print convergence information. PetscPrintf() produces a single
173: print statement from all processes that share a communicator.
174: */
175: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g System %D: iterations %D\n",(double)norm,k,its);
176: }
178: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: Clean up
180: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
181: /*
182: Free work space. All PETSc objects should be destroyed when they
183: are no longer needed.
184: */
185: KSPDestroy(&ksp);
186: VecDestroy(&u); VecDestroy(&x);
187: VecDestroy(&b); MatDestroy(&A);
189: /*
190: Always call PetscFinalize() before exiting a program. This routine
191: - finalizes the PETSc libraries as well as MPI
192: - provides summary and diagnostic information if certain runtime
193: options are chosen (e.g., -log_summary).
194: */
195: PetscFinalize();
196: return 0;
197: }