Actual source code: loaded_string.c
slepc-3.6.1 2015-09-03
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2015, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
21: /*
22: This example implements one of the problems found at
23: NLEVP: A Collection of Nonlinear Eigenvalue Problems,
24: The University of Manchester.
25: The details of the collection can be found at:
26: [1] T. Betcke et al., "NLEVP: A Collection of Nonlinear Eigenvalue
27: Problems", ACM Trans. Math. Software 39(2), Article 7, 2013.
29: The loaded_string problem is a rational eigenvalue problem for the
30: finite element model of a loaded vibrating string.
31: */
33: static char help[] = "NLEVP problem: loaded_string.\n\n"
34: "The command line options are:\n"
35: " -n <n>, dimension of the matrices.\n"
36: " -kappa <kappa>, stiffness of elastic spring.\n"
37: " -mass <m>, mass of the attached load.\n\n";
39: #include <slepcnep.h>
41: #define NMAT 3
45: int main(int argc,char **argv)
46: {
47: Mat A[NMAT]; /* problem matrices */
48: FN f[NMAT]; /* functions to define the nonlinear operator */
49: NEP nep; /* nonlinear eigensolver context */
50: PetscInt n=20,Istart,Iend,i;
51: PetscReal kappa=1.0,m=1.0;
52: PetscScalar sigma,numer[2],denom[2];
53: PetscBool terse;
56: SlepcInitialize(&argc,&argv,(char*)0,help);
58: PetscOptionsGetInt(NULL,"-n",&n,NULL);
59: PetscOptionsGetReal(NULL,"-kappa",&kappa,NULL);
60: PetscOptionsGetReal(NULL,"-mass",&m,NULL);
61: sigma = kappa/m;
62: PetscPrintf(PETSC_COMM_WORLD,"Loaded vibrating string, n=%D kappa=%g m=%g\n\n",n,(double)kappa,(double)m);
64: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65: Build the problem matrices
66: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68: /* initialize matrices */
69: for (i=0;i<NMAT;i++) {
70: MatCreate(PETSC_COMM_WORLD,&A[i]);
71: MatSetSizes(A[i],PETSC_DECIDE,PETSC_DECIDE,n,n);
72: MatSetFromOptions(A[i]);
73: MatSetUp(A[i]);
74: }
75: MatGetOwnershipRange(A[0],&Istart,&Iend);
77: /* A0 */
78: for (i=Istart;i<Iend;i++) {
79: MatSetValue(A[0],i,i,(i==n-1)?1.0*n:2.0*n,INSERT_VALUES);
80: if (i>0) { MatSetValue(A[0],i,i-1,-1.0*n,INSERT_VALUES); }
81: if (i<n-1) { MatSetValue(A[0],i,i+1,-1.0*n,INSERT_VALUES); }
82: }
84: /* A1 */
86: for (i=Istart;i<Iend;i++) {
87: MatSetValue(A[1],i,i,(i==n-1)?2.0/(6.0*n):4.0/(6.0*n),INSERT_VALUES);
88: if (i>0) { MatSetValue(A[1],i,i-1,1.0/(6.0*n),INSERT_VALUES); }
89: if (i<n-1) { MatSetValue(A[1],i,i+1,1.0/(6.0*n),INSERT_VALUES); }
90: }
92: /* A2 */
93: if (Istart<=n-1 && n-1<Iend) {
94: MatSetValue(A[2],n-1,n-1,kappa,INSERT_VALUES);
95: }
97: /* assemble matrices */
98: for (i=0;i<NMAT;i++) {
99: MatAssemblyBegin(A[i],MAT_FINAL_ASSEMBLY);
100: }
101: for (i=0;i<NMAT;i++) {
102: MatAssemblyEnd(A[i],MAT_FINAL_ASSEMBLY);
103: }
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Create the problem functions
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: /* f1=1 */
110: FNCreate(PETSC_COMM_WORLD,&f[0]);
111: FNSetType(f[0],FNRATIONAL);
112: numer[0] = 1.0;
113: FNRationalSetNumerator(f[0],1,numer);
115: /* f2=-lambda */
116: FNCreate(PETSC_COMM_WORLD,&f[1]);
117: FNSetType(f[1],FNRATIONAL);
118: numer[0] = -1.0; numer[1] = 0.0;
119: FNRationalSetNumerator(f[1],2,numer);
121: /* f3=lambda/(lambda-sigma) */
122: FNCreate(PETSC_COMM_WORLD,&f[2]);
123: FNSetType(f[2],FNRATIONAL);
124: numer[0] = 1.0; numer[1] = 0.0;
125: denom[0] = 1.0; denom[1] = -sigma;
126: FNRationalSetNumerator(f[2],2,numer);
127: FNRationalSetDenominator(f[2],2,denom);
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Create the eigensolver and solve the problem
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: NEPCreate(PETSC_COMM_WORLD,&nep);
134: NEPSetSplitOperator(nep,3,A,f,SUBSET_NONZERO_PATTERN);
135: NEPSetFromOptions(nep);
136: NEPSolve(nep);
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Display solution and clean up
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141:
142: /* show detailed info unless -terse option is given by user */
143: PetscOptionsHasName(NULL,"-terse",&terse);
144: if (terse) {
145: NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
146: } else {
147: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
148: NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
149: NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
150: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
151: }
152: NEPDestroy(&nep);
153: for (i=0;i<NMAT;i++) {
154: MatDestroy(&A[i]);
155: FNDestroy(&f[i]);
156: }
157: SlepcFinalize();
158: return 0;
159: }