Actual source code: test2.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2013, 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 QEP 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 <slepcqep.h>
43: int main(int argc,char **argv)
44: {
45: Mat M,C,K; /* problem matrices */
46: QEP qep; /* quadratic eigenproblem solver context */
47: QEPType 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: QEPCreate(PETSC_COMM_WORLD,&qep);
124: /*
125: Set matrices, the problem type and other settings
126: */
127: QEPSetOperators(qep,M,C,K);
128: QEPSetProblemType(qep,QEP_GENERAL);
129: QEPSetTolerances(qep,PETSC_SMALL,PETSC_DEFAULT);
130: QEPSetFromOptions(qep);
132: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: Solve the eigensystem
134: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136: QEPSolve(qep);
138: /*
139: Optional: Get some information from the solver and display it
140: */
141: QEPGetType(qep,&type);
142: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
143: QEPGetDimensions(qep,&nev,NULL,NULL);
144: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
145: QEPGetTolerances(qep,NULL,&maxit);
146: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: maxit=%D\n",maxit);
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Display solution and clean up
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: QEPPrintSolution(qep,NULL);
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Free work space
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158: QEPDestroy(&qep);
159: MatDestroy(&M);
160: MatDestroy(&C);
161: MatDestroy(&K);
162: SlepcFinalize();
163: return 0;
164: }