Actual source code: ex18.c

slepc-main 2024-11-15
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

 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";

 17: #include <slepceps.h>

 19: /*
 20:    User-defined routines
 21: */

 23: PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx);
 24: PetscErrorCode MatMarkovModel(PetscInt m,Mat A);

 26: int main(int argc,char **argv)
 27: {
 28:   Mat            A;               /* operator matrix */
 29:   EPS            eps;             /* eigenproblem solver context */
 30:   EPSType        type;
 31:   PetscScalar    target=0.5;
 32:   PetscInt       N,m=15,nev;
 33:   PetscBool      terse;
 34:   PetscViewer    viewer;
 35:   char           str[50];

 37:   PetscFunctionBeginUser;
 38:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));

 40:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
 41:   N = m*(m+1)/2;
 42:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n",N,m));
 43:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-target",&target,NULL));
 44:   PetscCall(SlepcSNPrintfScalar(str,sizeof(str),target,PETSC_FALSE));
 45:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Searching closest eigenvalues to the right of %s.\n\n",str));

 47:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 48:      Compute the operator matrix that defines the eigensystem, Ax=kx
 49:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 51:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
 52:   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
 53:   PetscCall(MatSetFromOptions(A));
 54:   PetscCall(MatMarkovModel(m,A));

 56:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 57:                 Create the eigensolver and set various options
 58:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 60:   /*
 61:      Create eigensolver context
 62:   */
 63:   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));

 65:   /*
 66:      Set operators. In this case, it is a standard eigenvalue problem
 67:   */
 68:   PetscCall(EPSSetOperators(eps,A,NULL));
 69:   PetscCall(EPSSetProblemType(eps,EPS_NHEP));

 71:   /*
 72:      Set the custom comparing routine in order to obtain the eigenvalues
 73:      closest to the target on the right only
 74:   */
 75:   PetscCall(EPSSetEigenvalueComparison(eps,MyEigenSort,&target));

 77:   /*
 78:      Set solver parameters at runtime
 79:   */
 80:   PetscCall(EPSSetFromOptions(eps));

 82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 83:                       Solve the eigensystem
 84:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 86:   PetscCall(EPSSolve(eps));

 88:   /*
 89:      Optional: Get some information from the solver and display it
 90:   */
 91:   PetscCall(EPSGetType(eps,&type));
 92:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
 93:   PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
 94:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));

 96:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 97:                     Display solution and clean up
 98:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

100:   /* show detailed info unless -terse option is given by user */
101:   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
102:   if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
103:   else {
104:     PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
105:     PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
106:     PetscCall(EPSConvergedReasonView(eps,viewer));
107:     PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,viewer));
108:     PetscCall(PetscViewerPopFormat(viewer));
109:   }
110:   PetscCall(EPSDestroy(&eps));
111:   PetscCall(MatDestroy(&A));
112:   PetscCall(SlepcFinalize());
113:   return 0;
114: }

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

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.

132:     Note: the code will actually compute the transpose of the stochastic matrix
133:     that contains the transition probabilities.
134: */
135: PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
136: {
137:   const PetscReal cst = 0.5/(PetscReal)(m-1);
138:   PetscReal       pd,pu;
139:   PetscInt        Istart,Iend,i,j,jmax,ix=0;

141:   PetscFunctionBeginUser;
142:   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
143:   for (i=1;i<=m;i++) {
144:     jmax = m-i+1;
145:     for (j=1;j<=jmax;j++) {
146:       ix = ix + 1;
147:       if (ix-1<Istart || ix>Iend) continue;  /* compute only owned rows */
148:       if (j!=jmax) {
149:         pd = cst*(PetscReal)(i+j-1);
150:         /* north */
151:         if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
152:         else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
153:         /* east */
154:         if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
155:         else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
156:       }
157:       /* south */
158:       pu = 0.5 - cst*(PetscReal)(i+j-3);
159:       if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
160:       /* west */
161:       if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
162:     }
163:   }
164:   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
165:   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
166:   PetscFunctionReturn(PETSC_SUCCESS);
167: }

169: /*
170:     Function for user-defined eigenvalue ordering criterion.

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: PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx)
178: {
179:   PetscScalar target = *(PetscScalar*)ctx;
180:   PetscReal   da,db;
181:   PetscBool   aisright,bisright;

183:   PetscFunctionBeginUser;
184:   if (PetscRealPart(target) < PetscRealPart(ar)) aisright = PETSC_TRUE;
185:   else aisright = PETSC_FALSE;
186:   if (PetscRealPart(target) < PetscRealPart(br)) bisright = PETSC_TRUE;
187:   else bisright = PETSC_FALSE;
188:   if (aisright == bisright) {
189:     /* both are on the same side of the target */
190:     da = SlepcAbsEigenvalue(ar-target,ai);
191:     db = SlepcAbsEigenvalue(br-target,bi);
192:     if (da < db) *r = -1;
193:     else if (da > db) *r = 1;
194:     else *r = 0;
195:   } else if (aisright && !bisright) *r = -1; /* 'a' is on the right */
196:   else *r = 1;  /* 'b' is on the right */
197:   PetscFunctionReturn(PETSC_SUCCESS);
198: }

200: /*TEST

202:    test:
203:       suffix: 1
204:       args: -eps_nev 4 -terse
205:       requires: !single
206:       filter: sed -e "s/[+-]0\.0*i//g"

208: TEST*/