Actual source code: test12.c

slepc-3.17.2 2022-08-09
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[] = "Test some NLEIGS interface functions.\n\n"
 12:   "Based on ex27.c. The command line options are:\n"
 13:   "  -n <n>, where <n> = matrix dimension.\n";

 15: /*
 16:    Solve T(lambda)x=0 using NLEIGS solver
 17:       with T(lambda) = -D+sqrt(lambda)*I
 18:       where D is the Laplacian operator in 1 dimension
 19:       and with the interpolation interval [.01,16]
 20: */

 22: #include <slepcnep.h>

 24: /*
 25:    User-defined routines
 26: */
 27: PetscErrorCode ComputeSingularities(NEP,PetscInt*,PetscScalar*,void*);

 29: int main(int argc,char **argv)
 30: {
 31:   NEP            nep;             /* nonlinear eigensolver context */
 32:   Mat            A[2];
 33:   PetscInt       n=100,Istart,Iend,i,ns,nsin;
 34:   PetscBool      terse,fb;
 35:   RG             rg;
 36:   FN             f[2];
 37:   PetscScalar    coeffs,shifts[]={1.06,1.1,1.12,1.15},*rkshifts,val;
 38:   PetscErrorCode (*fsing)(NEP,PetscInt*,PetscScalar*,void*);

 40:   SlepcInitialize(&argc,&argv,(char*)0,help);
 41:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 42:   PetscPrintf(PETSC_COMM_WORLD,"\nSquare root eigenproblem, n=%" PetscInt_FMT "\n\n",n);

 44:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 45:      Create nonlinear eigensolver and set some options
 46:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 48:   NEPCreate(PETSC_COMM_WORLD,&nep);
 49:   NEPSetType(nep,NEPNLEIGS);
 50:   NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,NULL);
 51:   NEPGetRG(nep,&rg);
 52:   RGSetType(rg,RGINTERVAL);
 53: #if defined(PETSC_USE_COMPLEX)
 54:   RGIntervalSetEndpoints(rg,0.01,16.0,-0.001,0.001);
 55: #else
 56:   RGIntervalSetEndpoints(rg,0.01,16.0,0,0);
 57: #endif
 58:   NEPSetTarget(nep,1.1);

 60:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61:      Define the nonlinear problem in split form
 62:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 64:   /* Create matrices */
 65:   MatCreate(PETSC_COMM_WORLD,&A[0]);
 66:   MatSetSizes(A[0],PETSC_DECIDE,PETSC_DECIDE,n,n);
 67:   MatSetFromOptions(A[0]);
 68:   MatSetUp(A[0]);
 69:   MatGetOwnershipRange(A[0],&Istart,&Iend);
 70:   for (i=Istart;i<Iend;i++) {
 71:     if (i>0) MatSetValue(A[0],i,i-1,1.0,INSERT_VALUES);
 72:     if (i<n-1) MatSetValue(A[0],i,i+1,1.0,INSERT_VALUES);
 73:     MatSetValue(A[0],i,i,-2.0,INSERT_VALUES);
 74:   }
 75:   MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY);
 76:   MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY);

 78:   MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,1.0,&A[1]);

 80:   /* Define functions */
 81:   FNCreate(PETSC_COMM_WORLD,&f[0]);
 82:   FNSetType(f[0],FNRATIONAL);
 83:   coeffs = 1.0;
 84:   FNRationalSetNumerator(f[0],1,&coeffs);
 85:   FNCreate(PETSC_COMM_WORLD,&f[1]);
 86:   FNSetType(f[1],FNSQRT);
 87:   NEPSetSplitOperator(nep,2,A,f,SUBSET_NONZERO_PATTERN);

 89:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 90:                         Set some options
 91:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 93:   NEPNLEIGSSetFullBasis(nep,PETSC_FALSE);
 94:   NEPNLEIGSSetRKShifts(nep,4,shifts);
 95:   NEPSetFromOptions(nep);

 97:   NEPNLEIGSGetFullBasis(nep,&fb);
 98:   PetscPrintf(PETSC_COMM_WORLD," Using full basis = %s\n",fb?"true":"false");
 99:   NEPNLEIGSGetRKShifts(nep,&ns,&rkshifts);
100:   if (ns) {
101:     PetscPrintf(PETSC_COMM_WORLD," Using %" PetscInt_FMT " RK shifts =",ns);
102:     for (i=0;i<ns;i++) PetscPrintf(PETSC_COMM_WORLD," %g",(double)PetscRealPart(rkshifts[i]));
103:     PetscPrintf(PETSC_COMM_WORLD,"\n");
104:     PetscFree(rkshifts);
105:   }
106:   NEPNLEIGSGetSingularitiesFunction(nep,&fsing,NULL);
107:   nsin = 1;
108:   (*fsing)(nep,&nsin,&val,NULL);
109:   PetscPrintf(PETSC_COMM_WORLD," First returned singularity = %g\n",(double)PetscRealPart(val));

111:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112:                       Solve the eigensystem
113:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114:   NEPSolve(nep);

116:   /* show detailed info unless -terse option is given by user */
117:   PetscOptionsHasName(NULL,NULL,"-terse",&terse);
118:   if (terse) NEPErrorView(nep,NEP_ERROR_BACKWARD,NULL);
119:   else {
120:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
121:     NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
122:     NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
123:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
124:   }
125:   NEPDestroy(&nep);
126:   MatDestroy(&A[0]);
127:   MatDestroy(&A[1]);
128:   FNDestroy(&f[0]);
129:   FNDestroy(&f[1]);
130:   SlepcFinalize();
131:   return 0;
132: }

134: /* ------------------------------------------------------------------- */
135: /*
136:    ComputeSingularities - Computes maxnp points (at most) in the complex plane where
137:    the function T(.) is not analytic.

139:    In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
140: */
141: PetscErrorCode ComputeSingularities(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *pt)
142: {
143:   PetscReal h;
144:   PetscInt  i;

147:   h = 11.0/(*maxnp-1);
148:   xi[0] = -1e-5; xi[*maxnp-1] = -1e+6;
149:   for (i=1;i<*maxnp-1;i++) xi[i] = -PetscPowReal(10,-5+h*i);
150:   PetscFunctionReturn(0);
151: }

153: /*TEST

155:    test:
156:       suffix: 1
157:       args: -nep_nev 3 -nep_nleigs_interpolation_degree 20 -terse -nep_view
158:       requires: double
159:       filter: grep -v tolerance | sed -e "s/[+-]0\.0*i//g"

161: TEST*/