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