Actual source code: test7.c
slepc-3.8.2 2017-12-01
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: }