Actual source code: ex11.c
petsc-3.8.4 2018-03-24
2: static char help[] = "Solves a linear system in parallel with KSP.\n\n";
4: /*T
5: Concepts: KSP^solving a Helmholtz equation
6: Concepts: complex numbers;
7: Concepts: Helmholtz equation
8: Processors: n
9: T*/
11: /*
12: Description: Solves a complex linear system in parallel with KSP.
14: The model problem:
15: Solve Helmholtz equation on the unit square: (0,1) x (0,1)
16: -delta u - sigma1*u + i*sigma2*u = f,
17: where delta = Laplace operator
18: Dirichlet b.c.'s on all sides
19: Use the 2-D, five-point finite difference stencil.
21: Compiling the code:
22: This code uses the complex numbers version of PETSc, so configure
23: must be run to enable this
24: */
26: /*
27: Include "petscksp.h" so that we can use KSP solvers. Note that this file
28: automatically includes:
29: petscsys.h - base PETSc routines petscvec.h - vectors
30: petscmat.h - matrices
31: petscis.h - index sets petscksp.h - Krylov subspace methods
32: petscviewer.h - viewers petscpc.h - preconditioners
33: */
34: #include <petscksp.h>
36: int main(int argc,char **args)
37: {
38: Vec x,b,u; /* approx solution, RHS, exact solution */
39: Mat A; /* linear system matrix */
40: KSP ksp; /* linear solver context */
41: PetscReal norm; /* norm of solution error */
42: PetscInt dim,i,j,Ii,J,Istart,Iend,n = 6,its,use_random;
44: PetscScalar v,none = -1.0,sigma2,pfive = 0.5,*xa;
45: PetscRandom rctx;
46: PetscReal h2,sigma1 = 100.0;
47: PetscBool flg = PETSC_FALSE;
49: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
50: PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);
51: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
52: dim = n*n;
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Compute the matrix and right-hand-side vector that define
56: the linear system, Ax = b.
57: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58: /*
59: Create parallel matrix, specifying only its global dimensions.
60: When using MatCreate(), the matrix format can be specified at
61: runtime. Also, the parallel partitioning of the matrix is
62: determined by PETSc at runtime.
63: */
64: MatCreate(PETSC_COMM_WORLD,&A);
65: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
66: MatSetFromOptions(A);
67: MatSetUp(A);
69: /*
70: Currently, all PETSc parallel matrix formats are partitioned by
71: contiguous chunks of rows across the processors. Determine which
72: rows of the matrix are locally owned.
73: */
74: MatGetOwnershipRange(A,&Istart,&Iend);
76: /*
77: Set matrix elements in parallel.
78: - Each processor needs to insert only elements that it owns
79: locally (but any non-local elements will be sent to the
80: appropriate processor during matrix assembly).
81: - Always specify global rows and columns of matrix entries.
82: */
84: PetscOptionsGetBool(NULL,NULL,"-norandom",&flg,NULL);
85: if (flg) use_random = 0;
86: else use_random = 1;
87: if (use_random) {
88: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
89: PetscRandomSetFromOptions(rctx);
90: PetscRandomSetInterval(rctx,0.0,PETSC_i);
91: } else {
92: sigma2 = 10.0*PETSC_i;
93: }
94: h2 = 1.0/((n+1)*(n+1));
95: for (Ii=Istart; Ii<Iend; Ii++) {
96: v = -1.0; i = Ii/n; j = Ii - i*n;
97: if (i>0) {
98: J = Ii-n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
99: }
100: if (i<n-1) {
101: J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
102: }
103: if (j>0) {
104: J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
105: }
106: if (j<n-1) {
107: J = Ii+1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
108: }
109: if (use_random) {PetscRandomGetValue(rctx,&sigma2);}
110: v = 4.0 - sigma1*h2 + sigma2*h2;
111: MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
112: }
113: if (use_random) {PetscRandomDestroy(&rctx);}
115: /*
116: Assemble matrix, using the 2-step process:
117: MatAssemblyBegin(), MatAssemblyEnd()
118: Computations can be done while messages are in transition
119: by placing code between these two statements.
120: */
121: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
122: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
124: /*
125: Create parallel vectors.
126: - When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
127: we specify only the vector's global
128: dimension; the parallel partitioning is determined at runtime.
129: - Note: We form 1 vector from scratch and then duplicate as needed.
130: */
131: VecCreate(PETSC_COMM_WORLD,&u);
132: VecSetSizes(u,PETSC_DECIDE,dim);
133: VecSetFromOptions(u);
134: VecDuplicate(u,&b);
135: VecDuplicate(b,&x);
137: /*
138: Set exact solution; then compute right-hand-side vector.
139: */
141: if (use_random) {
142: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
143: PetscRandomSetFromOptions(rctx);
144: VecSetRandom(u,rctx);
145: } else {
146: VecSet(u,pfive);
147: }
148: MatMult(A,u,b);
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Create the linear solver and set various options
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: /*
155: Create linear solver context
156: */
157: KSPCreate(PETSC_COMM_WORLD,&ksp);
159: /*
160: Set operators. Here the matrix that defines the linear system
161: also serves as the preconditioning matrix.
162: */
163: KSPSetOperators(ksp,A,A);
165: /*
166: Set runtime options, e.g.,
167: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
168: */
169: KSPSetFromOptions(ksp);
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Solve the linear system
173: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: KSPSolve(ksp,b,x);
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Check solution and clean up
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
181: /*
182: Print the first 3 entries of x; this demonstrates extraction of the
183: real and imaginary components of the complex vector, x.
184: */
185: flg = PETSC_FALSE;
186: PetscOptionsGetBool(NULL,NULL,"-print_x3",&flg,NULL);
187: if (flg) {
188: VecGetArray(x,&xa);
189: PetscPrintf(PETSC_COMM_WORLD,"The first three entries of x are:\n");
190: for (i=0; i<3; i++) {
191: PetscPrintf(PETSC_COMM_WORLD,"x[%D] = %g + %g i\n",i,(double)PetscRealPart(xa[i]),(double)PetscImaginaryPart(xa[i]));
192: }
193: VecRestoreArray(x,&xa);
194: }
196: /*
197: Check the error
198: */
199: VecAXPY(x,none,u);
200: VecNorm(x,NORM_2,&norm);
201: KSPGetIterationNumber(ksp,&its);
202: if (norm < 1.e-12) {
203: PetscPrintf(PETSC_COMM_WORLD,"Norm of error < 1.e-12 iterations %D\n",its);
204: } else {
205: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);
206: }
208: /*
209: Free work space. All PETSc objects should be destroyed when they
210: are no longer needed.
211: */
212: KSPDestroy(&ksp);
213: if (use_random) {PetscRandomDestroy(&rctx);}
214: VecDestroy(&u); VecDestroy(&x);
215: VecDestroy(&b); MatDestroy(&A);
216: PetscFinalize();
217: return ierr;
218: }