Actual source code: test7.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[] = "Test the NLEIGS solver with shell matrices.\n\n"
 12:   "This is based on ex27.\n"
 13:   "The command line options are:\n"
 14:   "  -n <n>, where <n> = matrix dimension.\n"
 15:   "  -split <0/1>, to select the split form in the problem definition (enabled by default)\n";

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

 24: #include <slepcnep.h>

 26: /* User-defined routines */
 27: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
 28: PetscErrorCode ComputeSingularities(NEP,PetscInt*,PetscScalar*,void*);
 29: PetscErrorCode MatMult_A0(Mat,Vec,Vec);
 30: PetscErrorCode MatGetDiagonal_A0(Mat,Vec);
 31: PetscErrorCode MatDuplicate_A0(Mat,MatDuplicateOption,Mat*);
 32: PetscErrorCode MatMult_A1(Mat,Vec,Vec);
 33: PetscErrorCode MatGetDiagonal_A1(Mat,Vec);
 34: PetscErrorCode MatDuplicate_A1(Mat,MatDuplicateOption,Mat*);
 35: PetscErrorCode MatMult_F(Mat,Vec,Vec);
 36: PetscErrorCode MatGetDiagonal_F(Mat,Vec);
 37: PetscErrorCode MatDuplicate_F(Mat,MatDuplicateOption,Mat*);
 38: PetscErrorCode MatDestroy_F(Mat);

 40: typedef struct {
 41:   PetscScalar t;  /* square root of lambda */
 42: } MatCtx;

 44: int main(int argc,char **argv)
 45: {
 46:   NEP            nep;
 47:   KSP            *ksp;
 48:   PC             pc;
 49:   Mat            F,A[2];
 50:   NEPType        type;
 51:   PetscInt       i,n=100,nev,its,nsolve;
 52:   PetscReal      keep,tol=PETSC_SQRT_MACHINE_EPSILON/10;
 54:   RG             rg;
 55:   FN             f[2];
 56:   PetscMPIInt    size;
 57:   PetscBool      terse,flg,lock,split=PETSC_TRUE;
 58:   PetscScalar    coeffs;
 59:   MatCtx         *ctx;

 61:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 62:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 63:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
 64:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 65:   PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL);
 66:   PetscPrintf(PETSC_COMM_WORLD,"\nSquare root eigenproblem, n=%D%s\n\n",n,split?" (in split form)":"");

 68:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 69:      Create NEP context, configure NLEIGS with appropriate options
 70:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 72:   NEPCreate(PETSC_COMM_WORLD,&nep);
 73:   NEPSetType(nep,NEPNLEIGS);
 74:   NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,NULL);
 75:   NEPGetRG(nep,&rg);
 76:   RGSetType(rg,RGINTERVAL);
 77: #if defined(PETSC_USE_COMPLEX)
 78:   RGIntervalSetEndpoints(rg,0.01,16.0,-0.001,0.001);
 79: #else
 80:   RGIntervalSetEndpoints(rg,0.01,16.0,0,0);
 81: #endif
 82:   NEPSetTarget(nep,1.1);
 83:   NEPNLEIGSGetKSPs(nep,&nsolve,&ksp);
 84:   for (i=0;i<nsolve;i++) {
 85:    KSPSetType(ksp[i],KSPBICG);
 86:    KSPGetPC(ksp[i],&pc);
 87:    PCSetType(pc,PCJACOBI);
 88:    KSPSetTolerances(ksp[i],tol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
 89:   }

 91:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92:      Define the nonlinear problem
 93:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 95:   if (split) {
 96:     /* Create matrix A0 (tridiagonal) */
 97:     MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,NULL,&A[0]);
 98:     MatShellSetOperation(A[0],MATOP_MULT,(void(*)(void))MatMult_A0);
 99:     MatShellSetOperation(A[0],MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A0);
100:     MatShellSetOperation(A[0],MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A0);
101:     MatShellSetOperation(A[0],MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A0);

103:     /* Create matrix A0 (identity) */
104:     MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,NULL,&A[1]);
105:     MatShellSetOperation(A[1],MATOP_MULT,(void(*)(void))MatMult_A1);
106:     MatShellSetOperation(A[1],MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A1);
107:     MatShellSetOperation(A[1],MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A1);
108:     MatShellSetOperation(A[1],MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A1);

