Actual source code: ex5.c

  2: static char help[] = "Eigenvalue problem associated with a Markov model of a random walk on a triangular grid. "
  3:   "It is a standard nonsymmetric eigenproblem with real eigenvalues and the rightmost eigenvalue is known to be 1.\n"
  4:   "This example illustrates how the user can set the initial vector.\n\n"
  5:   "The command line options are:\n"
  6:   "  -m <m>, where <m> = number of grid subdivisions in each dimension.\n\n";

 8:  #include slepceps.h

 10: /* 
 11:    User-defined routines
 12: */
 13: PetscErrorCode MatMarkovModel( PetscInt m, Mat A );

 17: int main( int argc, char **argv )
 18: {
 19:   Vec                  v0;                  /* initial vector */
 20:   Mat                  A;                  /* operator matrix */
 21:   EPS                  eps;                  /* eigenproblem solver context */
 22:   EPSType              type;
 23:   PetscReal            error, tol, re, im;
 24:   PetscScalar          kr, ki;
 25:   PetscInt             N, m=15;
 26:   int                  nev, maxit, i, its, nconv;
 28: 
 29:   SlepcInitialize(&argc,&argv,(char*)0,help);

 31:   PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
 32:   N = m*(m+1)/2;
 33:   PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%d (m=%d)\n\n",N,m);

 35:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 36:      Compute the operator matrix that defines the eigensystem, Ax=kx
 37:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 39:   MatCreate(PETSC_COMM_WORLD,&A);
 40:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
 41:   MatSetFromOptions(A);
 42:   MatMarkovModel( m, A );

 44:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 45:                 Create the eigensolver and set various options
 46:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 48:   /* 
 49:      Create eigensolver context
 50:   */
 51:   EPSCreate(PETSC_COMM_WORLD,&eps);

 53:   /* 
 54:      Set operators. In this case, it is a standard eigenvalue problem
 55:   */
 56:   EPSSetOperators(eps,A,PETSC_NULL);
 57:   EPSSetProblemType(eps,EPS_NHEP);

 59:   /*
 60:      Set solver parameters at runtime
 61:   */
 62:   EPSSetFromOptions(eps);

 64:   /*
 65:      Set the initial vector. This is optional, if not done the initial
 66:      vector is set to random values
 67:   */
 68:   MatGetVecs(A,&v0,PETSC_NULL);
 69:   VecSet(v0,1.0);
 70:   EPSSetInitialVector(eps,v0);

 72:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 73:                       Solve the eigensystem
 74:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 76:   EPSSolve(eps);
 77:   EPSGetIterationNumber(eps, &its);
 78:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);

 80:   /*
 81:      Optional: Get some information from the solver and display it
 82:   */
 83:   EPSGetType(eps,&type);
 84:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
 85:   EPSGetDimensions(eps,&nev,PETSC_NULL);
 86:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev);
 87:   EPSGetTolerances(eps,&tol,&maxit);
 88:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);

 90:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 91:                     Display solution and clean up
 92:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 94:   /* 
 95:      Get number of converged approximate eigenpairs
 96:   */
 97:   EPSGetConverged(eps,&nconv);
 98:   PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %d\n\n",nconv);

100:   if (nconv>0) {
101:     /*
102:        Display eigenvalues and relative errors
103:     */
104:     PetscPrintf(PETSC_COMM_WORLD,
105:          "           k          ||Ax-kx||/||kx||\n"
106:          "   ----------------- ------------------\n" );

108:     for( i=0; i<nconv; i++ ) {
109:       /* 
110:         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
111:         ki (imaginary part)
112:       */
113:       EPSGetEigenpair(eps,i,&kr,&ki,PETSC_NULL,PETSC_NULL);
114:       /*
115:          Compute the relative error associated to each eigenpair
116:       */
117:       EPSComputeRelativeError(eps,i,&error);

119: #ifdef PETSC_USE_COMPLEX
120:       re = PetscRealPart(kr);
121:       im = PetscImaginaryPart(kr);
122: #else
123:       re = kr;
124:       im = ki;
125: #endif 
126:       if (im!=0.0) {
127:         PetscPrintf(PETSC_COMM_WORLD," %9f%+9f j %12g\n",re,im,error);
128:       } else {
129:         PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12g\n",re,error);
130:       }
131:     }
132:     PetscPrintf(PETSC_COMM_WORLD,"\n" );
133:   }
134: 
135:   /* 
136:      Free work space
137:   */
138:   EPSDestroy(eps);
139:   MatDestroy(A);
140:   VecDestroy(v0);
141:   SlepcFinalize();
142:   return 0;
143: }

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

150:     This subroutine generates a test matrix that models a random walk on a 
151:     triangular grid. This test example was used by G. W. Stewart ["{SRRIT} - a 
152:     FORTRAN subroutine to calculate the dominant invariant subspaces of a real
153:     matrix", Tech. report. TR-514, University of Maryland (1978).] and in a few
154:     papers on eigenvalue problems by Y. Saad [see e.g. LAA, vol. 34, pp. 269-295
155:     (1980) ]. These matrices provide reasonably easy test problems for eigenvalue
156:     algorithms. The transpose of the matrix  is stochastic and so it is known 
157:     that one is an exact eigenvalue. One seeks the eigenvector of the transpose 
158:     associated with the eigenvalue unity. The problem is to calculate the steady
159:     state probability distribution of the system, which is the eigevector 
160:     associated with the eigenvalue one and scaled in such a way that the sum all
161:     the components is equal to one.

163:     Note: the code will actually compute the transpose of the stochastic matrix
164:     that contains the transition probabilities.
165: */
166: PetscErrorCode MatMarkovModel( PetscInt m, Mat A )
167: {
168:   const PetscReal cst = 0.5/(PetscReal)(m-1);
169:   PetscReal       pd, pu;
170:   PetscErrorCode  ierr;
171:   PetscInt        Istart, Iend, i, j, jmax, ix=0;

174:   MatGetOwnershipRange(A,&Istart,&Iend);
175:   for( i=1; i<=m; i++ ) {
176:     jmax = m-i+1;
177:     for( j=1; j<=jmax; j++ ) {
178:       ix = ix + 1;
179:       if( ix-1<Istart || ix>Iend ) continue;  /* compute only owned rows */
180:       if( j!=jmax ) {
181:         pd = cst*(PetscReal)(i+j-1);
182:         /* north */
183:         if( i==1 ) {
184:           MatSetValue( A, ix-1, ix, 2*pd, INSERT_VALUES );
185:         }
186:         else {
187:           MatSetValue( A, ix-1, ix, pd, INSERT_VALUES );
188:         }
189:         /* east */
190:         if( j==1 ) {
191:           MatSetValue( A, ix-1, ix+jmax-1, 2*pd, INSERT_VALUES );
192:         }
193:         else {
194:           MatSetValue( A, ix-1, ix+jmax-1, pd, INSERT_VALUES );
195:         }
196:       }
197:       /* south */
198:       pu = 0.5 - cst*(PetscReal)(i+j-3);
199:       if( j>1 ) {
200:         MatSetValue( A, ix-1, ix-2, pu, INSERT_VALUES );
201:       }
202:       /* west */
203:       if( i>1 ) {
204:         MatSetValue( A, ix-1, ix-jmax-2, pu, INSERT_VALUES );
205:       }
206:     }
207:   }
208:   MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY );
209:   MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY );
210:   return(0);
211: }