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 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";
15 :
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 : */
22 :
23 : #include <slepcnep.h>
24 :
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*);
31 :
32 20 : int main(int argc,char **argv)
33 : {
34 20 : NEP nep; /* nonlinear eigensolver context */
35 20 : Mat F,J,A[2];
36 20 : NEPType type;
37 20 : PetscInt n=100,nev,Istart,Iend,i;
38 20 : PetscBool terse,split=PETSC_TRUE;
39 20 : RG rg;
40 20 : FN f[2];
41 20 : PetscScalar coeffs;
42 :
43 20 : PetscFunctionBeginUser;
44 20 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
45 20 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
46 20 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL));
47 31 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSquare root eigenproblem, n=%" PetscInt_FMT "%s\n\n",n,split?" (in split form)":""));
48 :
49 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50 : Create nonlinear eigensolver context
51 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52 :
53 20 : PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
54 :
55 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56 : Select the NLEIGS solver and set required options for it
57 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58 :
59 20 : PetscCall(NEPSetType(nep,NEPNLEIGS));
60 20 : PetscCall(NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,NULL));
61 20 : PetscCall(NEPGetRG(nep,&rg));
62 20 : PetscCall(RGSetType(rg,RGINTERVAL));
63 : #if defined(PETSC_USE_COMPLEX)
64 20 : PetscCall(RGIntervalSetEndpoints(rg,0.01,16.0,-0.001,0.001));
65 : #else
66 : PetscCall(RGIntervalSetEndpoints(rg,0.01,16.0,0,0));
67 : #endif
68 20 : PetscCall(NEPSetTarget(nep,1.1));
69 :
70 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71 : Define the nonlinear problem
72 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73 :
74 20 : if (split) {
75 : /*
76 : Create matrices for the split form
77 : */
78 9 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A[0]));
79 9 : PetscCall(MatSetSizes(A[0],PETSC_DECIDE,PETSC_DECIDE,n,n));
80 9 : PetscCall(MatSetFromOptions(A[0]));
81 9 : PetscCall(MatGetOwnershipRange(A[0],&Istart,&Iend));
82 529 : for (i=Istart;i<Iend;i++) {
83 520 : if (i>0) PetscCall(MatSetValue(A[0],i,i-1,1.0,INSERT_VALUES));
84 520 : if (i<n-1) PetscCall(MatSetValue(A[0],i,i+1,1.0,INSERT_VALUES));
85 520 : PetscCall(MatSetValue(A[0],i,i,-2.0,INSERT_VALUES));
86 : }
87 9 : PetscCall(MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY));
88 9 : PetscCall(MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY));
89 :
90 9 : PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,1.0,&A[1]));
91 :
92 : /*
93 : Define functions for the split form
94 : */
95 9 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f[0]));
96 9 : PetscCall(FNSetType(f[0],FNRATIONAL));
97 9 : coeffs = 1.0;
98 9 : PetscCall(FNRationalSetNumerator(f[0],1,&coeffs));
99 9 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f[1]));
100 9 : PetscCall(FNSetType(f[1],FNSQRT));
101 9 : PetscCall(NEPSetSplitOperator(nep,2,A,f,SUBSET_NONZERO_PATTERN));
102 :
103 : } else {
104 : /*
105 : Callback form: create matrix and set Function evaluation routine
106 : */
107 11 : PetscCall(MatCreate(PETSC_COMM_WORLD,&F));
108 11 : PetscCall(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n));
109 11 : PetscCall(MatSetFromOptions(F));
110 11 : PetscCall(MatSeqAIJSetPreallocation(F,3,NULL));
111 11 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
112 11 : PetscCall(NEPSetFunction(nep,F,F,FormFunction,NULL));
113 :
114 11 : PetscCall(MatCreate(PETSC_COMM_WORLD,&J));
115 11 : PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n));
116 11 : PetscCall(MatSetFromOptions(J));
117 11 : PetscCall(MatSeqAIJSetPreallocation(J,1,NULL));
118 11 : PetscCall(MatMPIAIJSetPreallocation(J,1,NULL,1,NULL));
119 11 : PetscCall(NEPSetJacobian(nep,J,FormJacobian,NULL));
120 : }
121 :
122 : /*
123 : Set solver parameters at runtime
124 : */
125 20 : PetscCall(NEPSetFromOptions(nep));
126 :
127 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128 : Solve the eigensystem
129 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130 20 : PetscCall(NEPSolve(nep));
131 20 : PetscCall(NEPGetType(nep,&type));
132 20 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type));
133 20 : PetscCall(NEPGetDimensions(nep,&nev,NULL,NULL));
134 20 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
135 :
136 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137 : Display solution and clean up
138 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139 :
140 : /* show detailed info unless -terse option is given by user */
141 20 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
142 20 : if (terse) PetscCall(NEPErrorView(nep,NEP_ERROR_BACKWARD,NULL));
143 : else {
144 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
145 0 : PetscCall(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD));
146 0 : PetscCall(NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD));
147 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
148 : }
149 20 : PetscCall(NEPDestroy(&nep));
150 20 : if (split) {
151 9 : PetscCall(MatDestroy(&A[0]));
152 9 : PetscCall(MatDestroy(&A[1]));
153 9 : PetscCall(FNDestroy(&f[0]));
154 9 : PetscCall(FNDestroy(&f[1]));
155 : } else {
156 11 : PetscCall(MatDestroy(&F));
157 11 : PetscCall(MatDestroy(&J));
158 : }
159 20 : PetscCall(SlepcFinalize());
160 : return 0;
161 : }
162 :
163 : /* ------------------------------------------------------------------- */
164 : /*
165 : FormFunction - Computes Function matrix T(lambda)
166 : */
167 564 : PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
168 : {
169 564 : PetscInt i,n,col[3],Istart,Iend;
170 564 : PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
171 564 : PetscScalar value[3],t;
172 :
173 564 : PetscFunctionBeginUser;
174 : /*
175 : Compute Function entries and insert into matrix
176 : */
177 564 : t = PetscSqrtScalar(lambda);
178 564 : PetscCall(MatGetSize(fun,&n,NULL));
179 564 : PetscCall(MatGetOwnershipRange(fun,&Istart,&Iend));
180 564 : if (Istart==0) FirstBlock=PETSC_TRUE;
181 564 : if (Iend==n) LastBlock=PETSC_TRUE;
182 564 : value[0]=1.0; value[1]=t-2.0; value[2]=1.0;
183 45350 : for (i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++) {
184 44786 : col[0]=i-1; col[1]=i; col[2]=i+1;
185 44786 : PetscCall(MatSetValues(fun,1,&i,3,col,value,INSERT_VALUES));
186 : }
187 564 : if (LastBlock) {
188 457 : i=n-1; col[0]=n-2; col[1]=n-1;
189 457 : PetscCall(MatSetValues(fun,1,&i,2,col,value,INSERT_VALUES));
190 : }
191 564 : if (FirstBlock) {
192 457 : i=0; col[0]=0; col[1]=1; value[0]=t-2.0; value[1]=1.0;
193 457 : PetscCall(MatSetValues(fun,1,&i,2,col,value,INSERT_VALUES));
194 : }
195 :
196 : /*
197 : Assemble matrix
198 : */
199 564 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
200 564 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
201 564 : if (fun != B) {
202 0 : PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY));
203 0 : PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY));
204 : }
205 564 : PetscFunctionReturn(PETSC_SUCCESS);
206 : }
207 :
208 : /* ------------------------------------------------------------------- */
209 : /*
210 : FormJacobian - Computes Jacobian matrix T'(lambda)
211 : */
212 192 : PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
213 : {
214 192 : Vec d;
215 :
216 192 : PetscFunctionBeginUser;
217 192 : PetscCall(MatCreateVecs(jac,&d,NULL));
218 192 : PetscCall(VecSet(d,0.5/PetscSqrtScalar(lambda)));
219 192 : PetscCall(MatDiagonalSet(jac,d,INSERT_VALUES));
220 192 : PetscCall(VecDestroy(&d));
221 192 : PetscFunctionReturn(PETSC_SUCCESS);
222 : }
223 :
224 : /* ------------------------------------------------------------------- */
225 : /*
226 : ComputeSingularities - Computes maxnp points (at most) in the complex plane where
227 : the function T(.) is not analytic.
228 :
229 : In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
230 : */
231 5 : PetscErrorCode ComputeSingularities(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *pt)
232 : {
233 5 : PetscReal h;
234 5 : PetscInt i;
235 :
236 5 : PetscFunctionBeginUser;
237 5 : h = 11.0/(*maxnp-1);
238 5 : xi[0] = -1e-5; xi[*maxnp-1] = -1e+6;
239 49995 : for (i=1;i<*maxnp-1;i++) xi[i] = -PetscPowReal(10,-5+h*i);
240 5 : PetscFunctionReturn(PETSC_SUCCESS);
241 : }
242 :
243 : /*TEST
244 :
245 : testset:
246 : args: -nep_nev 3 -terse
247 : output_file: output/ex27_1.out
248 : requires: !single
249 : filter: sed -e "s/[+-]0\.0*i//g"
250 : test:
251 : suffix: 1
252 : args: -nep_nleigs_interpolation_degree 90
253 : test:
254 : suffix: 3
255 : args: -nep_tol 1e-8 -nep_nleigs_rk_shifts 1.06,1.1,1.12,1.15 -nep_conv_norm -nep_nleigs_interpolation_degree 20
256 : test:
257 : suffix: 5
258 : args: -mat_type aijcusparse
259 : requires: cuda
260 :
261 : testset:
262 : args: -split 0 -nep_nev 3 -terse
263 : output_file: output/ex27_2.out
264 : filter: sed -e "s/[+-]0\.0*i//g"
265 : test:
266 : suffix: 2
267 : args: -nep_nleigs_interpolation_degree 90
268 : requires: !single
269 : test:
270 : suffix: 4
271 : args: -nep_nleigs_rk_shifts 1.06,1.1,1.12,1.15 -nep_nleigs_interpolation_degree 20
272 : requires: double
273 : test:
274 : suffix: 6
275 : args: -mat_type aijcusparse
276 : requires: cuda !single
277 :
278 : testset:
279 : args: -split 0 -nep_type ciss -nep_ciss_extraction {{ritz hankel caa}} -rg_type ellipse -rg_ellipse_center 8 -rg_ellipse_radius .7 -nep_ciss_moments 4 -rg_ellipse_vscale 0.1 -terse
280 : requires: complex !single
281 : output_file: output/ex27_7.out
282 : timeoutfactor: 2
283 : test:
284 : suffix: 7
285 : test:
286 : suffix: 7_par
287 : nsize: 2
288 : args: -nep_ciss_partitions 2
289 :
290 : testset:
291 : args: -nep_type ciss -rg_type ellipse -rg_ellipse_center 8 -rg_ellipse_radius .7 -rg_ellipse_vscale 0.1 -terse
292 : requires: complex
293 : filter: sed -e "s/ (in split form)//" | sed -e "s/56925/56924/" | sed -e "s/60753/60754/" | sed -e "s/92630/92629/" | sed -e "s/24705/24706/"
294 : output_file: output/ex27_7.out
295 : timeoutfactor: 2
296 : test:
297 : suffix: 8
298 : test:
299 : suffix: 8_parallel
300 : nsize: 4
301 : args: -nep_ciss_partitions 4 -ds_parallel distributed
302 : test:
303 : suffix: 8_hpddm
304 : args: -nep_ciss_ksp_type hpddm
305 : requires: hpddm
306 :
307 : test:
308 : suffix: 9
309 : args: -nep_nev 4 -n 20 -terse
310 : requires: !single
311 : filter: sed -e "s/[+-]0\.0*i//g"
312 :
313 : TEST*/
|
