Actual source code: ex42.c
slepc-3.17.1 2022-04-11
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: */
10: /*
11: This example implements one of the problems found at
12: NLEVP: A Collection of Nonlinear Eigenvalue Problems,
13: The University of Manchester.
14: The details of the collection can be found at:
15: [1] T. Betcke et al., "NLEVP: A Collection of Nonlinear Eigenvalue
16: Problems", ACM Trans. Math. Software 39(2), Article 7, 2013.
18: The loaded_string problem is a rational eigenvalue problem for the
19: finite element model of a loaded vibrating string.
20: */
22: static char help[] = "Illustrates computation of left eigenvectors and resolvent.\n\n"
23: "This is based on loaded_string from the NLEVP collection.\n"
24: "The command line options are:\n"
25: " -n <n>, dimension of the matrices.\n"
26: " -kappa <kappa>, stiffness of elastic spring.\n"
27: " -mass <m>, mass of the attached load.\n\n";
29: #include <slepcnep.h>
31: #define NMAT 3
33: int main(int argc,char **argv)
34: {
35: Mat A[NMAT]; /* problem matrices */
36: FN f[NMAT]; /* functions to define the nonlinear operator */
37: NEP nep; /* nonlinear eigensolver context */
38: RG rg;
39: Vec v,r,z,w;
40: PetscInt n=100,Istart,Iend,i,nconv;
41: PetscReal kappa=1.0,m=1.0,nrm,tol;
42: PetscScalar lambda,sigma,numer[2],denom[2],omega1,omega2;
44: SlepcInitialize(&argc,&argv,(char*)0,help);
46: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
47: PetscOptionsGetReal(NULL,NULL,"-kappa",&kappa,NULL);
48: PetscOptionsGetReal(NULL,NULL,"-mass",&m,NULL);
49: sigma = kappa/m;
50: PetscPrintf(PETSC_COMM_WORLD,"Loaded vibrating string, n=%" PetscInt_FMT " kappa=%g m=%g\n\n",n,(double)kappa,(double)m);
52: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: Build the problem matrices
54: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56: /* initialize matrices */
57: for (i=0;i<NMAT;i++) {
58: MatCreate(PETSC_COMM_WORLD,&A[i]);
59: MatSetSizes(A[i],PETSC_DECIDE,PETSC_DECIDE,n,n);
60: MatSetFromOptions(A[i]);
61: MatSetUp(A[i]);
62: }
63: MatGetOwnershipRange(A[0],&Istart,&Iend);
65: /* A0 */
66: for (i=Istart;i<Iend;i++) {
67: MatSetValue(A[0],i,i,(i==n-1)?1.0*n:2.0*n,INSERT_VALUES);
68: if (i>0) MatSetValue(A[0],i,i-1,-1.0*n,INSERT_VALUES);
69: if (i<n-1) MatSetValue(A[0],i,i+1,-1.0*n,INSERT_VALUES);
70: }
72: /* A1 */
73: for (i=Istart;i<Iend;i++) {
74: MatSetValue(A[1],i,i,(i==n-1)?2.0/(6.0*n):4.0/(6.0*n),INSERT_VALUES);
75: if (i>0) MatSetValue(A[1],i,i-1,1.0/(6.0*n),INSERT_VALUES);
76: if (i<n-1) MatSetValue(A[1],i,i+1,1.0/(6.0*n),INSERT_VALUES);
77: }
79: /* A2 */
80: if (Istart<=n-1 && n-1<Iend) MatSetValue(A[2],n-1,n-1,kappa,INSERT_VALUES);
82: /* assemble matrices */
83: for (i=0;i<NMAT;i++) MatAssemblyBegin(A[i],MAT_FINAL_ASSEMBLY);
84: for (i=0;i<NMAT;i++) MatAssemblyEnd(A[i],MAT_FINAL_ASSEMBLY);
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Create the problem functions
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: /* f1=1 */
91: FNCreate(PETSC_COMM_WORLD,&f[0]);
92: FNSetType(f[0],FNRATIONAL);
93: numer[0] = 1.0;
94: FNRationalSetNumerator(f[0],1,numer);
96: /* f2=-lambda */
97: FNCreate(PETSC_COMM_WORLD,&f[1]);
98: FNSetType(f[1],FNRATIONAL);
99: numer[0] = -1.0; numer[1] = 0.0;
100: FNRationalSetNumerator(f[1],2,numer);
102: /* f3=lambda/(lambda-sigma) */
103: FNCreate(PETSC_COMM_WORLD,&f[2]);
104: FNSetType(f[2],FNRATIONAL);
105: numer[0] = 1.0; numer[1] = 0.0;
106: denom[0] = 1.0; denom[1] = -sigma;
107: FNRationalSetNumerator(f[2],2,numer);
108: FNRationalSetDenominator(f[2],2,denom);
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Create the eigensolver and solve the problem
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: NEPCreate(PETSC_COMM_WORLD,&nep);
115: NEPSetSplitOperator(nep,3,A,f,SUBSET_NONZERO_PATTERN);
116: NEPSetProblemType(nep,NEP_RATIONAL);
117: NEPSetDimensions(nep,8,PETSC_DEFAULT,PETSC_DEFAULT);
119: /* set two-sided NLEIGS solver */
120: NEPSetType(nep,NEPNLEIGS);
121: NEPNLEIGSSetFullBasis(nep,PETSC_TRUE);
122: NEPSetTwoSided(nep,PETSC_TRUE);
123: NEPGetRG(nep,&rg);
124: RGSetType(rg,RGINTERVAL);
125: #if defined(PETSC_USE_COMPLEX)
126: RGIntervalSetEndpoints(rg,4.0,700.0,-0.001,0.001);
127: #else
128: RGIntervalSetEndpoints(rg,4.0,700.0,0,0);
129: #endif
130: NEPSetTarget(nep,5.0);
132: NEPSetFromOptions(nep);
133: NEPSolve(nep);
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Check left residual
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: MatCreateVecs(A[0],&v,&r);
139: VecDuplicate(v,&w);
140: VecDuplicate(v,&z);
141: NEPGetConverged(nep,&nconv);
142: NEPGetTolerances(nep,&tol,NULL);
143: for (i=0;i<nconv;i++) {
144: NEPGetEigenpair(nep,i,&lambda,NULL,NULL,NULL);
145: NEPGetLeftEigenvector(nep,i,v,NULL);
146: NEPApplyAdjoint(nep,lambda,v,w,r,NULL,NULL);
147: VecNorm(r,NORM_2,&nrm);
148: if (nrm>tol*PetscAbsScalar(lambda)) PetscPrintf(PETSC_COMM_WORLD,"Left residual i=%" PetscInt_FMT " is above tolerance --> %g\n",i,(double)(nrm/PetscAbsScalar(lambda)));
149: }
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: Operate with resolvent
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: omega1 = 20.0;
155: omega2 = 150.0;
156: VecSet(v,0.0);
157: VecSetValue(v,0,-1.0,INSERT_VALUES);
158: VecSetValue(v,1,3.0,INSERT_VALUES);
159: VecAssemblyBegin(v);
160: VecAssemblyEnd(v);
161: NEPApplyResolvent(nep,NULL,omega1,v,r);
162: VecNorm(r,NORM_2,&nrm);
163: PetscPrintf(PETSC_COMM_WORLD,"resolvent, omega=%g: norm of computed vector=%g\n",(double)PetscRealPart(omega1),(double)nrm);
164: NEPApplyResolvent(nep,NULL,omega2,v,r);
165: VecNorm(r,NORM_2,&nrm);
166: PetscPrintf(PETSC_COMM_WORLD,"resolvent, omega=%g: norm of computed vector=%g\n",(double)PetscRealPart(omega2),(double)nrm);
167: VecSet(v,1.0);
168: NEPApplyResolvent(nep,NULL,omega1,v,r);
169: VecNorm(r,NORM_2,&nrm);
170: PetscPrintf(PETSC_COMM_WORLD,"resolvent, omega=%g: norm of computed vector=%g\n",(double)PetscRealPart(omega1),(double)nrm);
171: NEPApplyResolvent(nep,NULL,omega2,v,r);
172: VecNorm(r,NORM_2,&nrm);
173: PetscPrintf(PETSC_COMM_WORLD,"resolvent, omega=%g: norm of computed vector=%g\n",(double)PetscRealPart(omega2),(double)nrm);
175: /* clean up */
176: NEPDestroy(&nep);
177: for (i=0;i<NMAT;i++) {
178: MatDestroy(&A[i]);
179: FNDestroy(&f[i]);
180: }
181: VecDestroy(&v);
182: VecDestroy(&r);
183: VecDestroy(&w);
184: VecDestroy(&z);
185: SlepcFinalize();
186: return 0;
187: }
189: /*TEST
191: test:
192: suffix: 1
193: requires: !single
195: TEST*/