Actual source code: test21.c

slepc-3.19.2 2023-09-05
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[] = "Illustrates region filtering. "
 12:   "Based on ex5.\n"
 13:   "The command line options are:\n"
 14:   "  -m <m>, where <m> = number of grid subdivisions in each dimension.\n\n";

 16: #include <slepceps.h>

 18: /*
 19:    User-defined routines
 20: */
 21: PetscErrorCode MatMarkovModel(PetscInt m,Mat A);

 23: int main(int argc,char **argv)
 24: {
 25:   Mat            A;
 26:   EPS            eps;
 27:   ST             st;
 28:   RG             rg;
 29:   PetscReal      radius,tol=PETSC_SMALL;
 30:   PetscScalar    target=0.5,kr,ki;
 31:   PetscInt       N,m=15,nev,i,nconv;
 32:   PetscBool      checkfile;
 33:   char           filename[PETSC_MAX_PATH_LEN];
 34:   PetscViewer    viewer;
 35:   char           str[50];

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

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

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

 56:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 57:                 Create a standalone spectral transformation
 58:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 60:   PetscCall(STCreate(PETSC_COMM_WORLD,&st));
 61:   PetscCall(STSetType(st,STSINVERT));
 62:   PetscCall(STSetShift(st,target));

 64:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 65:                     Create a region for filtering
 66:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 68:   PetscCall(RGCreate(PETSC_COMM_WORLD,&rg));
 69:   PetscCall(RGSetType(rg,RGELLIPSE));
 70:   radius = (1.0-PetscRealPart(target))/2.0;
 71:   PetscCall(RGEllipseSetParameters(rg,target+radius,radius,0.1));

 73:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 74:                 Create the eigensolver and set various options
 75:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 77:   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
 78:   PetscCall(EPSSetST(eps,st));
 79:   PetscCall(EPSSetRG(eps,rg));
 80:   PetscCall(EPSSetOperators(eps,A,NULL));
 81:   PetscCall(EPSSetProblemType(eps,EPS_NHEP));
 82:   PetscCall(EPSSetTolerances(eps,tol,PETSC_DEFAULT));
 83:   PetscCall(EPSSetTarget(eps,target));
 84:   PetscCall(EPSSetFromOptions(eps));

 86:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 87:                Solve the eigensystem and display solution
 88:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 90:   PetscCall(EPSSolve(eps));
 91:   PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
 92:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
 93:   PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));

 95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:                    Check file containing the eigenvalues
 97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 98:   PetscCall(PetscOptionsGetString(NULL,NULL,"-checkfile",filename,sizeof(filename),&checkfile));
 99:   if (checkfile) {
100: #if defined(PETSC_HAVE_COMPLEX)
101:     PetscComplex *eigs,eval;

103:     PetscCall(EPSGetConverged(eps,&nconv));
104:     PetscCall(PetscMalloc1(nconv,&eigs));
105:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
106:     PetscCall(PetscViewerBinaryRead(viewer,eigs,nconv,NULL,PETSC_COMPLEX));
107:     PetscCall(PetscViewerDestroy(&viewer));
108:     for (i=0;i<nconv;i++) {
109:       PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,NULL,NULL));
110: #if defined(PETSC_USE_COMPLEX)
111:       eval = kr;
112: #else
113:       eval = PetscCMPLX(kr,ki);
114: #endif
115:       PetscCheck(eval==eigs[i],PETSC_COMM_WORLD,PETSC_ERR_FILE_UNEXPECTED,"Eigenvalues in the file do not match");
116:     }
117:     PetscCall(PetscFree(eigs));
118: #else
119:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"The -checkfile option requires C99 complex numbers");
120: #endif
121:   }

123:   PetscCall(EPSDestroy(&eps));
124:   PetscCall(STDestroy(&st));
125:   PetscCall(RGDestroy(&rg));
126:   PetscCall(MatDestroy(&A));
127:   PetscCall(SlepcFinalize());
128:   return 0;
129: }

131: PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
132: {
133:   const PetscReal cst = 0.5/(PetscReal)(m-1);
134:   PetscReal       pd,pu;
135:   PetscInt        Istart,Iend,i,j,jmax,ix=0;

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

165: /*TEST

167:    test:
168:       suffix: 1
169:       args: -eps_nev 4 -eps_ncv 20 -eps_view_values binary:myvalues.bin -checkfile myvalues.bin
170:       output_file: output/test11_1.out
171:       requires: !single c99_complex

173: TEST*/