Actual source code: test2.c
slepc-3.7.4 2017-05-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2016, 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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
21: Example based on spring problem in NLEVP collection [1]. See the parameters
22: meaning at Example 2 in [2].
24: [1] T. Betcke, N. J. Higham, V. Mehrmann, C. Schroder, and F. Tisseur,
25: NLEVP: A Collection of Nonlinear Eigenvalue Problems, MIMS EPrint
26: 2010.98, November 2010.
27: [2] F. Tisseur, Backward error and condition of polynomial eigenvalue
28: problems, Linear Algebra and its Applications, 309 (2000), pp. 339--361,
29: April 2000.
30: */
32: static char help[] = "Test the solution of a PEP from a finite element model of "
33: "damped mass-spring system (problem from NLEVP collection).\n\n"
34: "The command line options are:\n"
35: " -n <n>, where <n> = number of grid subdivisions.\n"
36: " -tau <tau>, where <tau> = tau parameter.\n"
37: " -kappa <kappa>, where <kappa> = kappa parameter.\n\n";
39: #include <slepcpep.h>
43: int main(int argc,char **argv)
44: {
45: Mat M,C,K,A[3]; /* problem matrices */
46: PEP pep; /* polynomial eigenproblem solver context */
47: PEPType type;
49: PetscInt n=30,Istart,Iend,i,maxit,nev;
50: PetscScalar mu=1.0,tau=10.0,kappa=5.0;
52: SlepcInitialize(&argc,&argv,(char*)0,help);
54: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
55: PetscOptionsGetScalar(NULL,NULL,"-mu",&mu,NULL);
56: PetscOptionsGetScalar(NULL,NULL,"-tau",&tau,NULL);
57: PetscOptionsGetScalar(NULL,NULL,"-kappa",&kappa,NULL);
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: /* K is a tridiagonal */
64: MatCreate(PETSC_COMM_WORLD,&K);
65: MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n);
66: MatSetFromOptions(K);
67: MatSetUp(K);
69: MatGetOwnershipRange(K,&Istart,&Iend);
70: for (i=Istart;i<Iend;i++) {
71: if (i>0) {
72: MatSetValue(K,i,i-1,-kappa,INSERT_VALUES);
73: }
74: MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES);
75: if (i<n-1) {
76: MatSetValue(K,i,i+1,-kappa,INSERT_VALUES);
77: }
78: }
80: MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
81: MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);
83: /* C is a tridiagonal */
84: MatCreate(PETSC_COMM_WORLD,&C);
85: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
86: MatSetFromOptions(C);
87: MatSetUp(C);
89: MatGetOwnershipRange(C,&Istart,&Iend);
90: for (i=Istart;i<Iend;i++) {
91: if (i>0) {
92: MatSetValue(C,i,i-1,-tau,INSERT_VALUES);
93: }
94: MatSetValue(C,i,i,tau*3.0,INSERT_VALUES);
95: if (i<n-1) {
96: MatSetValue(C,i,i+1,-tau,INSERT_VALUES);
97: }
98: }
100: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
101: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
103: /* M is a diagonal matrix */
104: MatCreate(PETSC_COMM_WORLD,&M);
105: MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n);
106: MatSetFromOptions(M);
107: MatSetUp(M);
108: MatGetOwnershipRange(M,&Istart,&Iend);
109: for (i=Istart;i<Iend;i++) {
110: MatSetValue(M,i,i,mu,INSERT_VALUES);
111: }
112: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
113: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Create the eigensolver and set various options
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: PEPCreate(PETSC_COMM_WORLD,&pep);
120: A[0] = K; A[1] = C; A[2] = M;
121: PEPSetOperators(pep,3,A);
122: PEPSetProblemType(pep,PEP_GENERAL);
123: PEPSetTolerances(pep,PETSC_SMALL,PETSC_DEFAULT);
124: PEPSetFromOptions(pep);
126: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: Solve the eigensystem
128: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: PEPSolve(pep);
132: /*
133: Optional: Get some information from the solver and display it
134: */
135: PEPGetType(pep,&type);
136: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
137: PEPGetDimensions(pep,&nev,NULL,NULL);
138: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
139: PEPGetTolerances(pep,NULL,&maxit);
140: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: maxit=%D\n",maxit);
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Display solution and clean up
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146: PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
147: PEPDestroy(&pep);
148: MatDestroy(&M);
149: MatDestroy(&C);
150: MatDestroy(&K);
151: SlepcFinalize();
152: return ierr;
153: }