Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : static char help[] = "Simple 1-D nonlinear eigenproblem.\n\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = number of grid subdivisions.\n"
14 : " -draw_sol, to draw the computed solution.\n\n";
15 :
16 : /*
17 : Solve 1-D PDE
18 : -u'' = lambda*u
19 : on [0,1] subject to
20 : u(0)=0, u'(1)=u(1)*lambda*kappa/(kappa-lambda)
21 : */
22 :
23 : #include <slepcnep.h>
24 :
25 : /*
26 : User-defined routines
27 : */
28 : PetscErrorCode FormInitialGuess(Vec);
29 : PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
30 : PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);
31 : PetscErrorCode CheckSolution(PetscScalar,Vec,PetscReal*,void*);
32 : PetscErrorCode FixSign(Vec);
33 :
34 : /*
35 : User-defined application context
36 : */
37 : typedef struct {
38 : PetscScalar kappa; /* ratio between stiffness of spring and attached mass */
39 : PetscReal h; /* mesh spacing */
40 : } ApplicationCtx;
41 :
42 5 : int main(int argc,char **argv)
43 : {
44 5 : NEP nep; /* nonlinear eigensolver context */
45 5 : Vec x; /* eigenvector */
46 5 : PetscScalar lambda; /* eigenvalue */
47 5 : Mat F,J; /* Function and Jacobian matrices */
48 5 : ApplicationCtx ctx; /* user-defined context */
49 5 : NEPType type;
50 5 : PetscInt n=128,nev,i,its,maxit,nconv;
51 5 : PetscReal re,im,tol,norm,error;
52 5 : PetscBool draw_sol=PETSC_FALSE;
53 :
54 5 : PetscFunctionBeginUser;
55 5 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
56 5 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
57 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%" PetscInt_FMT "\n\n",n));
58 5 : ctx.h = 1.0/(PetscReal)n;
59 5 : ctx.kappa = 1.0;
60 5 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-draw_sol",&draw_sol,NULL));
61 :
62 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63 : Create nonlinear eigensolver context
64 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65 :
66 5 : PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
67 :
68 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69 : Create matrix data structure; set Function evaluation routine
70 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71 :
72 5 : PetscCall(MatCreate(PETSC_COMM_WORLD,&F));
73 5 : PetscCall(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n));
74 5 : PetscCall(MatSetFromOptions(F));
75 5 : PetscCall(MatSeqAIJSetPreallocation(F,3,NULL));
76 5 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
77 :
78 : /*
79 : Set Function matrix data structure and default Function evaluation
80 : routine
81 : */
82 5 : PetscCall(NEPSetFunction(nep,F,F,FormFunction,&ctx));
83 :
84 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85 : Create matrix data structure; set Jacobian evaluation routine
86 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87 :
88 5 : PetscCall(MatCreate(PETSC_COMM_WORLD,&J));
89 5 : PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n));
90 5 : PetscCall(MatSetFromOptions(J));
91 5 : PetscCall(MatSeqAIJSetPreallocation(J,3,NULL));
92 5 : PetscCall(MatMPIAIJSetPreallocation(J,3,NULL,1,NULL));
93 :
94 : /*
95 : Set Jacobian matrix data structure and default Jacobian evaluation
96 : routine
97 : */
98 5 : PetscCall(NEPSetJacobian(nep,J,FormJacobian,&ctx));
99 :
100 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101 : Customize nonlinear solver; set runtime options
102 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103 :
104 5 : PetscCall(NEPSetTolerances(nep,1e-9,PETSC_CURRENT));
105 5 : PetscCall(NEPSetDimensions(nep,1,PETSC_DETERMINE,PETSC_DETERMINE));
106 :
107 : /*
108 : Set solver parameters at runtime
109 : */
110 5 : PetscCall(NEPSetFromOptions(nep));
111 :
112 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113 : Initialize application
114 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115 :
116 : /*
117 : Evaluate initial guess
118 : */
119 5 : PetscCall(MatCreateVecs(F,&x,NULL));
120 5 : PetscCall(FormInitialGuess(x));
121 5 : PetscCall(NEPSetInitialSpace(nep,1,&x));
122 :
123 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124 : Solve the eigensystem
125 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126 :
127 5 : PetscCall(NEPSolve(nep));
128 5 : PetscCall(NEPGetIterationNumber(nep,&its));
129 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of NEP iterations = %" PetscInt_FMT "\n\n",its));
130 :
131 : /*
132 : Optional: Get some information from the solver and display it
133 : */
134 5 : PetscCall(NEPGetType(nep,&type));
135 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type));
136 5 : PetscCall(NEPGetDimensions(nep,&nev,NULL,NULL));
137 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
138 5 : PetscCall(NEPGetTolerances(nep,&tol,&maxit));
139 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit));
140 :
141 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142 : Display solution and clean up
143 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144 :
145 : /*
146 : Get number of converged approximate eigenpairs
147 : */
148 5 : PetscCall(NEPGetConverged(nep,&nconv));
149 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %" PetscInt_FMT "\n\n",nconv));
150 :
151 5 : if (nconv>0) {
152 : /*
153 : Display eigenvalues and relative errors
154 : */
155 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,
156 : " k ||T(k)x|| error\n"
157 : " ----------------- ------------------ ------------------\n"));
158 20 : for (i=0;i<nconv;i++) {
159 : /*
160 : Get converged eigenpairs (in this example they are always real)
161 : */
162 15 : PetscCall(NEPGetEigenpair(nep,i,&lambda,NULL,x,NULL));
163 15 : PetscCall(FixSign(x));
164 : /*
165 : Compute residual norm and error
166 : */
167 15 : PetscCall(NEPComputeError(nep,i,NEP_ERROR_RELATIVE,&norm));
168 15 : PetscCall(CheckSolution(lambda,x,&error,&ctx));
169 :
170 : #if defined(PETSC_USE_COMPLEX)
171 : re = PetscRealPart(lambda);
172 : im = PetscImaginaryPart(lambda);
173 : #else
174 15 : re = lambda;
175 15 : im = 0.0;
176 : #endif
177 15 : if (im!=0.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12g %12g\n",(double)re,(double)im,(double)norm,(double)error));
178 15 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD," %12f %12g %12g\n",(double)re,(double)norm,(double)error));
179 15 : if (draw_sol) {
180 0 : PetscCall(PetscViewerDrawSetPause(PETSC_VIEWER_DRAW_WORLD,-1));
181 15 : PetscCall(VecView(x,PETSC_VIEWER_DRAW_WORLD));
182 : }
183 : }
184 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
185 : }
186 :
187 5 : PetscCall(NEPDestroy(&nep));
188 5 : PetscCall(MatDestroy(&F));
189 5 : PetscCall(MatDestroy(&J));
190 5 : PetscCall(VecDestroy(&x));
191 5 : PetscCall(SlepcFinalize());
192 : return 0;
193 : }
194 :
195 : /* ------------------------------------------------------------------- */
196 : /*
197 : FormInitialGuess - Computes initial guess.
198 :
199 : Input/Output Parameter:
200 : . x - the solution vector
201 : */
202 5 : PetscErrorCode FormInitialGuess(Vec x)
203 : {
204 5 : PetscFunctionBeginUser;
205 5 : PetscCall(VecSet(x,1.0));
206 5 : PetscFunctionReturn(PETSC_SUCCESS);
207 : }
208 :
209 : /* ------------------------------------------------------------------- */
210 : /*
211 : FormFunction - Computes Function matrix T(lambda)
212 :
213 : Input Parameters:
214 : . nep - the NEP context
215 : . lambda - the scalar argument
216 : . ctx - optional user-defined context, as set by NEPSetFunction()
217 :
218 : Output Parameters:
219 : . fun - Function matrix
220 : . B - optionally different preconditioning matrix
221 : */
222 415 : PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
223 : {
224 415 : ApplicationCtx *user = (ApplicationCtx*)ctx;
225 415 : PetscScalar A[3],c,d;
226 415 : PetscReal h;
227 415 : PetscInt i,n,j[3],Istart,Iend;
228 415 : PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
229 :
230 415 : PetscFunctionBeginUser;
231 : /*
232 : Compute Function entries and insert into matrix
233 : */
234 415 : PetscCall(MatGetSize(fun,&n,NULL));
235 415 : PetscCall(MatGetOwnershipRange(fun,&Istart,&Iend));
236 415 : if (Istart==0) FirstBlock=PETSC_TRUE;
237 415 : if (Iend==n) LastBlock=PETSC_TRUE;
238 415 : h = user->h;
239 415 : c = user->kappa/(lambda-user->kappa);
240 415 : d = n;
241 :
242 : /*
243 : Interior grid points
244 : */
245 62049 : for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
246 61634 : j[0] = i-1; j[1] = i; j[2] = i+1;
247 61634 : A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
248 61634 : PetscCall(MatSetValues(fun,1,&i,3,j,A,INSERT_VALUES));
249 : }
250 :
251 : /*
252 : Boundary points
253 : */
254 415 : if (FirstBlock) {
255 415 : i = 0;
256 415 : j[0] = 0; j[1] = 1;
257 415 : A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
258 415 : PetscCall(MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES));
259 : }
260 :
261 415 : if (LastBlock) {
262 415 : i = n-1;
263 415 : j[0] = n-2; j[1] = n-1;
264 415 : A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
265 415 : PetscCall(MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES));
266 : }
267 :
268 : /*
269 : Assemble matrix
270 : */
271 415 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
272 415 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
273 415 : if (fun != B) {
274 0 : PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY));
275 0 : PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY));
276 : }
277 415 : PetscFunctionReturn(PETSC_SUCCESS);
278 : }
279 :
280 : /* ------------------------------------------------------------------- */
281 : /*
282 : FormJacobian - Computes Jacobian matrix T'(lambda)
283 :
284 : Input Parameters:
285 : . nep - the NEP context
286 : . lambda - the scalar argument
287 : . ctx - optional user-defined context, as set by NEPSetJacobian()
288 :
289 : Output Parameters:
290 : . jac - Jacobian matrix
291 : */
292 79 : PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
293 : {
294 79 : ApplicationCtx *user = (ApplicationCtx*)ctx;
295 79 : PetscScalar A[3],c;
296 79 : PetscReal h;
297 79 : PetscInt i,n,j[3],Istart,Iend;
298 79 : PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
299 :
300 79 : PetscFunctionBeginUser;
301 : /*
302 : Compute Jacobian entries and insert into matrix
303 : */
304 79 : PetscCall(MatGetSize(jac,&n,NULL));
305 79 : PetscCall(MatGetOwnershipRange(jac,&Istart,&Iend));
306 79 : if (Istart==0) FirstBlock=PETSC_TRUE;
307 79 : if (Iend==n) LastBlock=PETSC_TRUE;
308 79 : h = user->h;
309 79 : c = user->kappa/(lambda-user->kappa);
310 :
311 : /*
312 : Interior grid points
313 : */
314 18737 : for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
315 18658 : j[0] = i-1; j[1] = i; j[2] = i+1;
316 18658 : A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
317 18658 : PetscCall(MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES));
318 : }
319 :
320 : /*
321 : Boundary points
322 : */
323 79 : if (FirstBlock) {
324 79 : i = 0;
325 79 : j[0] = 0; j[1] = 1;
326 79 : A[0] = -2.0*h/3.0; A[1] = -h/6.0;
327 79 : PetscCall(MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES));
328 : }
329 :
330 79 : if (LastBlock) {
331 79 : i = n-1;
332 79 : j[0] = n-2; j[1] = n-1;
333 79 : A[0] = -h/6.0; A[1] = -h/3.0-c*c;
334 79 : PetscCall(MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES));
335 : }
336 :
337 : /*
338 : Assemble matrix
339 : */
340 79 : PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY));
341 79 : PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY));
342 79 : PetscFunctionReturn(PETSC_SUCCESS);
343 : }
344 :
345 : /* ------------------------------------------------------------------- */
346 : /*
347 : CheckSolution - Given a computed solution (lambda,x) check if it
348 : satisfies the analytic solution.
349 :
350 : Input Parameters:
351 : + lambda - the computed eigenvalue
352 : - y - the computed eigenvector
353 :
354 : Output Parameter:
355 : . error - norm of difference between the computed and exact eigenvector
356 : */
357 15 : PetscErrorCode CheckSolution(PetscScalar lambda,Vec y,PetscReal *error,void *ctx)
358 : {
359 15 : PetscScalar *uu;
360 15 : PetscInt i,n,Istart,Iend;
361 15 : PetscReal nu,x;
362 15 : Vec u;
363 15 : ApplicationCtx *user = (ApplicationCtx*)ctx;
364 :
365 15 : PetscFunctionBeginUser;
366 15 : nu = PetscSqrtReal(PetscRealPart(lambda));
367 15 : PetscCall(VecDuplicate(y,&u));
368 15 : PetscCall(VecGetSize(u,&n));
369 15 : PetscCall(VecGetOwnershipRange(y,&Istart,&Iend));
370 15 : PetscCall(VecGetArray(u,&uu));
371 2191 : for (i=Istart;i<Iend;i++) {
372 2176 : x = (i+1)*user->h;
373 2176 : uu[i-Istart] = PetscSinReal(nu*x);
374 : }
375 15 : PetscCall(VecRestoreArray(u,&uu));
376 15 : PetscCall(VecNormalize(u,NULL));
377 15 : PetscCall(VecAXPY(u,-1.0,y));
378 15 : PetscCall(VecNorm(u,NORM_2,error));
379 15 : PetscCall(VecDestroy(&u));
380 15 : PetscFunctionReturn(PETSC_SUCCESS);
381 : }
382 :
383 : /* ------------------------------------------------------------------- */
384 : /*
385 : FixSign - Force the eigenfunction to be real and positive, since
386 : some eigensolvers may return the eigenvector multiplied by a
387 : complex number of modulus one.
388 :
389 : Input/Output Parameter:
390 : . x - the computed vector
391 : */
392 15 : PetscErrorCode FixSign(Vec x)
393 : {
394 15 : PetscMPIInt rank;
395 15 : PetscScalar sign=0.0;
396 15 : const PetscScalar *xx;
397 :
398 15 : PetscFunctionBeginUser;
399 15 : PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
400 15 : if (!rank) {
401 15 : PetscCall(VecGetArrayRead(x,&xx));
402 15 : sign = *xx/PetscAbsScalar(*xx);
403 15 : PetscCall(VecRestoreArrayRead(x,&xx));
404 : }
405 30 : PetscCallMPI(MPI_Bcast(&sign,1,MPIU_SCALAR,0,PETSC_COMM_WORLD));
406 15 : PetscCall(VecScale(x,1.0/sign));
407 15 : PetscFunctionReturn(PETSC_SUCCESS);
408 : }
409 :
410 : /*TEST
411 :
412 : build:
413 : requires: !single
414 :
415 : test:
416 : suffix: 1
417 : args: -nep_target 4
418 : filter: sed -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
419 : requires: !single
420 :
421 : testset:
422 : args: -nep_type nleigs -nep_target 10 -nep_nev 4
423 : test:
424 : suffix: 2
425 : filter: sed -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
426 : args: -rg_interval_endpoints 4,200
427 : requires: !single !complex
428 : test:
429 : suffix: 2_complex
430 : filter: sed -e "s/[+-]0.[0-9]*i//" -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
431 : output_file: output/ex20_2.out
432 : args: -rg_interval_endpoints 4,200,-.1,.1 -nep_target_real
433 : requires: !single complex
434 :
435 : testset:
436 : args: -nep_type nleigs -nep_target 10 -nep_nev 4 -nep_ncv 18 -nep_two_sided {{0 1}} -nep_nleigs_full_basis
437 : test:
438 : suffix: 3
439 : filter: sed -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
440 : output_file: output/ex20_2.out
441 : args: -rg_interval_endpoints 4,200
442 : requires: !single !complex
443 : test:
444 : suffix: 3_complex
445 : filter: sed -e "s/[+-]0.[0-9]*i//" -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
446 : output_file: output/ex20_2.out
447 : args: -rg_interval_endpoints 4,200,-.1,.1
448 : requires: !single complex
449 :
450 : test:
451 : suffix: 4
452 : args: -n 256 -nep_nev 2 -nep_target 10
453 : filter: sed -e "s/[+-]0\.0*i//g" -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
454 : requires: !single
455 :
456 : TEST*/
|
