Actual source code: test7.c

slepc-3.18.0 2022-10-01
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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

205: /* -------------------------------- A0 ----------------------------------- */

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

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

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

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

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

237: PetscErrorCode MatGetDiagonal_A0(Mat A,Vec diag)
238: {
240:   VecSet(diag,-2.0);
241:   return 0;
242: }

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

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

260: /* -------------------------------- A1 ----------------------------------- */

262: PetscErrorCode MatMult_A1(Mat A,Vec x,Vec y)
263: {
265:   VecCopy(x,y);
266:   return 0;
267: }

269: PetscErrorCode MatGetDiagonal_A1(Mat A,Vec diag)
270: {
272:   VecSet(diag,1.0);
273:   return 0;
274: }

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

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

292: /* -------------------------------- F ----------------------------------- */

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

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

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

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

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

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

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

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

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

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

362:   MatShellGetContext(A,&ctx);
363:   PetscFree(ctx);
364:   return 0;
365: }

367: /*TEST

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

381: TEST*/