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[] = "Solves the same problem as in ex5, but with a user-defined sorting criterion."
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 a custom spectrum selection.\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 :
23 : PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx);
24 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A);
25 :
26 5 : int main(int argc,char **argv)
27 : {
28 5 : Vec v0; /* initial vector */
29 5 : Mat A; /* operator matrix */
30 5 : EPS eps; /* eigenproblem solver context */
31 5 : ST st; /* spectral transformation associated */
32 5 : PetscReal tol=PETSC_SMALL;
33 5 : PetscScalar target=0.5;
34 5 : PetscInt N,m=15,nev;
35 5 : char str[50];
36 :
37 5 : PetscFunctionBeginUser;
38 5 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
39 :
40 5 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
41 5 : N = m*(m+1)/2;
42 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n",N,m));
43 5 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-target",&target,NULL));
44 5 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),target,PETSC_FALSE));
45 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Searching closest eigenvalues to the right of %s.\n\n",str));
46 :
47 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48 : Compute the operator matrix that defines the eigensystem, Ax=kx
49 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50 :
51 5 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
52 5 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
53 5 : PetscCall(MatSetFromOptions(A));
54 5 : PetscCall(MatMarkovModel(m,A));
55 :
56 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57 : Create the eigensolver and set various options
58 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59 :
60 : /*
61 : Create eigensolver context
62 : */
63 5 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
64 :
65 : /*
66 : Set operators. In this case, it is a standard eigenvalue problem
67 : */
68 5 : PetscCall(EPSSetOperators(eps,A,NULL));
69 5 : PetscCall(EPSSetProblemType(eps,EPS_NHEP));
70 5 : PetscCall(EPSSetTolerances(eps,tol,PETSC_CURRENT));
71 :
72 : /*
73 : Set the custom comparing routine in order to obtain the eigenvalues
74 : closest to the target on the right only
75 : */
76 5 : PetscCall(EPSSetEigenvalueComparison(eps,MyEigenSort,&target));
77 :
78 : /*
79 : Set solver parameters at runtime
80 : */
81 5 : PetscCall(EPSSetFromOptions(eps));
82 :
83 : /*
84 : Set the preconditioner based on A - target * I
85 : */
86 5 : PetscCall(EPSGetST(eps,&st));
87 5 : PetscCall(STSetShift(st,target));
88 :
89 : /*
90 : Set the initial vector. This is optional, if not done the initial
91 : vector is set to random values
92 : */
93 5 : PetscCall(MatCreateVecs(A,&v0,NULL));
94 5 : PetscCall(VecSet(v0,1.0));
95 5 : PetscCall(EPSSetInitialSpace(eps,1,&v0));
96 :
97 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98 : Solve the eigensystem
99 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100 :
101 5 : PetscCall(EPSSolve(eps));
102 5 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
103 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
104 :
105 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106 : Display solution and clean up
107 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108 :
109 5 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
110 5 : PetscCall(EPSDestroy(&eps));
111 5 : PetscCall(MatDestroy(&A));
112 5 : PetscCall(VecDestroy(&v0));
113 5 : PetscCall(SlepcFinalize());
114 : return 0;
115 : }
116 :
117 5 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
118 : {
119 5 : const PetscReal cst = 0.5/(PetscReal)(m-1);
120 5 : PetscReal pd,pu;
121 5 : PetscInt Istart,Iend,i,j,jmax,ix=0;
122 :
123 5 : PetscFunctionBeginUser;
124 5 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
125 80 : for (i=1;i<=m;i++) {
126 75 : jmax = m-i+1;
127 675 : for (j=1;j<=jmax;j++) {
128 600 : ix = ix + 1;
129 600 : if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
130 600 : if (j!=jmax) {
131 525 : pd = cst*(PetscReal)(i+j-1);
132 : /* north */
133 525 : if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
134 455 : else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
135 : /* east */
136 525 : if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
137 455 : else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
138 : }
139 : /* south */
140 600 : pu = 0.5 - cst*(PetscReal)(i+j-3);
141 600 : if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
142 : /* west */
143 600 : if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
144 : }
145 : }
146 5 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
147 5 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
148 5 : PetscFunctionReturn(PETSC_SUCCESS);
149 : }
150 :
151 : /*
152 : Function for user-defined eigenvalue ordering criterion.
153 :
154 : Given two eigenvalues ar+i*ai and br+i*bi, the subroutine must choose
155 : one of them as the preferred one according to the criterion.
156 : In this example, the preferred value is the one closest to the target,
157 : but on the right side.
158 : */
159 18909 : PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx)
160 : {
161 18909 : PetscScalar target = *(PetscScalar*)ctx;
162 18909 : PetscReal da,db;
163 18909 : PetscBool aisright,bisright;
164 :
165 18909 : PetscFunctionBeginUser;
166 18909 : if (PetscRealPart(target) < PetscRealPart(ar)) aisright = PETSC_TRUE;
167 11217 : else aisright = PETSC_FALSE;
168 18909 : if (PetscRealPart(target) < PetscRealPart(br)) bisright = PETSC_TRUE;
169 16225 : else bisright = PETSC_FALSE;
170 18909 : if (aisright == bisright) {
171 : /* both are on the same side of the target */
172 13063 : da = SlepcAbsEigenvalue(ar-target,ai);
173 13063 : db = SlepcAbsEigenvalue(br-target,bi);
174 13063 : if (da < db) *r = -1;
175 3591 : else if (da > db) *r = 1;
176 0 : else *r = 0;
177 5846 : } else if (aisright && !bisright) *r = -1; /* 'a' is on the right */
178 419 : else *r = 1; /* 'b' is on the right */
179 18909 : PetscFunctionReturn(PETSC_SUCCESS);
180 : }
181 :
182 : /*TEST
183 :
184 : testset:
185 : args: -eps_nev 4
186 : requires: !single
187 : output_file: output/test11_1.out
188 : test:
189 : suffix: 1
190 : args: -eps_type {{krylovschur arnoldi lapack}} -st_type sinvert
191 : test:
192 : suffix: 1_ks_cayley
193 : args: -st_type cayley -st_cayley_antishift 1
194 :
195 : test:
196 : suffix: 2
197 : args: -target 0.77 -eps_type gd -eps_nev 4 -eps_tol 1e-7 -eps_gd_krylov_start -eps_gd_blocksize 3
198 : requires: double
199 : filter: sed -e "s/[+-]0\.0*i//g"
200 :
201 : TEST*/
|
