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[] = "Illustrates the computation of left eigenvectors.\n\n"
12 : "The problem is the Markov model as in ex5.c.\n"
13 : "The command line options are:\n"
14 : " -m <m>, where <m> = number of grid subdivisions in each dimension.\n\n";
15 :
16 : #include <slepceps.h>
17 :
18 : /*
19 : User-defined routines
20 : */
21 : PetscErrorCode MatMarkovModel(PetscInt,Mat);
22 : PetscErrorCode ComputeResidualNorm(Mat,PetscBool,PetscScalar,PetscScalar,Vec,Vec,Vec,PetscReal*);
23 :
24 4 : int main(int argc,char **argv)
25 : {
26 4 : Vec v0,w0; /* initial vectors */
27 4 : Mat A; /* operator matrix */
28 4 : EPS eps; /* eigenproblem solver context */
29 4 : EPSType type;
30 4 : PetscInt i,N,m=15,nconv;
31 4 : PetscBool twosided;
32 4 : PetscReal nrmr,nrml=0.0,re,im,lev;
33 4 : PetscScalar *kr,*ki;
34 4 : Vec t,*xr,*xi,*yr,*yi;
35 4 : PetscMPIInt rank;
36 :
37 4 : PetscFunctionBeginUser;
38 4 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
39 :
40 4 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
41 4 : N = m*(m+1)/2;
42 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",N,m));
43 :
44 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45 : Compute the operator matrix that defines the eigensystem, Ax=kx
46 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
47 :
48 4 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
49 4 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
50 4 : PetscCall(MatSetFromOptions(A));
51 4 : PetscCall(MatMarkovModel(m,A));
52 4 : PetscCall(MatCreateVecs(A,NULL,&t));
53 :
54 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55 : Create the eigensolver and set various options
56 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57 :
58 4 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
59 4 : PetscCall(EPSSetOperators(eps,A,NULL));
60 4 : PetscCall(EPSSetProblemType(eps,EPS_NHEP));
61 :
62 : /* use a two-sided algorithm to compute left eigenvectors as well */
63 4 : PetscCall(EPSSetTwoSided(eps,PETSC_TRUE));
64 :
65 : /* allow user to change settings at run time */
66 4 : PetscCall(EPSSetFromOptions(eps));
67 4 : PetscCall(EPSGetTwoSided(eps,&twosided));
68 :
69 : /*
70 : Set the initial vectors. This is optional, if not done the initial
71 : vectors are set to random values
72 : */
73 4 : PetscCall(MatCreateVecs(A,&v0,&w0));
74 4 : PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
75 4 : if (!rank) {
76 4 : PetscCall(VecSetValue(v0,0,1.0,INSERT_VALUES));
77 4 : PetscCall(VecSetValue(v0,1,1.0,INSERT_VALUES));
78 4 : PetscCall(VecSetValue(v0,2,1.0,INSERT_VALUES));
79 4 : PetscCall(VecSetValue(w0,0,2.0,INSERT_VALUES));
80 4 : PetscCall(VecSetValue(w0,2,0.5,INSERT_VALUES));
81 : }
82 4 : PetscCall(VecAssemblyBegin(v0));
83 4 : PetscCall(VecAssemblyBegin(w0));
84 4 : PetscCall(VecAssemblyEnd(v0));
85 4 : PetscCall(VecAssemblyEnd(w0));
86 4 : PetscCall(EPSSetInitialSpace(eps,1,&v0));
87 4 : PetscCall(EPSSetLeftInitialSpace(eps,1,&w0));
88 :
89 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90 : Solve the eigensystem
91 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92 :
93 4 : PetscCall(EPSSolve(eps));
94 :
95 : /*
96 : Optional: Get some information from the solver and display it
97 : */
98 4 : PetscCall(EPSGetType(eps,&type));
99 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
100 :
101 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102 : Display solution and clean up
103 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104 :
105 : /*
106 : Get number of converged approximate eigenpairs
107 : */
108 4 : PetscCall(EPSGetConverged(eps,&nconv));
109 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %" PetscInt_FMT "\n\n",nconv));
110 4 : PetscCall(PetscMalloc2(nconv,&kr,nconv,&ki));
111 4 : PetscCall(VecDuplicateVecs(t,nconv,&xr));
112 4 : PetscCall(VecDuplicateVecs(t,nconv,&xi));
113 4 : if (twosided) {
114 4 : PetscCall(VecDuplicateVecs(t,nconv,&yr));
115 4 : PetscCall(VecDuplicateVecs(t,nconv,&yi));
116 : }
117 :
118 4 : if (nconv>0) {
119 : /*
120 : Display eigenvalues and relative errors
121 : */
122 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,
123 : " k ||Ax-kx|| ||y'A-y'k||\n"
124 : " ---------------- ------------------ ------------------\n"));
125 :
126 21 : for (i=0;i<nconv;i++) {
127 : /*
128 : Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
129 : ki (imaginary part)
130 : */
131 17 : PetscCall(EPSGetEigenpair(eps,i,&kr[i],&ki[i],xr[i],xi[i]));
132 17 : if (twosided) PetscCall(EPSGetLeftEigenvector(eps,i,yr[i],yi[i]));
133 : /*
134 : Compute the residual norms associated to each eigenpair
135 : */
136 17 : PetscCall(ComputeResidualNorm(A,PETSC_FALSE,kr[i],ki[i],xr[i],xi[i],t,&nrmr));
137 17 : if (twosided) PetscCall(ComputeResidualNorm(A,PETSC_TRUE,kr[i],ki[i],yr[i],yi[i],t,&nrml));
138 :
139 : #if defined(PETSC_USE_COMPLEX)
140 17 : re = PetscRealPart(kr[i]);
141 17 : im = PetscImaginaryPart(kr[i]);
142 : #else
143 : re = kr[i];
144 : im = ki[i];
145 : #endif
146 17 : if (im!=0.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD," %8f%+8fi %12g %12g\n",(double)re,(double)im,(double)nrmr,(double)nrml));
147 17 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD," %12f %12g %12g\n",(double)re,(double)nrmr,(double)nrml));
148 : }
149 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
150 : /*
151 : Check bi-orthogonality of eigenvectors
152 : */
153 4 : if (twosided) {
154 4 : PetscCall(VecCheckOrthogonality(xr,nconv,yr,nconv,NULL,NULL,&lev));
155 4 : if (lev<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors < 100*eps\n\n"));
156 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors: %g\n\n",(double)lev));
157 : }
158 : }
159 :
160 4 : PetscCall(EPSDestroy(&eps));
161 4 : PetscCall(MatDestroy(&A));
162 4 : PetscCall(VecDestroy(&v0));
163 4 : PetscCall(VecDestroy(&w0));
164 4 : PetscCall(VecDestroy(&t));
165 4 : PetscCall(PetscFree2(kr,ki));
166 4 : PetscCall(VecDestroyVecs(nconv,&xr));
167 4 : PetscCall(VecDestroyVecs(nconv,&xi));
168 4 : if (twosided) {
169 4 : PetscCall(VecDestroyVecs(nconv,&yr));
170 4 : PetscCall(VecDestroyVecs(nconv,&yi));
171 : }
172 4 : PetscCall(SlepcFinalize());
173 : return 0;
174 : }
175 :
176 : /*
177 : Matrix generator for a Markov model of a random walk on a triangular grid.
178 :
179 : This subroutine generates a test matrix that models a random walk on a
180 : triangular grid. This test example was used by G. W. Stewart ["{SRRIT} - a
181 : FORTRAN subroutine to calculate the dominant invariant subspaces of a real
182 : matrix", Tech. report. TR-514, University of Maryland (1978).] and in a few
183 : papers on eigenvalue problems by Y. Saad [see e.g. LAA, vol. 34, pp. 269-295
184 : (1980) ]. These matrices provide reasonably easy test problems for eigenvalue
185 : algorithms. The transpose of the matrix is stochastic and so it is known
186 : that one is an exact eigenvalue. One seeks the eigenvector of the transpose
187 : associated with the eigenvalue unity. The problem is to calculate the steady
188 : state probability distribution of the system, which is the eigevector
189 : associated with the eigenvalue one and scaled in such a way that the sum all
190 : the components is equal to one.
191 :
192 : Note: the code will actually compute the transpose of the stochastic matrix
193 : that contains the transition probabilities.
194 : */
195 4 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
196 : {
197 4 : const PetscReal cst = 0.5/(PetscReal)(m-1);
198 4 : PetscReal pd,pu;
199 4 : PetscInt Istart,Iend,i,j,jmax,ix=0;
200 :
201 4 : PetscFunctionBeginUser;
202 4 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
203 64 : for (i=1;i<=m;i++) {
204 60 : jmax = m-i+1;
205 540 : for (j=1;j<=jmax;j++) {
206 480 : ix = ix + 1;
207 480 : if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
208 480 : if (j!=jmax) {
209 420 : pd = cst*(PetscReal)(i+j-1);
210 : /* north */
211 420 : if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
212 364 : else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
213 : /* east */
214 420 : if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
215 364 : else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
216 : }
217 : /* south */
218 480 : pu = 0.5 - cst*(PetscReal)(i+j-3);
219 480 : if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
220 : /* west */
221 480 : if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
222 : }
223 : }
224 4 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
225 4 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
226 4 : PetscFunctionReturn(PETSC_SUCCESS);
227 : }
228 :
229 : /*
230 : ComputeResidualNorm - Computes the norm of the residual vector
231 : associated with an eigenpair.
232 :
233 : Input Parameters:
234 : trans - whether A' must be used instead of A
235 : kr,ki - eigenvalue
236 : xr,xi - eigenvector
237 : u - work vector
238 : */
239 34 : PetscErrorCode ComputeResidualNorm(Mat A,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec u,PetscReal *norm)
240 : {
241 : #if !defined(PETSC_USE_COMPLEX)
242 : PetscReal ni,nr;
243 : #endif
244 34 : PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultTranspose: MatMult;
245 :
246 34 : PetscFunctionBegin;
247 : #if !defined(PETSC_USE_COMPLEX)
248 : if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
249 : #endif
250 34 : PetscCall((*matmult)(A,xr,u));
251 34 : if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) PetscCall(VecAXPY(u,-kr,xr));
252 34 : PetscCall(VecNorm(u,NORM_2,norm));
253 : #if !defined(PETSC_USE_COMPLEX)
254 : } else {
255 : PetscCall((*matmult)(A,xr,u));
256 : if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
257 : PetscCall(VecAXPY(u,-kr,xr));
258 : PetscCall(VecAXPY(u,ki,xi));
259 : }
260 : PetscCall(VecNorm(u,NORM_2,&nr));
261 : PetscCall((*matmult)(A,xi,u));
262 : if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
263 : PetscCall(VecAXPY(u,-kr,xi));
264 : PetscCall(VecAXPY(u,-ki,xr));
265 : }
266 : PetscCall(VecNorm(u,NORM_2,&ni));
267 : *norm = SlepcAbsEigenvalue(nr,ni);
268 : }
269 : #endif
270 34 : PetscFunctionReturn(PETSC_SUCCESS);
271 : }
272 :
273 : /*TEST
274 :
275 : testset:
276 : args: -st_type sinvert -eps_target 1.1 -eps_nev 4
277 : filter: grep -v method | sed -e "s/[+-]0\.0*i//g" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
278 : requires: !single
279 : output_file: output/ex41_1.out
280 : test:
281 : suffix: 1
282 : args: -eps_type {{power krylovschur}}
283 : test:
284 : suffix: 1_balance
285 : args: -eps_balance {{oneside twoside}} -eps_ncv 17 -eps_krylovschur_locking 0
286 : requires: !__float128
287 :
288 : TEST*/
|
