Actual source code: ex12.c

  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2013, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.

  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 eigenproblem as in example ex5, but computing also left eigenvectors. "
 23:   "It is a Markov model of a random walk on a triangular grid. "
 24:   "A standard nonsymmetric eigenproblem with real eigenvalues. The rightmost eigenvalue is known to be 1.\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: */
 33: PetscErrorCode MatMarkovModel(PetscInt m,Mat A);

 37: int main (int argc,char **argv)
 38: {
 39:   Vec            v0,w0;           /* initial vector */
 40:   Vec            *X,*Y;           /* right and left eigenvectors */
 41:   Mat            A;               /* operator matrix */
 42:   EPS            eps;             /* eigenproblem solver context */
 43:   EPSType        type;
 44:   PetscReal      error1,error2,tol,re,im;
 45:   PetscScalar    kr,ki;
 46:   PetscInt       nev,maxit,i,its,nconv,N,m=15;

 49:   SlepcInitialize(&argc,&argv,(char*)0,help);

 51:   PetscOptionsGetInt(NULL,"-m",&m,NULL);
 52:   N = m*(m+1)/2;
 53:   PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%D (m=%D)\n\n",N,m);

 55:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 56:      Compute the operator matrix that defines the eigensystem, Ax=kx
 57:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 59:   MatCreate(PETSC_COMM_WORLD,&A);
 60:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
 61:   MatSetFromOptions(A);
 62:   MatSetUp(A);
 63:   MatMarkovModel(m,A);

 65:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 66:                 Create the eigensolver and set various options
 67:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 69:   /*
 70:      Create eigensolver context
 71:   */
 72:   EPSCreate(PETSC_COMM_WORLD,&eps);

 74:   /*
 75:      Set operators. In this case, it is a standard eigenvalue problem.
 76:      Request also left eigenvectors
 77:   */
 78:   EPSSetOperators(eps,A,NULL);
 79:   EPSSetProblemType(eps,EPS_NHEP);
 80:   EPSSetLeftVectorsWanted(eps,PETSC_TRUE);

 82:   /*
 83:      Set solver parameters at runtime
 84:   */
 85:   EPSSetFromOptions(eps);

 87:   /*
 88:      Set the initial vector. This is optional, if not done the initial
 89:      vector is set to random values
 90:   */
 91:   MatGetVecs(A,&v0,&w0);
 92:   VecSet(v0,1.0);
 93:   MatMult(A,v0,w0);
 94:   EPSSetInitialSpace(eps,1,&v0);
 95:   EPSSetInitialSpaceLeft(eps,1,&w0);

 97:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 98:                       Solve the eigensystem
 99:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

101:   EPSSolve(eps);
102:   EPSGetIterationNumber(eps,&its);
103:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);

105:   /*
106:      Optional: Get some information from the solver and display it
107:   */
108:   EPSGetType(eps,&type);
109:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
110:   EPSGetDimensions(eps,&nev,NULL,NULL);
111:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
112:   EPSGetTolerances(eps,&tol,&maxit);
113:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4G, maxit=%D\n",tol,maxit);

115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116:                     Display solution and clean up
117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

119:   /*
120:      Get number of converged approximate eigenpairs
121:   */
122:   EPSGetConverged(eps,&nconv);
123:   PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %D\n\n",nconv);