110:     /* Define funcions for the split form */
111:     FNCreate(PETSC_COMM_WORLD,&f[0]);
112:     FNSetType(f[0],FNRATIONAL);
113:     coeffs = 1.0;
114:     FNRationalSetNumerator(f[0],1,&coeffs);
115:     FNCreate(PETSC_COMM_WORLD,&f[1]);
116:     FNSetType(f[1],FNSQRT);
117:     NEPSetSplitOperator(nep,2,A,f,SUBSET_NONZERO_PATTERN);
118:   } else {
119:     /* Callback form: create shell matrix for F=A0+sqrt(lambda)*A1  */
120:     PetscNew(&ctx);
121:     MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,(void*)ctx,&F);
122:     MatShellSetOperation(F,MATOP_MULT,(void(*)(void))MatMult_F);
123:     MatShellSetOperation(F,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_F);
124:     MatShellSetOperation(F,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_F);
125:     MatShellSetOperation(F,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_F);
126:     MatShellSetOperation(F,MATOP_DESTROY,(void(*)(void))MatDestroy_F);
127:     /* Set Function evaluation routine */
128:     NEPSetFunction(nep,F,F,FormFunction,NULL);
129:   }

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

134:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135:                       Solve the eigensystem
136:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137:   NEPSolve(nep);
138:   NEPGetType(nep,&type);
139:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
140:   NEPGetDimensions(nep,&nev,NULL,NULL);
141:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
142:   PetscObjectTypeCompare((PetscObject)nep,NEPNLEIGS,&flg);
143:   if (flg) {
144:     NEPNLEIGSGetRestart(nep,&keep);
145:     NEPNLEIGSGetLocking(nep,&lock);
146:     NEPNLEIGSGetInterpolation(nep,&tol,&its);
147:     PetscPrintf(PETSC_COMM_WORLD," Restart factor is %3.2f",(double)keep);
148:     if (lock) { PetscPrintf(PETSC_COMM_WORLD," (locking activated)"); }
149:     PetscPrintf(PETSC_COMM_WORLD,"\n Divided diferences with tol=%6.2g maxit=%D\n",(double)tol,its);
150:     PetscPrintf(PETSC_COMM_WORLD,"\n");
151:   }

153:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154:                     Display solution and clean up
155:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

157:   /* show detailed info unless -terse option is given by user */
158:   PetscOptionsHasName(NULL,NULL,"-terse",&terse);
159:   if (terse) {
160:     NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
161:   } else {
162:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
163:     NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
164:     NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
165:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
166:   }
167:   NEPDestroy(&nep);
168:   if (split) {
169:     MatDestroy(&A[0]);
170:     MatDestroy(&A[1]);
171:     FNDestroy(&f[0]);
172:     FNDestroy(&f[1]);
173:   } else {
174:     MatDestroy(&F);
175:   }
176:   SlepcFinalize();
177:   return ierr;
178: }

180: /*
181:    FormFunction - Computes Function matrix  T(lambda)
182: */
183: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
184: {
186:   MatCtx         *ctxF;

189:   MatShellGetContext(fun,(void**)&ctxF);
190:   ctxF->t = PetscSqrtScalar(lambda);
191:   return(0);
192: }

194: /*
195:    ComputeSingularities - Computes maxnp points (at most) in the complex plane where
196:    the function T(.) is not analytic.

198:    In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
199: */
200: PetscErrorCode ComputeSingularities(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *pt)
201: {
202:   PetscReal h;
203:   PetscInt  i;

206:   h = 12.0/(*maxnp-1);
207:   xi[0] = -1e-6; xi[*maxnp-1] = -1e+6;
208:   for (i=1;i<*maxnp-1;i++) xi[i] = -PetscPowReal(10,-6+h*i);
209:   return(0);
210: }

212: /* -------------------------------- A0 ----------------------------------- */

214: PetscErrorCode MatMult_A0(Mat A,Vec x,Vec y)
215: {
216:   PetscErrorCode    ierr;
217:   PetscInt          i,n;
218:   const PetscScalar *px;
219:   PetscScalar       *py;

222:   VecGetArrayRead(x,&px);
223:   VecGetArray(y,&py);
224:   VecGetSize(x,&n);
225:   py[0] = -2.0*px[0]+px[1];
226:   for (i=1;i<n-1;i++) py[i] = px[i-1]-2.0*px[i]+px[i+1];
227:   py[n-1] = px[n-2]-2.0*px[n-1];
228:   VecRestoreArrayRead(x,&px);
229:   VecRestoreArray(y,&py);
230:   return(0);
231: }

