Actual source code: ex27.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: static char help[] = "Simple nonlinear eigenproblem using the NLEIGS solver.\n\n"
12: "The command line options are:\n"
13: " -n <n>, where <n> = matrix dimension.\n"
14: " -split <0/1>, to select the split form in the problem definition (enabled by default)\n";
16: /*
17: Solve T(lambda)x=0 using NLEIGS solver
18: with T(lambda) = -D+sqrt(lambda)*I
19: where D is the Laplacian operator in 1 dimension
20: and with the interpolation interval [.01,16]
21: */
23: #include <slepcnep.h>
25: /*
26: User-defined routines
27: */
28: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
29: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);
30: PetscErrorCode ComputeSingularities(NEP,PetscInt*,PetscScalar*,void*);
32: int main(int argc,char **argv)
33: {
34: NEP nep; /* nonlinear eigensolver context */
35: Mat F,J,A[2];
36: NEPType type;
37: PetscInt n=100,nev,Istart,Iend,i;
39: PetscBool terse,split=PETSC_TRUE;
40: RG rg;
41: FN f[2];
42: PetscScalar coeffs;
44: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
45: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
46: PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL);
47: PetscPrintf(PETSC_COMM_WORLD,"\nSquare root eigenproblem, n=%D%s\n\n",n,split?" (in split form)":"");
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Create nonlinear eigensolver context
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53: NEPCreate(PETSC_COMM_WORLD,&nep);
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Select the NLEIGS solver and set required options for it
57: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: NEPSetType(nep,NEPNLEIGS);
60: NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,NULL);
61: NEPGetRG(nep,&rg);
62: RGSetType(rg,RGINTERVAL);
63: #if defined(PETSC_USE_COMPLEX)
64: RGIntervalSetEndpoints(rg,0.01,16.0,-0.001,0.001);
65: #else
66: RGIntervalSetEndpoints(rg,0.01,16.0,0,0);
67: #endif
68: NEPSetTarget(nep,1.1);
70: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: Define the nonlinear problem
72: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74: if (split) {
75: /*
76: Create matrices for the split form
77: */
78: MatCreate(PETSC_COMM_WORLD,&A[0]);
79: MatSetSizes(A[0],PETSC_DECIDE,PETSC_DECIDE,n,n);
80: MatSetFromOptions(A[0]);
81: MatSetUp(A[0]);
82: MatGetOwnershipRange(A[0],&Istart,&Iend);
83: for (i=Istart;i<Iend;i++) {
84: if (i>0) { MatSetValue(A[0],i,i-1,1.0,INSERT_VALUES); }
85: if (i<n-1) { MatSetValue(A[0],i,i+1,1.0,INSERT_VALUES); }
86: MatSetValue(A[0],i,i,-2.0,INSERT_VALUES);
87: }
88: MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY);
89: MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY);
91: MatCreate(PETSC_COMM_WORLD,&A[1]);
92: MatSetSizes(A[1],PETSC_DECIDE,PETSC_DECIDE,n,n);
93: MatSetFromOptions(A[1]);
94: MatSetUp(A[1]);
95: MatAssemblyBegin(A[1],MAT_FINAL_ASSEMBLY);
96: MatAssemblyEnd(A[1],MAT_FINAL_ASSEMBLY);
97: MatShift(A[1],1.0);
99: /*
100: Define funcions for the split form
101: */
102: FNCreate(PETSC_COMM_WORLD,&f[0]);
103: FNSetType(f[0],FNRATIONAL);
104: coeffs = 1.0;
105: FNRationalSetNumerator(f[0],1,&coeffs);
106: FNCreate(PETSC_COMM_WORLD,&f[1]);
107: FNSetType(f[1],FNSQRT);
108: NEPSetSplitOperator(nep,2,A,f,SUBSET_NONZERO_PATTERN);
110: } else {
111: /*
112: Callback form: create matrix and set Function evaluation routine
113: */
114: MatCreate(PETSC_COMM_WORLD,&F);
115: MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n);
116: MatSetFromOptions(F);
117: MatSeqAIJSetPreallocation(F,3,NULL);
118: MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
119: MatSetUp(F);
120: NEPSetFunction(nep,F,F,FormFunction,NULL);
122: MatCreate(PETSC_COMM_WORLD,&J);
123: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
124: MatSetFromOptions(J);
125: MatSeqAIJSetPreallocation(J,1,NULL);
126: MatMPIAIJSetPreallocation(J,1,NULL,1,NULL);
127: MatSetUp(J);
128: NEPSetJacobian(nep,J,FormJacobian,NULL);
129: }
131: /*
132: Set solver parameters at runtime
133: */
134: NEPSetFromOptions(nep);
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Solve the eigensystem
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: NEPSolve(nep);
140: NEPGetType(nep,&type);
141: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
142: NEPGetDimensions(nep,&nev,NULL,NULL);
143: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Display solution and clean up
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: /* show detailed info unless -terse option is given by user */
150: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
151: if (terse) {
152: NEPErrorView(nep,NEP_ERROR_BACKWARD,NULL);
153: } else {
154: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
155: NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
156: NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
157: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
158: }
159: NEPDestroy(&nep);
160: if (split) {
161: MatDestroy(&A[0]);
162: MatDestroy(&A[1]);
163: FNDestroy(&f[0]);
164: FNDestroy(&f[1]);
165: } else {
166: MatDestroy(&F);
167: MatDestroy(&J);
168: }
169: SlepcFinalize();
170: return ierr;
171: }
173: /* ------------------------------------------------------------------- */
174: /*
175: FormFunction - Computes Function matrix T(lambda)
176: */
177: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
178: {
180: PetscInt i,n,col[3],Istart,Iend;
181: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
182: PetscScalar value[3],t;
185: /*
186: Compute Function entries and insert into matrix
187: */
188: t = PetscSqrtScalar(lambda);
189: MatGetSize(fun,&n,NULL);
190: MatGetOwnershipRange(fun,&Istart,&Iend);
191: if (Istart==0) FirstBlock=PETSC_TRUE;
192: if (Iend==n) LastBlock=PETSC_TRUE;
193: value[0]=1.0; value[1]=t-2.0; value[2]=1.0;
194: for (i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++) {
195: col[0]=i-1; col[1]=i; col[2]=i+1;
196: MatSetValues(fun,1,&i,3,col,value,INSERT_VALUES);
197: }
198: if (LastBlock) {
199: i=n-1; col[0]=n-2; col[1]=n-1;
200: MatSetValues(fun,1,&i,2,col,value,INSERT_VALUES);
201: }
202: if (FirstBlock) {
203: i=0; col[0]=0; col[1]=1; value[0]=t-2.0; value[1]=1.0;
204: MatSetValues(fun,1,&i,2,col,value,INSERT_VALUES);
205: }
207: /*
208: Assemble matrix
209: */
210: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
211: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
212: if (fun != B) {
213: MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
214: MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
215: }
216: return(0);
217: }
219: /* ------------------------------------------------------------------- */
220: /*
221: FormJacobian - Computes Jacobian matrix T'(lambda)
222: */
223: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
224: {
226: Vec d;
229: MatCreateVecs(jac,&d,NULL);
230: VecSet(d,0.5/PetscSqrtScalar(lambda));
231: MatDiagonalSet(jac,d,INSERT_VALUES);
232: VecDestroy(&d);
233: return(0);
234: }
236: /* ------------------------------------------------------------------- */
237: /*
238: ComputeSingularities - Computes maxnp points (at most) in the complex plane where
239: the function T(.) is not analytic.
241: In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
242: */
243: PetscErrorCode ComputeSingularities(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *pt)
244: {
245: PetscReal h;
246: PetscInt i;
249: h = 11.0/(*maxnp-1);
250: xi[0] = -1e-5; xi[*maxnp-1] = -1e+6;
251: for (i=1;i<*maxnp-1;i++) xi[i] = -PetscPowReal(10,-5+h*i);
252: return(0);
253: }
255: /*TEST
257: test:
258: suffix: 1
259: args: -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
261: test:
262: suffix: 2
263: args: -split 0 -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
265: test:
266: suffix: 3
267: args: -nep_nev 3 -nep_tol 1e-8 -nep_nleigs_rk_shifts 1.06,1.1,1.12,1.15 -nep_conv_norm -terse
268: requires: !single
269: output_file: output/ex27_1.out
271: test:
272: suffix: 4
273: args: -split 0 -nep_nev 3 -nep_nleigs_rk_shifts 1.06,1.1,1.12,1.15 -nep_nleigs_interpolation_degree 90 -terse
274: output_file: output/ex27_2.out
276: test:
277: suffix: 5
278: args: -nep_nev 3 -mat_type aijcusparse -terse
279: requires: cuda
280: output_file: output/ex27_1.out
282: test:
283: suffix: 6
284: args: -split 0 -nep_nev 3 -mat_type aijcusparse -terse
285: requires: cuda
286: output_file: output/ex27_2.out
288: test:
289: suffix: 7
290: args: -split 0 -nep_type ciss -rg_type ellipse -rg_ellipse_center 8 -rg_ellipse_radius .7 -rg_ellipse_vscale 0.1 -terse
291: requires: complex
293: test:
294: suffix: 8
295: args: -nep_type ciss -rg_type ellipse -rg_ellipse_center 8 -rg_ellipse_radius .7 -rg_ellipse_vscale 0.1 -terse
296: requires: complex
297: filter: sed -e "s/ (in split form)//"
298: output_file: output/ex27_7.out
300: test:
301: suffix: 9
302: args: -nep_nev 4 -n 20 -terse
303: requires: !single
305: TEST*/