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 : This example implements one of the problems found at
12 : NLEVP: A Collection of Nonlinear Eigenvalue Problems,
13 : The University of Manchester.
14 : The details of the collection can be found at:
15 : [1] T. Betcke et al., "NLEVP: A Collection of Nonlinear Eigenvalue
16 : Problems", ACM Trans. Math. Software 39(2), Article 7, 2013.
17 :
18 : The loaded_string problem is a rational eigenvalue problem for the
19 : finite element model of a loaded vibrating string.
20 : */
21 :
22 : static char help[] = "Finite element model of a loaded vibrating string.\n\n"
23 : "The command line options are:\n"
24 : " -n <n>, dimension of the matrices.\n"
25 : " -kappa <kappa>, stiffness of elastic spring.\n"
26 : " -mass <m>, mass of the attached load.\n\n";
27 :
28 : #include <slepcnep.h>
29 :
30 : #define NMAT 3
31 :
32 21 : int main(int argc,char **argv)
33 : {
34 21 : Mat A[NMAT]; /* problem matrices */
35 21 : FN f[NMAT]; /* functions to define the nonlinear operator */
36 21 : NEP nep; /* nonlinear eigensolver context */
37 21 : PetscInt n=100,Istart,Iend,i;
38 21 : PetscReal kappa=1.0,m=1.0;
39 21 : PetscScalar sigma,numer[2],denom[2];
40 21 : PetscBool terse;
41 :
42 21 : PetscFunctionBeginUser;
43 21 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
44 :
45 21 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
46 21 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-kappa",&kappa,NULL));
47 21 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-mass",&m,NULL));
48 21 : sigma = kappa/m;
49 21 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Loaded vibrating string, n=%" PetscInt_FMT " kappa=%g m=%g\n\n",n,(double)kappa,(double)m));
50 :
51 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52 : Build the problem matrices
53 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54 :
55 : /* initialize matrices */
56 84 : for (i=0;i<NMAT;i++) {
57 63 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A[i]));
58 63 : PetscCall(MatSetSizes(A[i],PETSC_DECIDE,PETSC_DECIDE,n,n));
59 63 : PetscCall(MatSetFromOptions(A[i]));
60 : }
61 :
62 : /* A0 */
63 21 : PetscCall(MatGetOwnershipRange(A[0],&Istart,&Iend));
64 1541 : for (i=Istart;i<Iend;i++) {
65 1520 : PetscCall(MatSetValue(A[0],i,i,(i==n-1)?1.0*n:2.0*n,INSERT_VALUES));
66 1520 : if (i>0) PetscCall(MatSetValue(A[0],i,i-1,-1.0*n,INSERT_VALUES));
67 1520 : if (i<n-1) PetscCall(MatSetValue(A[0],i,i+1,-1.0*n,INSERT_VALUES));
68 : }
69 :
70 : /* A1 */
71 21 : PetscCall(MatGetOwnershipRange(A[1],&Istart,&Iend));
72 1541 : for (i=Istart;i<Iend;i++) {
73 1520 : PetscCall(MatSetValue(A[1],i,i,(i==n-1)?2.0/(6.0*n):4.0/(6.0*n),INSERT_VALUES));
74 1520 : if (i>0) PetscCall(MatSetValue(A[1],i,i-1,1.0/(6.0*n),INSERT_VALUES));
75 1520 : if (i<n-1) PetscCall(MatSetValue(A[1],i,i+1,1.0/(6.0*n),INSERT_VALUES));
76 : }
77 :
78 : /* A2 */
79 21 : PetscCall(MatGetOwnershipRange(A[2],&Istart,&Iend));
80 21 : if (Istart<=n-1 && n-1<Iend) PetscCall(MatSetValue(A[2],n-1,n-1,kappa,INSERT_VALUES));
81 :
82 : /* assemble matrices */
83 84 : for (i=0;i<NMAT;i++) PetscCall(MatAssemblyBegin(A[i],MAT_FINAL_ASSEMBLY));
84 84 : for (i=0;i<NMAT;i++) PetscCall(MatAssemblyEnd(A[i],MAT_FINAL_ASSEMBLY));
85 :
86 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87 : Create the problem functions
88 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89 :
90 : /* f1=1 */
91 21 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f[0]));
92 21 : PetscCall(FNSetType(f[0],FNRATIONAL));
93 21 : numer[0] = 1.0;
94 21 : PetscCall(FNRationalSetNumerator(f[0],1,numer));
95 :
96 : /* f2=-lambda */
97 21 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f[1]));
98 21 : PetscCall(FNSetType(f[1],FNRATIONAL));
99 21 : numer[0] = -1.0; numer[1] = 0.0;
100 21 : PetscCall(FNRationalSetNumerator(f[1],2,numer));
101 :
102 : /* f3=lambda/(lambda-sigma) */
103 21 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f[2]));
104 21 : PetscCall(FNSetType(f[2],FNRATIONAL));
105 21 : numer[0] = 1.0; numer[1] = 0.0;
106 21 : denom[0] = 1.0; denom[1] = -sigma;
107 21 : PetscCall(FNRationalSetNumerator(f[2],2,numer));
108 21 : PetscCall(FNRationalSetDenominator(f[2],2,denom));
109 :
110 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111 : Create the eigensolver and solve the problem
112 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113 :
114 21 : PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
115 21 : PetscCall(NEPSetSplitOperator(nep,3,A,f,SUBSET_NONZERO_PATTERN));
116 21 : PetscCall(NEPSetProblemType(nep,NEP_RATIONAL));
117 21 : PetscCall(NEPSetFromOptions(nep));
118 21 : PetscCall(NEPSolve(nep));
119 :
120 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121 : Display solution and clean up
122 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123 :
124 : /* show detailed info unless -terse option is given by user */
125 21 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
126 21 : if (terse) PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
127 : else {
128 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
129 0 : PetscCall(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD));
130 0 : PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
131 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
132 : }
133 21 : PetscCall(NEPDestroy(&nep));
134 84 : for (i=0;i<NMAT;i++) {
135 63 : PetscCall(MatDestroy(&A[i]));
136 63 : PetscCall(FNDestroy(&f[i]));
137 : }
138 21 : PetscCall(SlepcFinalize());
139 : return 0;
140 : }
141 :
142 : /*TEST
143 :
144 : test:
145 : suffix: 1
146 : args: -nep_type rii -nep_target 4 -terse
147 : requires: !single
148 : filter: sed -e "s/[+-]0\.0*i//g"
149 :
150 : testset:
151 : args: -nep_type interpol -rg_type interval -rg_interval_endpoints 5,700,-.1,.1 -nep_nev 7 -nep_target 5 -nep_interpol_interpolation_degree 12 -nep_refine simple -terse
152 : requires: !single
153 : output_file: output/loaded_string_2.out
154 : test:
155 : suffix: 2
156 : args: -nep_refine_scheme {{schur explicit}}
157 : test:
158 : suffix: 2_mbe
159 : args: -nep_refine_scheme mbe -nep_refine_ksp_type preonly -nep_refine_pc_type lu
160 :
161 : testset:
162 : nsize: 2
163 : args: -nep_type interpol -rg_type interval -rg_interval_endpoints 5,700,-.1,.1 -nep_nev 7 -nep_target 5 -nep_interpol_interpolation_degree 12 -nep_refine simple -nep_refine_partitions 2 -nep_interpol_st_ksp_type bcgs -nep_interpol_st_pc_type bjacobi -terse
164 : requires: !single
165 : output_file: output/loaded_string_2.out
166 : timeoutfactor: 2
167 : test:
168 : suffix: 3_explicit
169 : args: -nep_refine_scheme explicit
170 : test:
171 : suffix: 3_mbe
172 : args: -nep_refine_scheme mbe -nep_refine_ksp_type preonly -nep_refine_pc_type cholesky
173 :
174 : test:
175 : suffix: 4
176 : nsize: 4
177 : args: -nep_type interpol -rg_type interval -rg_interval_endpoints 5,700,-.1,.1 -nep_nev 7 -nep_target 5 -nep_interpol_interpolation_degree 10 -nep_refine simple -nep_refine_partitions 2 -nep_refine_scheme explicit -nep_interpol_st_ksp_type bcgs -nep_interpol_st_pc_type bjacobi -terse -log_exclude nep,pep,fn
178 : requires: !single
179 : output_file: output/loaded_string_2.out
180 : timeoutfactor: 4
181 :
182 : test:
183 : suffix: 5
184 : args: -nep_type nleigs -rg_type interval -rg_interval_endpoints 4,700,-.1,.1 -nep_nev 8 -nep_target 5 -terse
185 : filter: sed -e "s/[+-]0\.0*i//g"
186 : requires: !single
187 :
188 : test:
189 : suffix: 6
190 : args: -nep_type nleigs -rg_type interval -rg_interval_endpoints 100,700 -nep_nev 5 -nep_tol 1e-9 -nep_target 140 -nep_nleigs_interpolation_degree 15 -nep_general -terse
191 : requires: !complex !single
192 :
193 : test:
194 : suffix: 6_complex
195 : args: -nep_type nleigs -rg_type interval -rg_interval_endpoints 100,700,-.1,.1 -nep_nev 5 -nep_tol 1e-9 -nep_target 140 -nep_nleigs_interpolation_degree 15 -nep_general -terse
196 : filter: sed -e "s/[+-]0\.0*i//g"
197 : requires: complex !single
198 : output_file: output/loaded_string_6.out
199 :
200 : test:
201 : suffix: 7
202 : args: -nep_type interpol -rg_type interval -rg_interval_endpoints 5,700 -nep_nev 5 -nep_target 100 -nep_interpol_interpolation_degree 20 -nep_ncv 20 -n 20 -nep_refine simple -nep_refine_its 1 -terse
203 : requires: !complex double
204 :
205 : test:
206 : suffix: 7_complex
207 : args: -nep_type interpol -rg_type interval -rg_interval_endpoints 5,700,-.1,.1 -nep_nev 5 -nep_target 100 -nep_interpol_interpolation_degree 20 -nep_ncv 20 -n 20 -nep_refine simple -nep_refine_its 1 -terse
208 : requires: complex double
209 : output_file: output/loaded_string_7.out
210 :
211 : testset:
212 : args: -nep_target 10 -nep_nev 3 -nep_tol 5e-10 -terse
213 : requires: !single
214 : output_file: output/loaded_string_8.out
215 : filter: sed -e "s/[+-]0\.0*i//g"
216 : test:
217 : suffix: 8
218 : args: -nep_type {{rii slp narnoldi}}
219 : test:
220 : suffix: 8_rii_thres
221 : args: -nep_type rii -nep_rii_deflation_threshold 5e-10
222 : test:
223 : suffix: 8_slp_thres
224 : args: -nep_type slp -nep_slp_deflation_threshold 5e-10
225 :
226 : test:
227 : suffix: 8_slp_two_thres
228 : args: -nep_type slp -nep_slp_deflation_threshold 5e-10 -nep_two_sided
229 :
230 : test:
231 : suffix: 9
232 : args: -nep_type ciss -rg_type ellipse -rg_ellipse_center 500 -rg_ellipse_radius 500 -rg_ellipse_vscale .1 -nep_ciss_moments 4 -nep_ciss_blocksize 5 -terse
233 : requires: complex double
234 :
235 : TEST*/
|