233: PetscErrorCode MatGetDiagonal_A0(Mat A,Vec diag)
234: {

238:   VecSet(diag,-2.0);
239:   return(0);
240: }

242: PetscErrorCode MatDuplicate_A0(Mat A,MatDuplicateOption op,Mat *B)
243: {
244:   PetscInt       n;
245:   MPI_Comm       comm;

249:   MatGetSize(A,&n,NULL);
250:   PetscObjectGetComm((PetscObject)A,&comm);
251:   MatCreateShell(comm,n,n,n,n,NULL,B);
252:   MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_A0);
253:   MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A0);
254:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A0);
255:   MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A0);
256:   return(0);
257: }

259: /* -------------------------------- A1 ----------------------------------- */

261: PetscErrorCode MatMult_A1(Mat A,Vec x,Vec y)
262: {

266:   VecCopy(x,y);
267:   return(0);
268: }

270: PetscErrorCode MatGetDiagonal_A1(Mat A,Vec diag)
271: {

275:   VecSet(diag,1.0);
276:   return(0);
277: }

279: PetscErrorCode MatDuplicate_A1(Mat A,MatDuplicateOption op,Mat *B)
280: {
281:   PetscInt       n;
282:   MPI_Comm       comm;

286:   MatGetSize(A,&n,NULL);
287:   PetscObjectGetComm((PetscObject)A,&comm);
288:   MatCreateShell(comm,n,n,n,n,NULL,B);
289:   MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_A1);
290:   MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A1);
291:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A1);
292:   MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A1);
293:   return(0);
294: }

296: /* -------------------------------- F ----------------------------------- */

298: PetscErrorCode MatMult_F(Mat A,Vec x,Vec y)
299: {
300:   PetscErrorCode    ierr;
301:   PetscInt          i,n;
302:   const PetscScalar *px;
303:   PetscScalar       *py,d;
304:   MatCtx            *ctx;

307:   MatShellGetContext(A,(void**)&ctx);
308:   VecGetArrayRead(x,&px);
309:   VecGetArray(y,&py);
310:   VecGetSize(x,&n);
311:   d = -2.0+ctx->t;
312:   py[0] = d*px[0]+px[1];
313:   for (i=1;i<n-1;i++) py[i] = px[i-1]+d*px[i]+px[i+1];
314:   py[n-1] = px[n-2]+d*px[n-1];
315:   VecRestoreArrayRead(x,&px);
316:   VecRestoreArray(y,&py);
317:   return(0);
318: }

320: PetscErrorCode MatGetDiagonal_F(Mat A,Vec diag)
321: {
323:   MatCtx         *ctx;

326:   MatShellGetContext(A,(void**)&ctx);
327:   VecSet(diag,-2.0+ctx->t);
328:   return(0);
329: }

331: PetscErrorCode MatDuplicate_F(Mat A,MatDuplicateOption op,Mat *B)
332: {
333:   MatCtx         *actx,*bctx;
334:   PetscInt       n;
335:   MPI_Comm       comm;

339:   MatShellGetContext(A,(void**)&actx);
340:   MatGetSize(A,&n,NULL);
341:   PetscNew(&bctx);
342:   bctx->t = actx->t;
343:   PetscObjectGetComm((PetscObject)A,&comm);
344:   MatCreateShell(comm,n,n,n,n,(void*)bctx,B);
345:   MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_F);
346:   MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_F);
347:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_F);
348:   MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_F);
349:   MatShellSetOperation(*B,MATOP_DESTROY,(void(*)(void))MatDestroy_F);
350:   return(0);
351: }

353: PetscErrorCode MatDestroy_F(Mat A)
354: {
355:   MatCtx         *ctx;

359:   MatShellGetContext(A,(void**)&ctx);
360:   PetscFree(ctx);
361:   return(0);
362: }

364: /*TEST

366:    test:
367:       suffix: 1
368:       args: -nep_nev 3 -nep_tol 1e-8 -nep_nleigs_locking 0 -nep_nleigs_interpolation_degree 90 -nep_nleigs_interpolation_tol 1e-8 -nep_nleigs_restart 0.4 -terse
369:       requires: !single

371:    test:
372:       suffix: 2
373:       args: -split 0 -nep_nev 3 -nep_tol 1e-8 -terse
374:       requires: !single

376: TEST*/