1
2 static char help[] = "Krylov methods to solve u'' = f in parallel with periodic boundary conditions,\n\
3 with a singular, inconsistent system.\n\n";
4
5 /*T
6 Concepts: KSP^basic parallel example
7 Concepts: periodic boundary conditions
8 Processors: n
9 T*/
10
11
12
13 /*
14
15 This tests solving singular inconsistent systems with GMRES
16
17 Default: Solves a symmetric system
18 -symmetric false: Solves a non-symmetric system where nullspace(A) != nullspace(A')
19
20 -ksp_pc_side left or right
21
22 See the KSPSolve() for a discussion of when right preconditioning with nullspace(A) != nullspace(A') can fail to produce the
23 norm minimizing solution.
24
25 Note that though this example does solve the system with right preconditioning and nullspace(A) != nullspace(A') it does not produce the
26 norm minimizing solution, that is the computed solution is not orthogonal to the nullspace(A).
27
28
29 Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
30 Include "petscksp.h" so that we can use KSP solvers. Note that this
31 file automatically includes:
32 petscsys.h - base PETSc routines petscvec.h - vectors
33 petscmat.h - matrices
34 petscis.h - index sets petscksp.h - Krylov subspace methods
35 petscviewer.h - viewers petscpc.h - preconditioners
36 petscksp.h - linear solvers
37 */
38
39 #include <petscdm.h>
40 #include <petscdmda.h>
41 #include <petscsnes.h>
42 #include <petsc/private/kspimpl.h>
43
44 PetscBool symmetric = PETSC_TRUE;
45
46 PetscErrorCode FormMatrix(Mat,void*);
47 PetscErrorCode FormRightHandSide(Vec,void*);
48
main(int argc,char ** argv)49 int main(int argc,char **argv)
50 {
51 KSP ksp;
52 Mat J;
53 DM da;
54 Vec x,r; /* vectors */
55 PetscErrorCode ierr;
56 PetscInt M = 10;
57 MatNullSpace constants,nconstants;
58
59 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
60 ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr);
61 ierr = PetscOptionsGetBool(NULL,NULL,"-symmetric",&symmetric,NULL);CHKERRQ(ierr);
62
63 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64 Create linear solver context
65 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66
67 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
68
69 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70 Create vector data structures; set function evaluation routine
71 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72
73 /*
74 Create distributed array (DMDA) to manage parallel grid and vectors
75 */
76 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,M,1,2,NULL,&da);CHKERRQ(ierr);
77 ierr = DMSetFromOptions(da);CHKERRQ(ierr);
78 ierr = DMSetUp(da);CHKERRQ(ierr);
79
80 /*
81 Extract global and local vectors from DMDA; then duplicate for remaining
82 vectors that are the same types
83 */
84 ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);
85 ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
86
87 /*
88 Set function evaluation routine and vector. Whenever the nonlinear
89 solver needs to compute the nonlinear function, it will call this
90 routine.
91 - Note that the final routine argument is the user-defined
92 context that provides application-specific data for the
93 function evaluation routine.
94 */
95 ierr = FormRightHandSide(r,da);CHKERRQ(ierr);
96
97 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98 Create matrix data structure;
99 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100 ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr);
101 ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants);CHKERRQ(ierr);
102 if (symmetric) {
103 ierr = MatSetOption(J,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr);
104 ierr = MatSetOption(J,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
105 } else {
106 Vec n;
107 PetscInt zero = 0;
108 PetscScalar zeros = 0.0;
109 ierr = VecDuplicate(x,&n);CHKERRQ(ierr);
110 /* the nullspace(A') is the constant vector but with a zero in the very first entry; hence nullspace(A') != nullspace(A) */
111 ierr = VecSet(n,1.0);CHKERRQ(ierr);
112 ierr = VecSetValues(n,1,&zero,&zeros,INSERT_VALUES);CHKERRQ(ierr);
113 ierr = VecAssemblyBegin(n);CHKERRQ(ierr);
114 ierr = VecAssemblyEnd(n);CHKERRQ(ierr);
115 ierr = VecNormalize(n,NULL);CHKERRQ(ierr);
116 ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,1,&n,&nconstants);CHKERRQ(ierr);
117 ierr = MatSetTransposeNullSpace(J,nconstants);CHKERRQ(ierr);
118 ierr = MatNullSpaceDestroy(&nconstants);CHKERRQ(ierr);
119 ierr = VecDestroy(&n);CHKERRQ(ierr);
120 }
121 ierr = MatSetNullSpace(J,constants);CHKERRQ(ierr);
122 ierr = FormMatrix(J,da);CHKERRQ(ierr);
123
124 ierr = KSPSetOperators(ksp,J,J);CHKERRQ(ierr);
125
126 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
127 ierr = KSPSolve(ksp,r,x);CHKERRQ(ierr);
128 ierr = KSPSolveTranspose(ksp,r,x);CHKERRQ(ierr);
129
130 ierr = VecDestroy(&x);CHKERRQ(ierr);
131 ierr = VecDestroy(&r);CHKERRQ(ierr);
132 ierr = MatDestroy(&J);CHKERRQ(ierr);
133 ierr = MatNullSpaceDestroy(&constants);CHKERRQ(ierr);
134 ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
135 ierr = DMDestroy(&da);CHKERRQ(ierr);
136 ierr = PetscFinalize();
137 return ierr;
138 }
139
140 /*
141
142 This intentionally includes something in the right hand side that is not in the range of the matrix A.
143 Since MatSetNullSpace() is called and the matrix is symmetric; the Krylov method automatically removes this
144 portion of the right hand side before solving the linear system.
145 */
FormRightHandSide(Vec f,void * ctx)146 PetscErrorCode FormRightHandSide(Vec f,void *ctx)
147 {
148 DM da = (DM) ctx;
149 PetscScalar *ff;
150 PetscErrorCode ierr;
151 PetscInt i,M,xs,xm;
152 PetscReal h;
153
154 PetscFunctionBeginUser;
155 ierr = DMDAVecGetArray(da,f,&ff);CHKERRQ(ierr);
156
157 /*
158 Get local grid boundaries (for 1-dimensional DMDA):
159 xs, xm - starting grid index, width of local grid (no ghost points)
160 */
161 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);
162 ierr = DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
163
164 /*
165 Compute function over locally owned part of the grid
166 Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
167 */
168 h = 1.0/M;
169 for (i=xs; i<xs+xm; i++) ff[i] = - PetscSinReal(2.0*PETSC_PI*i*h) + 1.0;
170
171 /*
172 Restore vectors
173 */
174 ierr = DMDAVecRestoreArray(da,f,&ff);CHKERRQ(ierr);
175 PetscFunctionReturn(0);
176 }
177 /* ------------------------------------------------------------------- */
FormMatrix(Mat jac,void * ctx)178 PetscErrorCode FormMatrix(Mat jac,void *ctx)
179 {
180 PetscScalar A[3];
181 PetscErrorCode ierr;
182 PetscInt i,M,xs,xm;
183 DM da = (DM) ctx;
184 MatStencil row,cols[3];
185 PetscReal h;
186
187 PetscFunctionBeginUser;
188 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);
189
190 /*
191 Get range of locally owned matrix
192 */
193 ierr = DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
194
195 ierr = MatZeroEntries(jac);CHKERRQ(ierr);
196 h = 1.0/M;
197 /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
198 if (symmetric) {
199 for (i=xs; i<xs+xm; i++) {
200 row.i = i;
201 cols[0].i = i - 1;
202 cols[1].i = i;
203 cols[2].i = i + 1;
204 A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
205 ierr = MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);CHKERRQ(ierr);
206 }
207 } else {
208 MatStencil *acols;
209 PetscScalar *avals;
210
211 /* only works for one process */
212 ierr = MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
213 row.i = 0;
214 ierr = PetscMalloc1(M,&acols);CHKERRQ(ierr);
215 ierr = PetscMalloc1(M,&avals);CHKERRQ(ierr);
216 for (i=0; i<M; i++) {
217 acols[i].i = i;
218 avals[i] = (i % 2) ? 1 : -1;
219 }
220 ierr = MatSetValuesStencil(jac,1,&row,M,acols,avals,ADD_VALUES);CHKERRQ(ierr);
221 ierr = PetscFree(acols);CHKERRQ(ierr);
222 ierr = PetscFree(avals);CHKERRQ(ierr);
223 row.i = 1;
224 cols[0].i = - 1;
225 cols[1].i = 1;
226 cols[2].i = 1 + 1;
227 A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
228 ierr = MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);CHKERRQ(ierr);
229 for (i=2; i<xs+xm-1; i++) {
230 row.i = i;
231 cols[0].i = i - 1;
232 cols[1].i = i;
233 cols[2].i = i + 1;
234 A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
235 ierr = MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);CHKERRQ(ierr);
236 }
237 row.i = M - 1 ;
238 cols[0].i = M-2;
239 cols[1].i = M-1;
240 cols[2].i = M+1;
241 A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
242 ierr = MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);CHKERRQ(ierr);
243 }
244 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
245 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
246 PetscFunctionReturn(0);
247 }
248
249
250 /*TEST
251
252 test:
253 suffix: nonsymmetric_left
254 args: -symmetric false -ksp_view -ksp_converged_reason -pc_type jacobi -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side left
255 filter: sed 's/ATOL/RTOL/g'
256 requires: !single
257
258 test:
259 suffix: nonsymmetric_right
260 args: -symmetric false -ksp_view -ksp_converged_reason -pc_type jacobi -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side right
261 filter: sed 's/ATOL/RTOL/g'
262 requires: !single
263
264 test:
265 suffix: symmetric_left
266 args: -ksp_view -ksp_converged_reason -pc_type sor -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side left
267 requires: !single
268
269 test:
270 suffix: symmetric_right
271 args: -ksp_view -ksp_converged_reason -pc_type sor -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side right
272 filter: sed 's/ATOL/RTOL/g'
273 requires: !single
274
275 test:
276 suffix: transpose_asm
277 args: -symmetric false -ksp_monitor -ksp_view -pc_type asm -sub_pc_type lu -sub_pc_factor_zeropivot 1.e-33 -ksp_converged_reason
278 filter: sed 's/ATOL/RTOL/g'
279
280 TEST*/
281