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 1 : int main(int argc,char **argv)
27 : {
28 1 : Mat A; /* operator matrix */
29 1 : EPS eps; /* eigenproblem solver context */
30 1 : EPSType type;
31 1 : PetscScalar target=0.5;
32 1 : PetscInt N,m=15,nev;
33 1 : PetscBool terse;
34 1 : PetscViewer viewer;
35 1 : char str[50];
36 :
37 1 : PetscFunctionBeginUser;
38 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
39 :
40 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
41 1 : N = m*(m+1)/2;
42 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n",N,m));
43 1 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-target",&target,NULL));
44 1 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),target,PETSC_FALSE));
45 1 : 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 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
52 1 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
53 1 : PetscCall(MatSetFromOptions(A));
54 1 : PetscCall(MatMarkovModel(m,A));
55 :
56 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57 : Create the eigensolver and set various options
58 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59 :
60 : /*
61 : Create eigensolver context
62 : */
63 1 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
64 :
65 : /*
66 : Set operators. In this case, it is a standard eigenvalue problem
67 : */
68 1 : PetscCall(EPSSetOperators(eps,A,NULL));
69 1 : PetscCall(EPSSetProblemType(eps,EPS_NHEP));
70 :
71 : /*
72 : Set the custom comparing routine in order to obtain the eigenvalues
73 : closest to the target on the right only
74 : */
75 1 : PetscCall(EPSSetEigenvalueComparison(eps,MyEigenSort,&target));
76 :
77 : /*
78 : Set solver parameters at runtime
79 : */
80 1 : PetscCall(EPSSetFromOptions(eps));
81 :
82 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83 : Solve the eigensystem
84 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85 :
86 1 : PetscCall(EPSSolve(eps));
87 :
88 : /*
89 : Optional: Get some information from the solver and display it
90 : */
91 1 : PetscCall(EPSGetType(eps,&type));
92 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
93 1 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
94 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
95 :
96 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97 : Display solution and clean up
98 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99 :
100 : /* show detailed info unless -terse option is given by user */
101 1 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
102 1 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
103 : else {
104 0 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
105 0 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
106 0 : PetscCall(EPSConvergedReasonView(eps,viewer));
107 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,viewer));
108 0 : PetscCall(PetscViewerPopFormat(viewer));
109 : }
110 1 : PetscCall(EPSDestroy(&eps));
111 1 : PetscCall(MatDestroy(&A));
112 1 : PetscCall(SlepcFinalize());
113 : return 0;
114 : }
115 :
116 : /*
117 : Matrix generator for a Markov model of a random walk on a triangular grid.
118 :
119 : This subroutine generates a test matrix that models a random walk on a
120 : triangular grid. This test example was used by G. W. Stewart ["{SRRIT} - a
121 : FORTRAN subroutine to calculate the dominant invariant subspaces of a real
122 : matrix", Tech. report. TR-514, University of Maryland (1978).] and in a few
123 : papers on eigenvalue problems by Y. Saad [see e.g. LAA, vol. 34, pp. 269-295
124 : (1980) ]. These matrices provide reasonably easy test problems for eigenvalue
125 : algorithms. The transpose of the matrix is stochastic and so it is known
126 : that one is an exact eigenvalue. One seeks the eigenvector of the transpose
127 : associated with the eigenvalue unity. The problem is to calculate the steady
128 : state probability distribution of the system, which is the eigevector
129 : associated with the eigenvalue one and scaled in such a way that the sum all
130 : the components is equal to one.
131 :
132 : Note: the code will actually compute the transpose of the stochastic matrix
133 : that contains the transition probabilities.
134 : */
135 1 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
136 : {
137 1 : const PetscReal cst = 0.5/(PetscReal)(m-1);
138 1 : PetscReal pd,pu;
139 1 : PetscInt Istart,Iend,i,j,jmax,ix=0;
140 :
141 1 : PetscFunctionBeginUser;
142 1 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
143 16 : for (i=1;i<=m;i++) {
144 15 : jmax = m-i+1;
145 135 : for (j=1;j<=jmax;j++) {
146 120 : ix = ix + 1;
147 120 : if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
148 120 : if (j!=jmax) {
149 105 : pd = cst*(PetscReal)(i+j-1);
150 : /* north */
151 105 : if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
152 91 : else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
153 : /* east */
154 105 : if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
155 91 : else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
156 : }
157 : /* south */
158 120 : pu = 0.5 - cst*(PetscReal)(i+j-3);
159 120 : if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
160 : /* west */
161 120 : if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
162 : }
163 : }
164 1 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
165 1 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
166 1 : PetscFunctionReturn(PETSC_SUCCESS);
167 : }
168 :
169 : /*
170 : Function for user-defined eigenvalue ordering criterion.
171 :
172 : Given two eigenvalues ar+i*ai and br+i*bi, the subroutine must choose
173 : one of them as the preferred one according to the criterion.
174 : In this example, the preferred value is the one closest to the target,
175 : but on the right side.
176 : */
177 5185 : PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx)
178 : {
179 5185 : PetscScalar target = *(PetscScalar*)ctx;
180 5185 : PetscReal da,db;
181 5185 : PetscBool aisright,bisright;
182 :
183 5185 : PetscFunctionBeginUser;
184 5185 : if (PetscRealPart(target) < PetscRealPart(ar)) aisright = PETSC_TRUE;
185 2950 : else aisright = PETSC_FALSE;
186 5185 : if (PetscRealPart(target) < PetscRealPart(br)) bisright = PETSC_TRUE;
187 2934 : else bisright = PETSC_FALSE;
188 5185 : if (aisright == bisright) {
189 : /* both are on the same side of the target */
190 4469 : da = SlepcAbsEigenvalue(ar-target,ai);
191 4469 : db = SlepcAbsEigenvalue(br-target,bi);
192 4469 : if (da < db) *r = -1;
193 3872 : else if (da > db) *r = 1;
194 0 : else *r = 0;
195 716 : } else if (aisright && !bisright) *r = -1; /* 'a' is on the right */
196 366 : else *r = 1; /* 'b' is on the right */
197 5185 : PetscFunctionReturn(PETSC_SUCCESS);
198 : }
199 :
200 : /*TEST
201 :
202 : test:
203 : suffix: 1
204 : args: -eps_nev 4 -terse
205 : requires: !single
206 : filter: sed -e "s/[+-]0\.0*i//g"
207 :
208 : TEST*/
|
