Actual source code: test7.c

slepc-3.8.3 2018-04-03
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2017, 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       n=100,nev,its;
 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,&ksp);
 84:   KSPSetType(ksp[0],KSPBICG);
 85:   KSPGetPC(ksp[0],&pc);
 86:   PCSetType(pc,PCJACOBI);
 87:   KSPSetTolerances(ksp[0],tol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);

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

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

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

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

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

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

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

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

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

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

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

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

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

210: /* -------------------------------- A0 ----------------------------------- */

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

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

231: PetscErrorCode MatGetDiagonal_A0(Mat A,Vec diag)
232: {

236:   VecSet(diag,-2.0);
237:   return(0);
238: }

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

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

257: /* -------------------------------- A1 ----------------------------------- */

259: PetscErrorCode MatMult_A1(Mat A,Vec x,Vec y)
260: {

264:   VecCopy(x,y);
265:   return(0);
266: }

268: PetscErrorCode MatGetDiagonal_A1(Mat A,Vec diag)
269: {

273:   VecSet(diag,1.0);
274:   return(0);
275: }

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

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

294: /* -------------------------------- F ----------------------------------- */

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

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

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

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

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

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

351: PetscErrorCode MatDestroy_F(Mat A)
352: {
353:   MatCtx         *ctx;

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