Actual source code: test7.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 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;
 53:   RG             rg;
 54:   FN             f[2];
 55:   PetscBool      terse,flg,lock,split=PETSC_TRUE;
 56:   PetscScalar    coeffs;
 57:   MatCtx         *ctx;

 59:   SlepcInitialize(&argc,&argv,(char*)0,help);
 60:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 61:   PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL);
 62:   PetscPrintf(PETSC_COMM_WORLD,"\nSquare root eigenproblem, n=%" PetscInt_FMT "%s\n\n",n,split?" (in split form)":"");

 64:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 65:      Create NEP context, configure NLEIGS with appropriate options
 66:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

 87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88:      Define the nonlinear problem
 89:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 91:   if (split) {
 92:     /* Create matrix A0 (tridiagonal) */
 93:     MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,NULL,&A[0]);
 94:     MatShellSetOperation(A[0],MATOP_MULT,(void(*)(void))MatMult_A0);
 95:     MatShellSetOperation(A[0],MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A0);
 96:     MatShellSetOperation(A[0],MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A0);
 97:     MatShellSetOperation(A[0],MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A0);

 99:     /* Create matrix A0 (identity) */
100:     MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,NULL,&A[1]);
101:     MatShellSetOperation(A[1],MATOP_MULT,(void(*)(void))MatMult_A1);
102:     MatShellSetOperation(A[1],MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A1);
103:     MatShellSetOperation(A[1],MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A1);
104:     MatShellSetOperation(A[1],MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A1);

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

127:   /* Set solver parameters at runtime */
128:   NEPSetFromOptions(nep);

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

149:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150:                     Display solution and clean up
151:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

173: /*
174:    FormFunction - Computes Function matrix  T(lambda)
175: */
176: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
177: {
178:   MatCtx         *ctxF;

181:   MatShellGetContext(fun,&ctxF);
182:   ctxF->t = PetscSqrtScalar(lambda);
183:   PetscFunctionReturn(0);
184: }

186: /*
187:    ComputeSingularities - Computes maxnp points (at most) in the complex plane where
188:    the function T(.) is not analytic.

190:    In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
191: */
192: PetscErrorCode ComputeSingularities(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *pt)
193: {
194:   PetscReal h;
195:   PetscInt  i;

198:   h = 12.0/(*maxnp-1);
199:   xi[0] = -1e-6; xi[*maxnp-1] = -1e+6;
200:   for (i=1;i<*maxnp-1;i++) xi[i] = -PetscPowReal(10,-6+h*i);
201:   PetscFunctionReturn(0);
202: }

204: /* -------------------------------- A0 ----------------------------------- */

206: PetscErrorCode MatMult_A0(Mat A,Vec x,Vec y)
207: {
208:   PetscInt          i,n;
209:   PetscMPIInt       rank,size,next,prev;
210:   const PetscScalar *px;
211:   PetscScalar       *py,upper=0.0,lower=0.0;
212:   MPI_Comm          comm;

215:   PetscObjectGetComm((PetscObject)A,&comm);
216:   MPI_Comm_size(comm,&size);
217:   MPI_Comm_rank(comm,&rank);
218:   next = rank==size-1? MPI_PROC_NULL: rank+1;
219:   prev = rank==0? MPI_PROC_NULL: rank-1;

221:   VecGetArrayRead(x,&px);
222:   VecGetArray(y,&py);
223:   VecGetLocalSize(x,&n);

225:   MPI_Sendrecv(px,1,MPIU_SCALAR,prev,0,&lower,1,MPIU_SCALAR,next,0,comm,MPI_STATUS_IGNORE);
226:   MPI_Sendrecv(px+n-1,1,MPIU_SCALAR,next,0,&upper,1,MPIU_SCALAR,prev,0,comm,MPI_STATUS_IGNORE);

228:   py[0] = upper-2.0*px[0]+px[1];
229:   for (i=1;i<n-1;i++) py[i] = px[i-1]-2.0*px[i]+px[i+1];
230:   py[n-1] = px[n-2]-2.0*px[n-1]+lower;
231:   VecRestoreArrayRead(x,&px);
232:   VecRestoreArray(y,&py);
233:   PetscFunctionReturn(0);
234: }

236: PetscErrorCode MatGetDiagonal_A0(Mat A,Vec diag)
237: {
239:   VecSet(diag,-2.0);
240:   PetscFunctionReturn(0);
241: }

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

248:   MatGetSize(A,&M,&N);
249:   MatGetLocalSize(A,&m,&n);
250:   PetscObjectGetComm((PetscObject)A,&comm);
251:   MatCreateShell(comm,m,n,M,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:   PetscFunctionReturn(0);
257: }

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

