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[] = "Eigenvalue problem associated with a Markov model of a random walk on a triangular grid. "
12 : "It is a standard nonsymmetric eigenproblem with real eigenvalues and the rightmost eigenvalue is known to be 1.\n"
13 : "This example illustrates how the user can set the initial vector.\n\n"
14 : "The command line options are:\n"
15 : " -m <m>, where <m> = number of grid subdivisions in each dimension.\n\n";
16 :
17 : #include <slepceps.h>
18 :
19 : /*
20 : User-defined routines
21 : */
22 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A);
23 : PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx);
24 :
25 : /*
26 : Check if computed eigenvectors have unit norm
27 : */
28 25 : PetscErrorCode CheckNormalizedVectors(EPS eps)
29 : {
30 25 : PetscInt i,nconv;
31 25 : Mat A;
32 25 : Vec xr,xi;
33 25 : PetscReal error=0.0,normr;
34 : #if !defined(PETSC_USE_COMPLEX)
35 : PetscReal normi;
36 : #endif
37 :
38 25 : PetscFunctionBeginUser;
39 25 : PetscCall(EPSGetConverged(eps,&nconv));
40 25 : if (nconv>0) {
41 25 : PetscCall(EPSGetOperators(eps,&A,NULL));
42 25 : PetscCall(MatCreateVecs(A,&xr,&xi));
43 244 : for (i=0;i<nconv;i++) {
44 219 : PetscCall(EPSGetEigenvector(eps,i,xr,xi));
45 : #if defined(PETSC_USE_COMPLEX)
46 219 : PetscCall(VecNorm(xr,NORM_2,&normr));
47 262 : error = PetscMax(error,PetscAbsReal(normr-PetscRealConstant(1.0)));
48 : #else
49 : PetscCall(VecNormBegin(xr,NORM_2,&normr));
50 : PetscCall(VecNormBegin(xi,NORM_2,&normi));
51 : PetscCall(VecNormEnd(xr,NORM_2,&normr));
52 : PetscCall(VecNormEnd(xi,NORM_2,&normi));
53 : error = PetscMax(error,PetscAbsReal(SlepcAbsEigenvalue(normr,normi)-PetscRealConstant(1.0)));
54 : #endif
55 : }
56 25 : PetscCall(VecDestroy(&xr));
57 25 : PetscCall(VecDestroy(&xi));
58 25 : if (error>100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Vectors are not normalized. Error=%g\n",(double)error));
59 : }
60 25 : PetscFunctionReturn(PETSC_SUCCESS);
61 : }
62 :
63 29 : int main(int argc,char **argv)
64 : {
65 29 : Vec v0; /* initial vector */
66 29 : Mat A; /* operator matrix */
67 29 : EPS eps; /* eigenproblem solver context */
68 29 : PetscReal tol=0.5*PETSC_SMALL;
69 29 : PetscInt N,m=15,nev;
70 29 : PetscScalar origin=0.0;
71 29 : PetscBool flg,delay,skipnorm=PETSC_FALSE;
72 :
73 29 : PetscFunctionBeginUser;
74 29 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
75 :
76 29 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
77 29 : N = m*(m+1)/2;
78 29 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",N,m));
79 29 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-skipnorm",&skipnorm,NULL));
80 :
81 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82 : Compute the operator matrix that defines the eigensystem, Ax=kx
83 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84 :
85 29 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
86 29 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
87 29 : PetscCall(MatSetFromOptions(A));
88 29 : PetscCall(MatMarkovModel(m,A));
89 :
90 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91 : Create the eigensolver and set various options
92 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93 :
94 : /*
95 : Create eigensolver context
96 : */
97 29 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
98 :
99 : /*
100 : Set operators. In this case, it is a standard eigenvalue problem
101 : */
102 29 : PetscCall(EPSSetOperators(eps,A,NULL));
103 29 : PetscCall(EPSSetProblemType(eps,EPS_NHEP));
104 29 : PetscCall(EPSSetTolerances(eps,tol,PETSC_DEFAULT));
105 :
106 : /*
107 : Set the custom comparing routine in order to obtain the eigenvalues
108 : closest to the target on the right only
109 : */
110 29 : PetscCall(EPSSetEigenvalueComparison(eps,MyEigenSort,&origin));
111 :
112 : /*
113 : Set solver parameters at runtime
114 : */
115 29 : PetscCall(EPSSetFromOptions(eps));
116 29 : PetscCall(PetscObjectTypeCompare((PetscObject)eps,EPSARNOLDI,&flg));
117 29 : if (flg) {
118 5 : PetscCall(EPSArnoldiGetDelayed(eps,&delay));
119 5 : if (delay) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Warning: delayed reorthogonalization may be unstable\n"));
120 : }
121 :
122 : /*
123 : Set the initial vector. This is optional, if not done the initial
124 : vector is set to random values
125 : */
126 29 : PetscCall(MatCreateVecs(A,&v0,NULL));
127 29 : PetscCall(VecSetValue(v0,0,-1.5,INSERT_VALUES));
128 29 : PetscCall(VecSetValue(v0,1,2.1,INSERT_VALUES));
129 29 : PetscCall(VecAssemblyBegin(v0));
130 29 : PetscCall(VecAssemblyEnd(v0));
131 29 : PetscCall(EPSSetInitialSpace(eps,1,&v0));
132 :
133 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134 : Solve the eigensystem
135 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136 :
137 29 : PetscCall(EPSSolve(eps));
138 29 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
139 29 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
140 :
141 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142 : Display solution and clean up
143 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144 :
145 29 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
146 29 : if (!skipnorm) PetscCall(CheckNormalizedVectors(eps));
147 29 : PetscCall(EPSDestroy(&eps));
148 29 : PetscCall(MatDestroy(&A));
149 29 : PetscCall(VecDestroy(&v0));
150 29 : PetscCall(SlepcFinalize());
151 : return 0;
152 : }
153 :
154 29 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
155 : {
156 29 : const PetscReal cst = 0.5/(PetscReal)(m-1);
157 29 : PetscReal pd,pu;
158 29 : PetscInt Istart,Iend,i,j,jmax,ix=0;
159 :
160 29 : PetscFunctionBeginUser;
161 29 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
162 464 : for (i=1;i<=m;i++) {
163 435 : jmax = m-i+1;
164 3915 : for (j=1;j<=jmax;j++) {
165 3480 : ix = ix + 1;
166 3480 : if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
167 3240 : if (j!=jmax) {
168 2835 : pd = cst*(PetscReal)(i+j-1);
169 : /* north */
170 2835 : if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
171 2457 : else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
172 : /* east */
173 2835 : if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
174 2457 : else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
175 : }
176 : /* south */
177 3240 : pu = 0.5 - cst*(PetscReal)(i+j-3);
178 3240 : if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
179 : /* west */
180 3480 : if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
181 : }
182 : }
183 29 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
184 29 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
185 29 : PetscFunctionReturn(PETSC_SUCCESS);
186 : }
187 :
188 : /*
189 : Function for user-defined eigenvalue ordering criterion.
190 :
191 : Given two eigenvalues ar+i*ai and br+i*bi, the subroutine must choose
192 : one of them as the preferred one according to the criterion.
193 : In this example, the preferred value is the one furthest away from the origin.
194 : */
195 83485 : PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx)
196 : {
197 83485 : PetscScalar origin = *(PetscScalar*)ctx;
198 83485 : PetscReal d;
199 :
200 83485 : PetscFunctionBeginUser;
201 83485 : d = (SlepcAbsEigenvalue(br-origin,bi) - SlepcAbsEigenvalue(ar-origin,ai))/PetscMax(SlepcAbsEigenvalue(ar-origin,ai),SlepcAbsEigenvalue(br-origin,bi));
202 83485 : *r = d > PETSC_SQRT_MACHINE_EPSILON ? 1 : (d < -PETSC_SQRT_MACHINE_EPSILON ? -1 : PetscSign(PetscRealPart(br)));
203 83485 : PetscFunctionReturn(PETSC_SUCCESS);
204 : }
205 :
206 : /*TEST
207 :
208 : testset:
209 : args: -eps_nev 4
210 : output_file: output/test9_1.out
211 : test:
212 : suffix: 1
213 : args: -eps_type {{krylovschur arnoldi lapack}} -eps_ncv 7 -eps_max_it 300
214 : test:
215 : suffix: 1_gd
216 : args: -eps_type gd -st_pc_type none
217 : test:
218 : suffix: 1_gd2
219 : args: -eps_type gd -eps_gd_double_expansion -st_pc_type none
220 :
221 : test:
222 : suffix: 2
223 : args: -eps_balance {{none oneside twoside}} -eps_krylovschur_locking {{0 1}} -eps_nev 4 -eps_max_it 1500
224 : requires: double
225 : output_file: output/test9_1.out
226 :
227 : test:
228 : suffix: 3
229 : nsize: 2
230 : args: -eps_type arnoldi -eps_arnoldi_delayed -eps_largest_real -eps_nev 3 -eps_tol 1e-7 -bv_orthog_refine {{never ifneeded}} -skipnorm
231 : requires: !single
232 : output_file: output/test9_3.out
233 :
234 : test:
235 : suffix: 4
236 : args: -eps_nev 4 -eps_true_residual
237 : requires: !single
238 : output_file: output/test9_1.out
239 :
240 : test:
241 : suffix: 5
242 : args: -eps_type jd -eps_nev 3 -eps_target .5 -eps_harmonic -st_ksp_type bicg -st_pc_type lu -eps_jd_minv 2
243 : filter: sed -e "s/[+-]0\.0*i//g"
244 : requires: !single
245 :
246 : test:
247 : suffix: 5_arpack
248 : args: -eps_nev 3 -st_type sinvert -eps_target .5 -eps_type arpack -eps_ncv 10
249 : requires: arpack !single
250 : output_file: output/test9_5.out
251 :
252 : testset:
253 : args: -eps_type ciss -eps_tol 1e-9 -rg_type ellipse -rg_ellipse_center 0.55 -rg_ellipse_radius 0.05 -rg_ellipse_vscale 0.1 -eps_ciss_usest 0 -eps_all
254 : requires: !single
255 : output_file: output/test9_6.out
256 : test:
257 : suffix: 6
258 : test:
259 : suffix: 6_hankel
260 : args: -eps_ciss_extraction hankel -eps_ciss_spurious_threshold 1e-6 -eps_ncv 64
261 : test:
262 : suffix: 6_cheby
263 : args: -eps_ciss_quadrule chebyshev
264 : test:
265 : suffix: 6_hankel_cheby
266 : args: -eps_ciss_extraction hankel -eps_ciss_quadrule chebyshev -eps_ncv 64
267 : test:
268 : suffix: 6_refine
269 : args: -eps_ciss_moments 4 -eps_ciss_blocksize 5 -eps_ciss_refine_inner 1 -eps_ciss_refine_blocksize 2
270 : test:
271 : suffix: 6_bcgs
272 : args: -eps_ciss_realmats -eps_ciss_ksp_type bcgs -eps_ciss_pc_type ilu -eps_ciss_integration_points 8
273 :
274 : test:
275 : suffix: 6_cheby_interval
276 : args: -eps_type ciss -eps_tol 1e-9 -rg_type interval -rg_interval_endpoints 0.5,0.6 -eps_ciss_quadrule chebyshev -eps_ciss_usest 0 -eps_all
277 : requires: !single
278 : output_file: output/test9_6.out
279 :
280 : testset:
281 : args: -eps_nev 4 -eps_two_sided -eps_view_vectors ::ascii_info -eps_view_values
282 : filter: sed -e "s/\(0x[0-9a-fA-F]*\)/objectid/"
283 : test:
284 : suffix: 7_real
285 : requires: !single !complex
286 : test:
287 : suffix: 7
288 : requires: !single complex
289 :
290 : test:
291 : suffix: 8
292 : args: -eps_nev 4 -eps_ncv 7 -eps_view_values draw -eps_monitor draw::draw_lg
293 : requires: x
294 : output_file: output/test9_1.out
295 :
296 : test:
297 : suffix: 5_ksphpddm
298 : args: -eps_nev 3 -st_type sinvert -eps_target .5 -st_ksp_type hpddm -st_ksp_hpddm_type gcrodr -eps_ncv 10
299 : requires: hpddm
300 : output_file: output/test9_5.out
301 :
302 : test:
303 : suffix: 5_pchpddm
304 : args: -eps_nev 3 -st_type sinvert -eps_target .5 -st_pc_type hpddm -st_pc_hpddm_coarse_pc_type lu -eps_ncv 10
305 : requires: hpddm
306 : output_file: output/test9_5.out
307 :
308 : TEST*/
|
