Actual source code: ex27.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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[] = "Simple nonlinear eigenproblem using the NLEIGS solver.\n\n"
 12:   "The command line options are:\n"
 13:   "  -n <n>, where <n> = matrix dimension.\n"
 14:   "  -split <0/1>, to select the split form in the problem definition (enabled by default)\n";

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

 23: #include <slepcnep.h>

 25: /*
 26:    User-defined routines
 27: */
 28: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
 29: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);
 30: PetscErrorCode ComputeSingularities(NEP,PetscInt*,PetscScalar*,void*);

 32: int main(int argc,char **argv)
 33: {
 34:   NEP            nep;             /* nonlinear eigensolver context */
 35:   Mat            F,J,A[2];
 36:   NEPType        type;
 37:   PetscInt       n=100,nev,Istart,Iend,i;
 39:   PetscBool      terse,split=PETSC_TRUE;
 40:   RG             rg;
 41:   FN             f[2];
 42:   PetscScalar    coeffs;

 44:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 45:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 46:   PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL);
 47:   PetscPrintf(PETSC_COMM_WORLD,"\nSquare root eigenproblem, n=%D%s\n\n",n,split?" (in split form)":"");

 49:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 50:      Create nonlinear eigensolver context
 51:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 53:   NEPCreate(PETSC_COMM_WORLD,&nep);

 55:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 56:      Select the NLEIGS solver and set required options for it
 57:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 59:   NEPSetType(nep,NEPNLEIGS);
 60:   NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,NULL);
 61:   NEPGetRG(nep,&rg);
 62:   RGSetType(rg,RGINTERVAL);
 63: #if defined(PETSC_USE_COMPLEX)
 64:   RGIntervalSetEndpoints(rg,0.01,16.0,-0.001,0.001);
 65: #else
 66:   RGIntervalSetEndpoints(rg,0.01,16.0,0,0);
 67: #endif
 68:   NEPSetTarget(nep,1.1);

 70:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 71:      Define the nonlinear problem
 72:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 74:   if (split) {
 75:     /*
 76:        Create matrices for the split form
 77:     */
 78:     MatCreate(PETSC_COMM_WORLD,&A[0]);
 79:     MatSetSizes(A[0],PETSC_DECIDE,PETSC_DECIDE,n,n);
 80:     MatSetFromOptions(A[0]);
 81:     MatSetUp(A[0]);
 82:     MatGetOwnershipRange(A[0],&Istart,&Iend);
 83:     for (i=Istart;i<Iend;i++) {
 84:       if (i>0) { MatSetValue(A[0],i,i-1,1.0,INSERT_VALUES); }
 85:       if (i<n-1) { MatSetValue(A[0],i,i+1,1.0,INSERT_VALUES); }
 86:       MatSetValue(A[0],i,i,-2.0,INSERT_VALUES);
 87:     }
 88:     MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY);
 89:     MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY);

 91:     MatCreate(PETSC_COMM_WORLD,&A[1]);
 92:     MatSetSizes(A[1],PETSC_DECIDE,PETSC_DECIDE,n,n);
 93:     MatSetFromOptions(A[1]);
 94:     MatSetUp(A[1]);
 95:     MatAssemblyBegin(A[1],MAT_FINAL_ASSEMBLY);
 96:     MatAssemblyEnd(A[1],MAT_FINAL_ASSEMBLY);
 97:     MatShift(A[1],1.0);

 99:     /*
100:        Define funcions for the split form
101:      */
102:     FNCreate(PETSC_COMM_WORLD,&f[0]);
103:     FNSetType(f[0],FNRATIONAL);
104:     coeffs = 1.0;
105:     FNRationalSetNumerator(f[0],1,&coeffs);
106:     FNCreate(PETSC_COMM_WORLD,&f[1]);
107:     FNSetType(f[1],FNSQRT);
108:     NEPSetSplitOperator(nep,2,A,f,SUBSET_NONZERO_PATTERN);

110:   } else {
111:     /*
112:        Callback form: create matrix and set Function evaluation routine
113:      */
114:     MatCreate(PETSC_COMM_WORLD,&F);
115:     MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n);
116:     MatSetFromOptions(F);
117:     MatSeqAIJSetPreallocation(F,3,NULL);
118:     MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
119:     MatSetUp(F);
120:     NEPSetFunction(nep,F,F,FormFunction,NULL);

122:     MatCreate(PETSC_COMM_WORLD,&J);
123:     MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
124:     MatSetFromOptions(J);
125:     MatSeqAIJSetPreallocation(J,1,NULL);
126:     MatMPIAIJSetPreallocation(J,1,NULL,1,NULL);
127:     MatSetUp(J);
128:     NEPSetJacobian(nep,J,FormJacobian,NULL);
129:   }

131:   /*
132:      Set solver parameters at runtime
133:   */
134:   NEPSetFromOptions(nep);

136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137:                       Solve the eigensystem
138:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139:   NEPSolve(nep);
140:   NEPGetType(nep,&type);
141:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
142:   NEPGetDimensions(nep,&nev,NULL,NULL);
143:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);

145:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146:                     Display solution and clean up
147:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