125:   if (nconv>0) {
126:     /*
127:        Display eigenvalues and relative errors
128:     */
129:     PetscPrintf(PETSC_COMM_WORLD,
130:          "           k          ||Ax-kx||/||kx||   ||y'A-ky'||/||ky||\n"
131:          "   ----------------- ------------------ --------------------\n");

133:     for (i=0;i<nconv;i++) {
134:       /*
135:         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
136:         ki (imaginary part)
137:       */
138:       EPSGetEigenvalue(eps,i,&kr,&ki);
139:       /*
140:          Compute the relative errors associated to both right and left eigenvectors
141:       */
142:       EPSComputeRelativeError(eps,i,&error1);
143:       EPSComputeRelativeErrorLeft(eps,i,&error2);

145: #if defined(PETSC_USE_COMPLEX)
146:       re = PetscRealPart(kr);
147:       im = PetscImaginaryPart(kr);
148: #else
149:       re = kr;
150:       im = ki;
151: #endif
152:       if (im!=0.0) {
153:         PetscPrintf(PETSC_COMM_WORLD," %9F%+9F j %12G%12G\n",re,im,error1,error2);
154:       } else {
155:         PetscPrintf(PETSC_COMM_WORLD,"   %12F       %12G       %12G\n",re,error1,error2);
156:       }
157:     }
158:     PetscPrintf(PETSC_COMM_WORLD,"\n");

160:     VecDuplicateVecs(v0,nconv,&X);
161:     VecDuplicateVecs(w0,nconv,&Y);
162:     for (i=0;i<nconv;i++) {
163:       EPSGetEigenvector(eps,i,X[i],NULL);
164:       EPSGetEigenvectorLeft(eps,i,Y[i],NULL);
165:     }
166:     PetscPrintf(PETSC_COMM_WORLD,
167:          "                   Bi-orthogonality <x,y>                   \n"
168:          "   ---------------------------------------------------------\n");

170:     SlepcCheckOrthogonality(X,nconv,Y,nconv,NULL,NULL,NULL);
171:     PetscPrintf(PETSC_COMM_WORLD,"\n");
172:     VecDestroyVecs(nconv,&X);
173:     VecDestroyVecs(nconv,&Y);

175:   }

177:   /*
178:      Free work space
179:   */
180:   VecDestroy(&v0);
181:   VecDestroy(&w0);
182:   EPSDestroy(&eps);
183:   MatDestroy(&A);
184:   SlepcFinalize();
185:   return 0;
186: }

190: /*
191:     Matrix generator for a Markov model of a random walk on a triangular grid.

193:     This subroutine generates a test matrix that models a random walk on a
194:     triangular grid. This test example was used by G. W. Stewart ["{SRRIT} - a
195:     FORTRAN subroutine to calculate the dominant invariant subspaces of a real
196:     matrix", Tech. report. TR-514, University of Maryland (1978).] and in a few
197:     papers on eigenvalue problems by Y. Saad [see e.g. LAA, vol. 34, pp. 269-295
198:     (1980) ]. These matrices provide reasonably easy test problems for eigenvalue
199:     algorithms. The transpose of the matrix  is stochastic and so it is known
200:     that one is an exact eigenvalue. One seeks the eigenvector of the transpose
201:     associated with the eigenvalue unity. The problem is to calculate the steady
202:     state probability distribution of the system, which is the eigevector
203:     associated with the eigenvalue one and scaled in such a way that the sum all
204:     the components is equal to one.
205:     Note: the code will actually compute the transpose of the stochastic matrix
206:     that contains the transition probabilities.
207: */
208: PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
209: {
210:   const PetscReal cst = 0.5/(PetscReal)(m-1);
211:   PetscReal       pd,pu;
212:   PetscInt        i,j,jmax,ix=0,Istart,Iend;
213:   PetscErrorCode  ierr;

216:   MatGetOwnershipRange(A,&Istart,&Iend);
217:   for (i=1;i<=m;i++) {
218:     jmax = m-i+1;
219:     for (j=1;j<=jmax;j++) {
220:       ix = ix + 1;
221:       if (ix-1<Istart || ix>Iend) continue;  /* compute only owned rows */
222:       if (j!=jmax) {
223:         pd = cst*(PetscReal)(i+j-1);
224:         /* north */
225:         if (i==1) {
226:           MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES);
227:         } else {
228:           MatSetValue(A,ix-1,ix,pd,INSERT_VALUES);
229:         }
230:         /* east */
231:         if (j==1) {
232:           MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES);
233:         } else {
234:           MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES);
235:         }
236:       }
237:       /* south */
238:       pu = 0.5 - cst*(PetscReal)(i+j-3);
239:       if (j>1) {
240:         MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES);
241:       }
242:       /* west */
243:       if (i>1) {
244:         MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES);
245:       }
246:     }
247:   }
248:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
249:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
250:   return(0);
251: }