261: PetscErrorCode MatMult_A1(Mat A,Vec x,Vec y)
262: {
264:   VecCopy(x,y);
265:   PetscFunctionReturn(0);
266: }

268: PetscErrorCode MatGetDiagonal_A1(Mat A,Vec diag)
269: {
271:   VecSet(diag,1.0);
272:   PetscFunctionReturn(0);
273: }

275: PetscErrorCode MatDuplicate_A1(Mat A,MatDuplicateOption op,Mat *B)
276: {
277:   PetscInt       m,n,M,N;
278:   MPI_Comm       comm;

280:   MatGetSize(A,&M,&N);
281:   MatGetLocalSize(A,&m,&n);
282:   PetscObjectGetComm((PetscObject)A,&comm);
283:   MatCreateShell(comm,m,n,M,N,NULL,B);
284:   MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_A1);
285:   MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A1);
286:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A1);
287:   MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A1);
288:   PetscFunctionReturn(0);
289: }

291: /* -------------------------------- F ----------------------------------- */

293: PetscErrorCode MatMult_F(Mat A,Vec x,Vec y)
294: {
295:   PetscInt          i,n;
296:   PetscMPIInt       rank,size,next,prev;
297:   const PetscScalar *px;
298:   PetscScalar       *py,d,upper=0.0,lower=0.0;
299:   MatCtx            *ctx;
300:   MPI_Comm          comm;

303:   PetscObjectGetComm((PetscObject)A,&comm);
304:   MPI_Comm_size(comm,&size);
305:   MPI_Comm_rank(comm,&rank);
306:   next = rank==size-1? MPI_PROC_NULL: rank+1;
307:   prev = rank==0? MPI_PROC_NULL: rank-1;

309:   MatShellGetContext(A,&ctx);
310:   VecGetArrayRead(x,&px);
311:   VecGetArray(y,&py);
312:   VecGetLocalSize(x,&n);

314:   MPI_Sendrecv(px,1,MPIU_SCALAR,prev,0,&lower,1,MPIU_SCALAR,next,0,comm,MPI_STATUS_IGNORE);
315:   MPI_Sendrecv(px+n-1,1,MPIU_SCALAR,next,0,&upper,1,MPIU_SCALAR,prev,0,comm,MPI_STATUS_IGNORE);

317:   d = -2.0+ctx->t;
318:   py[0] = upper+d*px[0]+px[1];
319:   for (i=1;i<n-1;i++) py[i] = px[i-1]+d*px[i]+px[i+1];
320:   py[n-1] = px[n-2]+d*px[n-1]+lower;
321:   VecRestoreArrayRead(x,&px);
322:   VecRestoreArray(y,&py);
323:   PetscFunctionReturn(0);
324: }

326: PetscErrorCode MatGetDiagonal_F(Mat A,Vec diag)
327: {
328:   MatCtx         *ctx;

331:   MatShellGetContext(A,&ctx);
332:   VecSet(diag,-2.0+ctx->t);
333:   PetscFunctionReturn(0);
334: }

336: PetscErrorCode MatDuplicate_F(Mat A,MatDuplicateOption op,Mat *B)
337: {
338:   MatCtx         *actx,*bctx;
339:   PetscInt       m,n,M,N;
340:   MPI_Comm       comm;

342:   MatShellGetContext(A,&actx);
343:   MatGetSize(A,&M,&N);
344:   MatGetLocalSize(A,&m,&n);
345:   PetscNew(&bctx);
346:   bctx->t = actx->t;
347:   PetscObjectGetComm((PetscObject)A,&comm);
348:   MatCreateShell(comm,m,n,M,N,(void*)bctx,B);
349:   MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_F);
350:   MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_F);
351:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_F);
352:   MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_F);
353:   MatShellSetOperation(*B,MATOP_DESTROY,(void(*)(void))MatDestroy_F);
354:   PetscFunctionReturn(0);
355: }

357: PetscErrorCode MatDestroy_F(Mat A)
358: {
359:   MatCtx         *ctx;

361:   MatShellGetContext(A,&ctx);
362:   PetscFree(ctx);
363:   PetscFunctionReturn(0);
364: }

366: /*TEST

368:    testset:
369:       nsize: {{1 2}}
370:       args: -nep_nev 3 -nep_tol 1e-8 -terse
371:       filter: sed -e "s/[+-]0\.0*i//g"
372:       requires: !single
373:       test:
374:          suffix: 1
375:          args: -nep_nleigs_locking 0 -nep_nleigs_interpolation_degree 90 -nep_nleigs_interpolation_tol 1e-8 -nep_nleigs_restart 0.4
376:       test:
377:          suffix: 2
378:          args: -split 0

380: TEST*/