Actual source code: test7.c
slepc-3.17.2 2022-08-09
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*/