Actual source code: test12.c

slepc-3.18.2 2023-01-26
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*);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

154: /*TEST

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

162: TEST*/