Actual source code: test2.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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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 paramter.\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,"-n",&n,NULL);
55: PetscOptionsGetScalar(NULL,"-mu",&mu,NULL);
56: PetscOptionsGetScalar(NULL,"-tau",&tau,NULL);
57: PetscOptionsGetScalar(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: /*
120: Create eigensolver context
121: */
122: PEPCreate(PETSC_COMM_WORLD,&pep);
124: /*
125: Set matrices, the problem type and other settings
126: */
127: A[0] = K; A[1] = C; A[2] = M;
128: PEPSetOperators(pep,3,A);
129: PEPSetProblemType(pep,PEP_GENERAL);
130: PEPSetTolerances(pep,PETSC_SMALL,PETSC_DEFAULT);
131: PEPSetFromOptions(pep);
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: Solve the eigensystem
135: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: PEPSolve(pep);
139: /*
140: Optional: Get some information from the solver and display it
141: */
142: PEPGetType(pep,&type);
143: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
144: PEPGetDimensions(pep,&nev,NULL,NULL);
145: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
146: PEPGetTolerances(pep,NULL,&maxit);
147: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: maxit=%D\n",maxit);
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Display solution and clean up
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
154: PEPDestroy(&pep);
155: MatDestroy(&M);
156: MatDestroy(&C);
157: MatDestroy(&K);
158: SlepcFinalize();
159: return 0;
160: }