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 24 : PetscErrorCode CheckNormalizedVectors(EPS eps)
29 : {
30 24 : PetscInt i,nconv;
31 24 : Mat A;
32 24 : Vec xr,xi;
33 24 : PetscReal error=0.0,normr;
34 : #if !defined(PETSC_USE_COMPLEX)
35 24 : PetscReal normi;
36 : #endif
37 :
38 24 : PetscFunctionBeginUser;
39 24 : PetscCall(EPSGetConverged(eps,&nconv));
40 24 : if (nconv>0) {
41 24 : PetscCall(EPSGetOperators(eps,&A,NULL));
42 24 : PetscCall(MatCreateVecs(A,&xr,&xi));
43 239 : for (i=0;i<nconv;i++) {
44 215 : PetscCall(EPSGetEigenvector(eps,i,xr,xi));
45 : #if defined(PETSC_USE_COMPLEX)
46 : PetscCall(VecNorm(xr,NORM_2,&normr));
47 : error = PetscMax(error,PetscAbsReal(normr-PetscRealConstant(1.0)));
48 : #else
49 215 : PetscCall(VecNormBegin(xr,NORM_2,&normr));
50 215 : PetscCall(VecNormBegin(xi,NORM_2,&normi));
51 215 : PetscCall(VecNormEnd(xr,NORM_2,&normr));
52 215 : PetscCall(VecNormEnd(xi,NORM_2,&normi));
53 215 : error = PetscMax(error,PetscAbsReal(SlepcAbsEigenvalue(normr,normi)-PetscRealConstant(1.0)));
54 : #endif
55 : }
56 24 : PetscCall(VecDestroy(&xr));
57 24 : PetscCall(VecDestroy(&xi));
58 24 : if (error>100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Vectors are not normalized. Error=%g\n",(double)error));
59 : }
60 24 : PetscFunctionReturn(PETSC_SUCCESS);
61 : }
62 :
63 28 : int main(int argc,char **argv)
64 : {
65 28 : Vec v0; /* initial vector */
66 28 : Mat A; /* operator matrix */
67 28 : EPS eps; /* eigenproblem solver context */
68 28 : PetscReal tol=0.5*PETSC_SMALL;
69 28 : PetscInt N,m=15,nev;
70 28 : PetscScalar origin=0.0;
71 28 : PetscBool flg,delay,skipnorm=PETSC_FALSE;
72 :
73 28 : PetscFunctionBeginUser;
74 28 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
75 :
76 28 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
77 28 : N = m*(m+1)/2;
78 28 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",N,m));
79 28 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-skipnorm",&skipnorm,NULL));
80 :
81 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82 : Compute the operator matrix that defines the eigensystem, Ax=kx
83 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84 :
85 28 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
86 28 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
87 28 : PetscCall(MatSetFromOptions(A));
88 28 : PetscCall(MatMarkovModel(m,A));
89 :
90 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91 : Create the eigensolver and set various options
92 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93 :
94 : /*
95 : Create eigensolver context
96 : */
97 28 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
98 :
99 : /*
100 : Set operators. In this case, it is a standard eigenvalue problem
101 : */
102 28 : PetscCall(EPSSetOperators(eps,A,NULL));
103 28 : PetscCall(EPSSetProblemType(eps,EPS_NHEP));
104 28 : PetscCall(EPSSetTolerances(eps,tol,PETSC_CURRENT));
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 28 : PetscCall(EPSSetEigenvalueComparison(eps,MyEigenSort,&origin));
111 :
112 : /*
113 : Set solver parameters at runtime
114 : */
115 28 : PetscCall(EPSSetFromOptions(eps));
116 28 : PetscCall(PetscObjectTypeCompare((PetscObject)eps,EPSARNOLDI,&flg));
117 28 : 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 28 : PetscCall(MatCreateVecs(A,&v0,NULL));
127 28 : PetscCall(VecSetValue(v0,0,-1.5,INSERT_VALUES));
128 28 : PetscCall(VecSetValue(v0,1,2.1,INSERT_VALUES));
129 28 : PetscCall(VecAssemblyBegin(v0));
130 28 : PetscCall(VecAssemblyEnd(v0));
131 28 : PetscCall(EPSSetInitialSpace(eps,1,&v0));
132 :
133 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134 : Solve the eigensystem
135 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136 :
137 28 : PetscCall(EPSSolve(eps));
138 28 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
139 28 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
140 :
141 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142 : Display solution and clean up
143 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144 :
145 28 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
146 28 : if (!skipnorm) PetscCall(CheckNormalizedVectors(eps));
147 28 : PetscCall(EPSDestroy(&eps));
148 28 : PetscCall(MatDestroy(&A));
149 28 : PetscCall(VecDestroy(&v0));
150 28 : PetscCall(SlepcFinalize());
151 : return 0;
152 : }
153 :
154 28 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
155 : {
156 28 : const PetscReal cst = 0.5/(PetscReal)(m-1);
157 28 : PetscReal pd,pu;
158 28 : PetscInt Istart,Iend,i,j,jmax,ix=0;
159 :
160 28 : PetscFunctionBeginUser;
161 28 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
162 448 : for (i=1;i<=m;i++) {
163 420 : jmax = m-i+1;
164 3780 : for (j=1;j<=jmax;j++) {
165 3360 : ix = ix + 1;
166 3360 : if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
167 3120 : if (j!=jmax) {
168 2730 : pd = cst*(PetscReal)(i+j-1);
169 : /* north */
170 2730 : if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
171 2366 : else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
172 : /* east */
173 2730 : if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
174 2366 : else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
175 : }
176 : /* south */
177 3120 : pu = 0.5 - cst*(PetscReal)(i+j-3);
178 3120 : if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
179 : /* west */
180 3360 : if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
181 : }
182 : }
183 28 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
184 28 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
185 28 : 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 66909 : PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx)
196 : {
197 66909 : PetscScalar origin = *(PetscScalar*)ctx;
198 66909 : PetscReal d;
199 :
200 66909 : PetscFunctionBeginUser;
201 66909 : d = (SlepcAbsEigenvalue(br-origin,bi) - SlepcAbsEigenvalue(ar-origin,ai))/PetscMax(SlepcAbsEigenvalue(ar-origin,ai),SlepcAbsEigenvalue(br-origin,bi));
202 66909 : *r = d > PETSC_SQRT_MACHINE_EPSILON ? 1 : (d < -PETSC_SQRT_MACHINE_EPSILON ? -1 : PetscSign(PetscRealPart(br)));
203 66909 : PetscFunctionReturn(PETSC_SUCCESS);
204 : }
205 :
206 : /*TEST
207 :
208 : testset:
209 : args: -eps_nev 4
210 : output_file: output/test9_1.out
211 : filter: sed -e "s/97136/97137/g"
212 : test:
213 : suffix: 1
214 : args: -eps_type {{krylovschur arnoldi lapack}} -eps_ncv 8 -eps_max_it 300
215 : test:
216 : suffix: 1_gd
217 : args: -eps_type gd -st_pc_type none
218 : test:
219 : suffix: 1_gd2
220 : args: -eps_type gd -eps_gd_double_expansion -st_pc_type none
221 :
222 : test:
223 : suffix: 2
224 : args: -eps_balance {{none oneside twoside}} -eps_krylovschur_locking {{0 1}} -eps_nev 4 -eps_max_it 1500
225 : requires: double
226 : output_file: output/test9_1.out
227 :
228 : test:
229 : suffix: 3
230 : nsize: 2
231 : args: -eps_type arnoldi -eps_arnoldi_delayed -eps_largest_real -eps_nev 3 -eps_tol 1e-7 -bv_orthog_refine {{never ifneeded}} -skipnorm
232 : requires: !single
233 : output_file: output/test9_3.out
234 :
235 : test:
236 : suffix: 4
237 : args: -eps_nev 4 -eps_true_residual
238 : requires: !single
239 : output_file: output/test9_1.out
240 :
241 : test:
242 : suffix: 5
243 : args: -eps_type jd -eps_nev 3 -eps_target .5 -eps_harmonic -st_ksp_type bicg -st_pc_type lu -eps_jd_minv 2
244 : filter: sed -e "s/[+-]0\.0*i//g"
245 : requires: !single
246 :
247 : test:
248 : suffix: 5_arpack
249 : args: -eps_nev 3 -st_type sinvert -eps_target .5 -eps_type arpack -eps_ncv 10
250 : requires: arpack !single
251 : output_file: output/test9_5.out
252 :
253 : testset:
254 : 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
255 : requires: !single
256 : output_file: output/test9_6.out
257 : test:
258 : suffix: 6
259 : test:
260 : suffix: 6_hankel
261 : args: -eps_ciss_extraction hankel -eps_ciss_spurious_threshold 1e-6 -eps_ncv 64
262 : test:
263 : suffix: 6_cheby
264 : args: -eps_ciss_quadrule chebyshev
265 : test:
266 : suffix: 6_hankel_cheby
267 : args: -eps_ciss_extraction hankel -eps_ciss_quadrule chebyshev -eps_ncv 64
268 : test:
269 : suffix: 6_refine
270 : args: -eps_ciss_moments 4 -eps_ciss_blocksize 5 -eps_ciss_refine_inner 1 -eps_ciss_refine_blocksize 2
271 : test:
272 : suffix: 6_bcgs
273 : args: -eps_ciss_realmats -eps_ciss_ksp_type bcgs -eps_ciss_pc_type ilu -eps_ciss_integration_points 8
274 :
275 : test:
276 : suffix: 6_cheby_interval
277 : 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
278 : requires: !single
279 : output_file: output/test9_6.out
280 :
281 : testset:
282 : args: -eps_nev 4 -eps_two_sided -eps_view_vectors ::ascii_info -eps_view_values
283 : filter: sed -e "s/\(0x[0-9a-fA-F]*\)/objectid/"
284 : test:
285 : suffix: 7_real
286 : requires: !single !complex
287 : test:
288 : suffix: 7
289 : requires: !single complex
290 :
291 : test:
292 : suffix: 8
293 : args: -eps_nev 4 -eps_ncv 7 -eps_view_values draw -eps_monitor draw::draw_lg
294 : requires: x
295 : output_file: output/test9_1.out
296 :
297 : test:
298 : suffix: 5_ksphpddm
299 : args: -eps_nev 3 -st_type sinvert -eps_target .5 -st_ksp_type hpddm -st_ksp_hpddm_type gcrodr -eps_ncv 10
300 : requires: hpddm
301 : output_file: output/test9_5.out
302 :
303 : test:
304 : suffix: 5_pchpddm
305 : args: -eps_nev 3 -st_type sinvert -eps_target .5 -st_pc_type hpddm -st_pc_hpddm_coarse_pc_type lu -eps_ncv 10
306 : requires: hpddm
307 : output_file: output/test9_5.out
308 :
309 : TEST*/
|
