Actual source code: ex18.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7:
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Solves the same problem as in ex5, but with a user-defined sorting criterion."
23: "It is a standard nonsymmetric eigenproblem with real eigenvalues and the rightmost eigenvalue is known to be 1.\n"
24: "This example illustrates how the user can set a custom spectrum selection.\n\n"
25: "The command line options are:\n"
26: " -m <m>, where <m> = number of grid subdivisions in each dimension.\n\n";
28: #include slepceps.h
30: /*
31: User-defined routines
32: */
34: PetscErrorCode MyEigenSort(EPS eps, PetscScalar ar, PetscScalar ai, PetscScalar br, PetscScalar bi, PetscInt *r, void *ctx);
36: PetscErrorCode MatMarkovModel( PetscInt m, Mat A );
40: int main( int argc, char **argv )
41: {
42: Vec v0; /* initial vector */
43: Mat A; /* operator matrix */
44: EPS eps; /* eigenproblem solver context */
45: const EPSType type;
46: PetscReal error, tol, re, im;
47: PetscScalar kr, ki, target=0.5;
48: PetscInt N, m=15, nev, maxit, i, its, nconv;
50:
51: SlepcInitialize(&argc,&argv,(char*)0,help);
53: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
54: N = m*(m+1)/2;
55: PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%d (m=%d)\n",N,m);
56: PetscOptionsGetScalar(PETSC_NULL,"-target",&target,PETSC_NULL);
57: PetscPrintf(PETSC_COMM_WORLD,"Searching closest eigenvalues to the right of %g.\n\n",target);
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Compute the operator matrix that defines the eigensystem, Ax=kx
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: MatCreate(PETSC_COMM_WORLD,&A);
64: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
65: MatSetFromOptions(A);
66: MatMarkovModel( m, A );
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Create the eigensolver and set various options
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: /*
73: Create eigensolver context
74: */
75: EPSCreate(PETSC_COMM_WORLD,&eps);
77: /*
78: Set operators. In this case, it is a standard eigenvalue problem
79: */
80: EPSSetOperators(eps,A,PETSC_NULL);
81: EPSSetProblemType(eps,EPS_NHEP);
83: /*
84: Set the custom comparing routine in order to obtain the eigenvalues
85: closest to the target on the right only
86: */
87: EPSSetEigenvalueComparison(eps,MyEigenSort,&target);
89: /*
90: Set solver parameters at runtime
91: */
92: EPSSetFromOptions(eps);
94: /*
95: Set the initial vector. This is optional, if not done the initial
96: vector is set to random values
97: */
98: MatGetVecs(A,&v0,PETSC_NULL);
99: VecSet(v0,1.0);
100: EPSSetInitialSpace(eps,1,&v0);
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Solve the eigensystem
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: EPSSolve(eps);
107: EPSGetIterationNumber(eps, &its);
108: PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);
110: /*
111: Optional: Get some information from the solver and display it
112: */
113: EPSGetType(eps,&type);
114: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
115: EPSGetDimensions(eps,&nev,PETSC_NULL,PETSC_NULL);
116: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev);
117: EPSGetTolerances(eps,&tol,&maxit);
118: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Display solution and clean up
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: /*
125: Get number of converged approximate eigenpairs
126: */
127: EPSGetConverged(eps,&nconv);
128: PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %d\n\n",nconv);
130: if (nconv>0) {
131: /*
132: Display eigenvalues and relative errors
133: */
134: PetscPrintf(PETSC_COMM_WORLD,
135: " k ||Ax-kx||/||kx||\n"
136: " ----------------- ------------------\n" );
138: for( i=0; i<nconv; i++ ) {
139: /*
140: Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
141: ki (imaginary part)
142: */
143: EPSGetEigenpair(eps,i,&kr,&ki,PETSC_NULL,PETSC_NULL);
144: /*
145: Compute the relative error associated to each eigenpair
146: */
147: EPSComputeRelativeError(eps,i,&error);
149: #ifdef PETSC_USE_COMPLEX
150: re = PetscRealPart(kr);
151: im = PetscImaginaryPart(kr);
152: #else
153: re = kr;
154: im = ki;
155: #endif
156: if (im!=0.0) {
157: PetscPrintf(PETSC_COMM_WORLD," %9f%+9f j %12g\n",re,im,error);
158: } else {
159: PetscPrintf(PETSC_COMM_WORLD," %12f %12g\n",re,error);
160: }
161: }
162: PetscPrintf(PETSC_COMM_WORLD,"\n" );
163: }
164:
165: /*
166: Free work space
167: */
168: EPSDestroy(eps);
169: MatDestroy(A);
170: VecDestroy(v0);
171: SlepcFinalize();
172: return 0;
173: }
177: /*
178: Matrix generator for a Markov model of a random walk on a triangular grid.
180: This subroutine generates a test matrix that models a random walk on a
181: triangular grid. This test example was used by G. W. Stewart ["{SRRIT} - a
182: FORTRAN subroutine to calculate the dominant invariant subspaces of a real
183: matrix", Tech. report. TR-514, University of Maryland (1978).] and in a few
184: papers on eigenvalue problems by Y. Saad [see e.g. LAA, vol. 34, pp. 269-295
185: (1980) ]. These matrices provide reasonably easy test problems for eigenvalue
186: algorithms. The transpose of the matrix is stochastic and so it is known
187: that one is an exact eigenvalue. One seeks the eigenvector of the transpose
188: associated with the eigenvalue unity. The problem is to calculate the steady
189: state probability distribution of the system, which is the eigevector
190: associated with the eigenvalue one and scaled in such a way that the sum all
191: the components is equal to one.
193: Note: the code will actually compute the transpose of the stochastic matrix
194: that contains the transition probabilities.
195: */
196: PetscErrorCode MatMarkovModel( PetscInt m, Mat A )
197: {
198: const PetscReal cst = 0.5/(PetscReal)(m-1);
199: PetscReal pd, pu;
200: PetscErrorCode ierr;
201: PetscInt Istart, Iend, i, j, jmax, ix=0;
204: MatGetOwnershipRange(A,&Istart,&Iend);
205: for( i=1; i<=m; i++ ) {
206: jmax = m-i+1;
207: for( j=1; j<=jmax; j++ ) {
208: ix = ix + 1;
209: if( ix-1<Istart || ix>Iend ) continue; /* compute only owned rows */
210: if( j!=jmax ) {
211: pd = cst*(PetscReal)(i+j-1);
212: /* north */
213: if( i==1 ) {
214: MatSetValue( A, ix-1, ix, 2*pd, INSERT_VALUES );
215: }
216: else {
217: MatSetValue( A, ix-1, ix, pd, INSERT_VALUES );
218: }
219: /* east */
220: if( j==1 ) {
221: MatSetValue( A, ix-1, ix+jmax-1, 2*pd, INSERT_VALUES );
222: }
223: else {
224: MatSetValue( A, ix-1, ix+jmax-1, pd, INSERT_VALUES );
225: }
226: }
227: /* south */
228: pu = 0.5 - cst*(PetscReal)(i+j-3);
229: if( j>1 ) {
230: MatSetValue( A, ix-1, ix-2, pu, INSERT_VALUES );
231: }
232: /* west */
233: if( i>1 ) {
234: MatSetValue( A, ix-1, ix-jmax-2, pu, INSERT_VALUES );
235: }
236: }
237: }
238: MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY );
239: MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY );
240: return(0);
241: }
245: /*
246: Function for user-defined eigenvalue ordering criterion.
248: Given two eigenvalues ar+i*ai and br+i*bi, the subroutine must choose
249: one of them as the preferred one according to the criterion.
250: In this example, the preferred value is the one closest to the target,
251: but on the right side.
252: */
253: PetscErrorCode MyEigenSort(EPS eps, PetscScalar ar, PetscScalar ai, PetscScalar br, PetscScalar bi, PetscInt *r, void *ctx)
254: {
255: PetscScalar target = *(PetscScalar*)ctx;
256: PetscReal da,db;
257: PetscTruth aisright, bisright;
261: if (PetscRealPart(target) < PetscRealPart(ar)) aisright = PETSC_TRUE;
262: else aisright = PETSC_FALSE;
263: if (PetscRealPart(target) < PetscRealPart(br)) bisright = PETSC_TRUE;
264: else bisright = PETSC_FALSE;
265: if (aisright == bisright) {
266: /* both are on the same side of the target */
267: da = SlepcAbsEigenvalue(ar-target,ai);
268: db = SlepcAbsEigenvalue(br-target,bi);
269: if (da < db) *r = -1;
270: else if (da > db) *r = 1;
271: else *r = 0;
272: } else if (aisright && !bisright)
273: *r = -1; /* 'a' is on the right */
274: else
275: *r = 1; /* 'b' is on the right */
277: return(0);
278: }