Actual source code: test7.c
slepc-3.11.2 2019-07-30
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*/