149:   /* show detailed info unless -terse option is given by user */
150:   PetscOptionsHasName(NULL,NULL,"-terse",&terse);
151:   if (terse) {
152:     NEPErrorView(nep,NEP_ERROR_BACKWARD,NULL);
153:   } else {
154:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
155:     NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
156:     NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
157:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
158:   }
159:   NEPDestroy(&nep);
160:   if (split) {
161:     MatDestroy(&A[0]);
162:     MatDestroy(&A[1]);
163:     FNDestroy(&f[0]);
164:     FNDestroy(&f[1]);
165:   } else {
166:     MatDestroy(&F);
167:     MatDestroy(&J);
168:   }
169:   SlepcFinalize();
170:   return ierr;
171: }

173: /* ------------------------------------------------------------------- */
174: /*
175:    FormFunction - Computes Function matrix  T(lambda)
176: */
177: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
178: {
180:   PetscInt       i,n,col[3],Istart,Iend;
181:   PetscBool      FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
182:   PetscScalar    value[3],t;

185:   /*
186:      Compute Function entries and insert into matrix
187:   */
188:   t = PetscSqrtScalar(lambda);
189:   MatGetSize(fun,&n,NULL);
190:   MatGetOwnershipRange(fun,&Istart,&Iend);
191:   if (Istart==0) FirstBlock=PETSC_TRUE;
192:   if (Iend==n) LastBlock=PETSC_TRUE;
193:   value[0]=1.0; value[1]=t-2.0; value[2]=1.0;
194:   for (i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++) {
195:     col[0]=i-1; col[1]=i; col[2]=i+1;
196:     MatSetValues(fun,1,&i,3,col,value,INSERT_VALUES);
197:   }
198:   if (LastBlock) {
199:     i=n-1; col[0]=n-2; col[1]=n-1;
200:     MatSetValues(fun,1,&i,2,col,value,INSERT_VALUES);
201:   }
202:   if (FirstBlock) {
203:     i=0; col[0]=0; col[1]=1; value[0]=t-2.0; value[1]=1.0;
204:     MatSetValues(fun,1,&i,2,col,value,INSERT_VALUES);
205:   }

207:   /*
208:      Assemble matrix
209:   */
210:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
211:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
212:   if (fun != B) {
213:     MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
214:     MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
215:   }
216:   return(0);
217: }

219: /* ------------------------------------------------------------------- */
220: /*
221:    FormJacobian - Computes Jacobian matrix  T'(lambda)
222: */
223: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
224: {
226:   Vec            d;

229:   MatCreateVecs(jac,&d,NULL);
230:   VecSet(d,0.5/PetscSqrtScalar(lambda));
231:   MatDiagonalSet(jac,d,INSERT_VALUES);
232:   VecDestroy(&d);
233:   return(0);
234: }

236: /* ------------------------------------------------------------------- */
237: /*
238:    ComputeSingularities - Computes maxnp points (at most) in the complex plane where
239:    the function T(.) is not analytic.

241:    In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
242: */
243: PetscErrorCode ComputeSingularities(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *pt)
244: {
245:   PetscReal h;
246:   PetscInt  i;

249:   h = 11.0/(*maxnp-1);
250:   xi[0] = -1e-5; xi[*maxnp-1] = -1e+6;
251:   for (i=1;i<*maxnp-1;i++) xi[i] = -PetscPowReal(10,-5+h*i);
252:   return(0);
253: }

255: /*TEST

257:    test:
258:       suffix: 1
259:       args: -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse

261:    test:
262:       suffix: 2
263:       args: -split 0 -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse

265:    test:
266:       suffix: 3
267:       args: -nep_nev 3 -nep_tol 1e-8 -nep_nleigs_rk_shifts 1.06,1.1,1.12,1.15 -nep_conv_norm -terse
268:       requires: !single
269:       output_file: output/ex27_1.out

271:    test:
272:       suffix: 4
273:       args: -split 0 -nep_nev 3 -nep_nleigs_rk_shifts 1.06,1.1,1.12,1.15 -nep_nleigs_interpolation_degree 90 -terse
274:       output_file: output/ex27_2.out

276:    test:
277:       suffix: 5
278:       args: -nep_nev 3 -mat_type aijcusparse -terse
279:       requires: cuda
280:       output_file: output/ex27_1.out

282:    test:
283:       suffix: 6
284:       args: -split 0 -nep_nev 3 -mat_type aijcusparse -terse
285:       requires: cuda
286:       output_file: output/ex27_2.out

288:    test:
289:       suffix: 7
290:       args: -split 0 -nep_type ciss -rg_type ellipse -rg_ellipse_center 8 -rg_ellipse_radius .7 -rg_ellipse_vscale 0.1 -terse
291:       requires: complex

293:    test:
294:       suffix: 8
295:       args: -nep_type ciss -rg_type ellipse -rg_ellipse_center 8 -rg_ellipse_radius .7 -rg_ellipse_vscale 0.1 -terse
296:       requires: complex
297:       filter: sed -e "s/ (in split form)//"
298:       output_file: output/ex27_7.out

300:    test:
301:       suffix: 9
302:       args: -nep_nev 4 -n 20 -terse
303:       requires: !single

305: TEST*/