Actual source code: test7.c
slepc-3.21.1 2024-04-26
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: PetscFunctionBeginUser;
60: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
61: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
62: PetscCall(PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL));
63: PetscCall(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: PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
70: PetscCall(NEPSetType(nep,NEPNLEIGS));
71: PetscCall(NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,NULL));
72: PetscCall(NEPGetRG(nep,&rg));
73: PetscCall(RGSetType(rg,RGINTERVAL));
74: #if defined(PETSC_USE_COMPLEX)
75: PetscCall(RGIntervalSetEndpoints(rg,0.01,16.0,-0.001,0.001));
76: #else
77: PetscCall(RGIntervalSetEndpoints(rg,0.01,16.0,0,0));
78: #endif
79: PetscCall(NEPSetTarget(nep,1.1));
80: PetscCall(NEPNLEIGSGetKSPs(nep,&nsolve,&ksp));
81: for (i=0;i<nsolve;i++) {
82: PetscCall(KSPSetType(ksp[i],KSPBICG));
83: PetscCall(KSPGetPC(ksp[i],&pc));
84: PetscCall(PCSetType(pc,PCJACOBI));
85: PetscCall(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: PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,NULL,&A[0]));
95: PetscCall(MatShellSetOperation(A[0],MATOP_MULT,(void(*)(void))MatMult_A0));
96: PetscCall(MatShellSetOperation(A[0],MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A0));
97: PetscCall(MatShellSetOperation(A[0],MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A0));
98: PetscCall(MatShellSetOperation(A[0],MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A0));
100: /* Create matrix A0 (identity) */
101: PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,NULL,&A[1]));
102: PetscCall(MatShellSetOperation(A[1],MATOP_MULT,(void(*)(void))MatMult_A1));
103: PetscCall(MatShellSetOperation(A[1],MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A1));
104: PetscCall(MatShellSetOperation(A[1],MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A1));
105: PetscCall(MatShellSetOperation(A[1],MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A1));
107: /* Define functions for the split form */
108: PetscCall(FNCreate(PETSC_COMM_WORLD,&f[0]));
109: PetscCall(FNSetType(f[0],FNRATIONAL));
110: coeffs = 1.0;
111: PetscCall(FNRationalSetNumerator(f[0],1,&coeffs));
112: PetscCall(FNCreate(PETSC_COMM_WORLD,&f[1]));
113: PetscCall(FNSetType(f[1],FNSQRT));
114: PetscCall(NEPSetSplitOperator(nep,2,A,f,SUBSET_NONZERO_PATTERN));
115: } else {
116: /* Callback form: create shell matrix for F=A0+sqrt(lambda)*A1 */
117: PetscCall(PetscNew(&ctx));
118: PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctx,&F));
119: PetscCall(MatShellSetOperation(F,MATOP_MULT,(void(*)(void))MatMult_F));
120: PetscCall(MatShellSetOperation(F,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_F));
121: PetscCall(MatShellSetOperation(F,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_F));
122: PetscCall(MatShellSetOperation(F,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_F));
123: PetscCall(MatShellSetOperation(F,MATOP_DESTROY,(void(*)(void))MatDestroy_F));
124: /* Set Function evaluation routine */
125: PetscCall(NEPSetFunction(nep,F,F,FormFunction,NULL));
126: }
128: /* Set solver parameters at runtime */
129: PetscCall(NEPSetFromOptions(nep));
131: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: Solve the eigensystem
133: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: PetscCall(NEPSolve(nep));
135: PetscCall(NEPGetType(nep,&type));
136: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type));
137: PetscCall(NEPGetDimensions(nep,&nev,NULL,NULL));
138: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
139: PetscCall(PetscObjectTypeCompare((PetscObject)nep,NEPNLEIGS,&flg));
140: if (flg) {
141: PetscCall(NEPNLEIGSGetRestart(nep,&keep));
142: PetscCall(NEPNLEIGSGetLocking(nep,&lock));
143: PetscCall(NEPNLEIGSGetInterpolation(nep,&tol,&its));
144: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Restart factor is %3.2f",(double)keep));
145: if (lock) PetscCall(PetscPrintf(PETSC_COMM_WORLD," (locking activated)"));
146: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n Divided diferences with tol=%6.2g maxit=%" PetscInt_FMT "\n",(double)tol,its));
147: PetscCall(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: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
156: if (terse) PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
157: else {
158: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
159: PetscCall(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD));
160: PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
161: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
162: }
163: PetscCall(NEPDestroy(&nep));
164: if (split) {
165: PetscCall(MatDestroy(&A[0]));
166: PetscCall(MatDestroy(&A[1]));
167: PetscCall(FNDestroy(&f[0]));
168: PetscCall(FNDestroy(&f[1]));
169: } else PetscCall(MatDestroy(&F));
170: PetscCall(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;
181: PetscFunctionBeginUser;
182: PetscCall(MatShellGetContext(fun,&ctxF));
183: ctxF->t = PetscSqrtScalar(lambda);
184: PetscFunctionReturn(PETSC_SUCCESS);
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;
198: PetscFunctionBeginUser;
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: PetscFunctionReturn(PETSC_SUCCESS);
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;
215: PetscFunctionBeginUser;
216: PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
217: PetscCallMPI(MPI_Comm_size(comm,&size));
218: PetscCallMPI(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: PetscCall(VecGetArrayRead(x,&px));
223: PetscCall(VecGetArray(y,&py));
224: PetscCall(VecGetLocalSize(x,&n));
226: PetscCallMPI(MPI_Sendrecv(px,1,MPIU_SCALAR,prev,0,&lower,1,MPIU_SCALAR,next,0,comm,MPI_STATUS_IGNORE));
227: PetscCallMPI(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: PetscCall(VecRestoreArrayRead(x,&px));
233: PetscCall(VecRestoreArray(y,&py));
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: PetscErrorCode MatGetDiagonal_A0(Mat A,Vec diag)
238: {
239: PetscFunctionBeginUser;
240: PetscCall(VecSet(diag,-2.0));
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: PetscErrorCode MatDuplicate_A0(Mat A,MatDuplicateOption op,Mat *B)
245: {
246: PetscInt m,n,M,N;
247: MPI_Comm comm;
249: PetscFunctionBegin;
250: PetscCall(MatGetSize(A,&M,&N));
251: PetscCall(MatGetLocalSize(A,&m,&n));
252: PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
253: PetscCall(MatCreateShell(comm,m,n,M,N,NULL,B));
254: PetscCall(MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_A0));
255: PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A0));
256: PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A0));
257: PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A0));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: /* -------------------------------- A1 ----------------------------------- */
263: PetscErrorCode MatMult_A1(Mat A,Vec x,Vec y)
264: {
265: PetscFunctionBeginUser;
266: PetscCall(VecCopy(x,y));
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: PetscErrorCode MatGetDiagonal_A1(Mat A,Vec diag)
271: {
272: PetscFunctionBeginUser;
273: PetscCall(VecSet(diag,1.0));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: PetscErrorCode MatDuplicate_A1(Mat A,MatDuplicateOption op,Mat *B)
278: {
279: PetscInt m,n,M,N;
280: MPI_Comm comm;
282: PetscFunctionBegin;
283: PetscCall(MatGetSize(A,&M,&N));
284: PetscCall(MatGetLocalSize(A,&m,&n));
285: PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
286: PetscCall(MatCreateShell(comm,m,n,M,N,NULL,B));
287: PetscCall(MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_A1));
288: PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_A1));
289: PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_A1));
290: PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_A1));
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
294: /* -------------------------------- F ----------------------------------- */
296: PetscErrorCode MatMult_F(Mat A,Vec x,Vec y)
297: {
298: PetscInt i,n;
299: PetscMPIInt rank,size,next,prev;
300: const PetscScalar *px;
301: PetscScalar *py,d,upper=0.0,lower=0.0;
302: MatCtx *ctx;
303: MPI_Comm comm;
305: PetscFunctionBeginUser;
306: PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
307: PetscCallMPI(MPI_Comm_size(comm,&size));
308: PetscCallMPI(MPI_Comm_rank(comm,&rank));
309: next = rank==size-1? MPI_PROC_NULL: rank+1;
310: prev = rank==0? MPI_PROC_NULL: rank-1;
312: PetscCall(MatShellGetContext(A,&ctx));
313: PetscCall(VecGetArrayRead(x,&px));
314: PetscCall(VecGetArray(y,&py));
315: PetscCall(VecGetLocalSize(x,&n));
317: PetscCallMPI(MPI_Sendrecv(px,1,MPIU_SCALAR,prev,0,&lower,1,MPIU_SCALAR,next,0,comm,MPI_STATUS_IGNORE));
318: PetscCallMPI(MPI_Sendrecv(px+n-1,1,MPIU_SCALAR,next,0,&upper,1,MPIU_SCALAR,prev,0,comm,MPI_STATUS_IGNORE));
320: d = -2.0+ctx->t;
321: py[0] = upper+d*px[0]+px[1];
322: for (i=1;i<n-1;i++) py[i] = px[i-1]+d*px[i]+px[i+1];
323: py[n-1] = px[n-2]+d*px[n-1]+lower;
324: PetscCall(VecRestoreArrayRead(x,&px));
325: PetscCall(VecRestoreArray(y,&py));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: PetscErrorCode MatGetDiagonal_F(Mat A,Vec diag)
330: {
331: MatCtx *ctx;
333: PetscFunctionBeginUser;
334: PetscCall(MatShellGetContext(A,&ctx));
335: PetscCall(VecSet(diag,-2.0+ctx->t));
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: PetscErrorCode MatDuplicate_F(Mat A,MatDuplicateOption op,Mat *B)
340: {
341: MatCtx *actx,*bctx;
342: PetscInt m,n,M,N;
343: MPI_Comm comm;
345: PetscFunctionBegin;
346: PetscCall(MatShellGetContext(A,&actx));
347: PetscCall(MatGetSize(A,&M,&N));
348: PetscCall(MatGetLocalSize(A,&m,&n));
349: PetscCall(PetscNew(&bctx));
350: bctx->t = actx->t;
351: PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
352: PetscCall(MatCreateShell(comm,m,n,M,N,(void*)bctx,B));
353: PetscCall(MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_F));
354: PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMult_F));
355: PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_F));
356: PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_F));
357: PetscCall(MatShellSetOperation(*B,MATOP_DESTROY,(void(*)(void))MatDestroy_F));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: PetscErrorCode MatDestroy_F(Mat A)
362: {
363: MatCtx *ctx;
365: PetscFunctionBegin;
366: PetscCall(MatShellGetContext(A,&ctx));
367: PetscCall(PetscFree(ctx));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: /*TEST
373: testset:
374: nsize: {{1 2}}
375: args: -nep_nev 3 -nep_tol 1e-8 -terse
376: filter: sed -e "s/[+-]0\.0*i//g"
377: requires: !single
378: test:
379: suffix: 1
380: args: -nep_nleigs_locking 0 -nep_nleigs_interpolation_degree 90 -nep_nleigs_interpolation_tol 1e-8 -nep_nleigs_restart 0.4
381: test:
382: suffix: 2
383: args: -split 0
385: TEST*/