eigen_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 #include "libmesh/libmesh_config.h"
20 
21 // Currently, the EigenSystem should only be available
22 // if SLEPc support is enabled.
23 #if defined(LIBMESH_HAVE_SLEPC)
24 
25 // C++ includes
26 
27 // Local includes
28 #include "libmesh/eigen_system.h"
30 #include "libmesh/sparse_matrix.h"
31 #include "libmesh/eigen_solver.h"
32 #include "libmesh/dof_map.h"
33 #include "libmesh/mesh_base.h"
34 
35 namespace libMesh
36 {
37 
38 
39 // ------------------------------------------------------------
40 // EigenSystem implementation
42  const std::string& name,
43  const unsigned int number
44  ) :
45  Parent (es, name, number),
46  matrix_A (NULL),
47  matrix_B (NULL),
48  eigen_solver (EigenSolver<Number>::build(es.comm())),
49  _n_converged_eigenpairs (0),
50  _n_iterations (0),
51  _is_generalized_eigenproblem (false),
52  _eigen_problem_type (NHEP)
53 {
54 }
55 
56 
57 
59 {
60  // clear data
61  this->clear();
62 }
63 
64 
65 
67 {
68  // Clear the parent data
69  Parent::clear();
70 
71  // delete the matricies
72  delete matrix_A;
73  delete matrix_B;
74 
75  // NULL-out the matricies.
76  matrix_A = NULL;
77  matrix_B = NULL;
78 
79  // clear the solver
80  eigen_solver->clear();
81 
82 }
83 
84 
86 {
87  _eigen_problem_type = ept;
88 
89  eigen_solver->set_eigenproblem_type(ept);
90 
91  // libMesh::out<< "The Problem type is set to be: "<<std::endl;
92 
93  switch (_eigen_problem_type)
94  {
95  case HEP: // libMesh::out<<"Hermitian"<<std::endl;
96  break;
97 
98  case NHEP: // libMesh::out<<"Non-Hermitian"<<std::endl;
99  break;
100 
101  case GHEP: // libMesh::out<<"Gerneralized Hermitian"<<std::endl;
102  break;
103 
104  case GNHEP: // libMesh::out<<"Generalized Non-Hermitian"<<std::endl;
105  break;
106 
107  default: // libMesh::out<<"not properly specified"<<std::endl;
108  libmesh_error();
109  break;
110 
111  }
112 }
113 
114 
115 
116 
118 {
119  // initialize parent data
121 
122  // define the type of eigenproblem
125 
126  // build the system matrix
127  matrix_A = SparseMatrix<Number>::build(this->comm()).release();
128 
129  this->init_matrices();
130 }
131 
132 
133 
135 {
136  DofMap& dof_map = this->get_dof_map();
137 
138  dof_map.attach_matrix(*matrix_A);
139 
140  // build matrix_B only in case of a
141  // generalized problem
143  {
144  matrix_B = SparseMatrix<Number>::build(this->comm()).release();
145  dof_map.attach_matrix(*matrix_B);
146  }
147 
148  dof_map.compute_sparsity(this->get_mesh());
149 
150  // initialize and zero system matrix
151  matrix_A->init();
152  matrix_A->zero();
153 
154  // eventually initialize and zero system matrix_B
156  {
157  matrix_B->init();
158  matrix_B->zero();
159  }
160 }
161 
162 
163 
165 {
166  // initialize parent data
167  Parent::reinit();
168 
169  // Clear the matrices
170  matrix_A->clear();
171 
173  matrix_B->clear();
174 
175  DofMap& dof_map = this->get_dof_map();
176 
177  // Clear the sparsity pattern
178  dof_map.clear_sparsity();
179 
180  // Compute the sparsity pattern for the current
181  // mesh and DOF distribution. This also updates
182  // both matrices, \p DofMap now knows them
183  dof_map.compute_sparsity(this->get_mesh());
184 
185  matrix_A->init();
186  matrix_A->zero();
187 
189  {
190  matrix_B->init();
191  matrix_B->zero();
192  }
193 }
194 
195 
196 
198 {
199 
200  // A reference to the EquationSystems
202 
203  // check that necessary parameters have been set
204  libmesh_assert (es.parameters.have_parameter<unsigned int>("eigenpairs"));
205  libmesh_assert (es.parameters.have_parameter<unsigned int>("basis vectors"));
206 
207  if (this->assemble_before_solve)
208  // Assemble the linear system
209  this->assemble ();
210 
211  // Get the tolerance for the solver and the maximum
212  // number of iterations. Here, we simply adopt the linear solver
213  // specific parameters.
214  const Real tol =
215  es.parameters.get<Real>("linear solver tolerance");
216 
217  const unsigned int maxits =
218  es.parameters.get<unsigned int>("linear solver maximum iterations");
219 
220  const unsigned int nev =
221  es.parameters.get<unsigned int>("eigenpairs");
222 
223  const unsigned int ncv =
224  es.parameters.get<unsigned int>("basis vectors");
225 
226  std::pair<unsigned int, unsigned int> solve_data;
227 
228  // call the solver depending on the type of eigenproblem
230  {
231  //in case of a generalized eigenproblem
232  solve_data = eigen_solver->solve_generalized (*matrix_A,*matrix_B, nev, ncv, tol, maxits);
233 
234  }
235 
236  else
237  {
239 
240  //in case of a standard eigenproblem
241  solve_data = eigen_solver->solve_standard (*matrix_A, nev, ncv, tol, maxits);
242 
243  }
244 
245  this->_n_converged_eigenpairs = solve_data.first;
246  this->_n_iterations = solve_data.second;
247 
248 }
249 
250 
252 {
253 
254  // Assemble the linear system
255  Parent::assemble ();
256 
257 }
258 
259 
260 std::pair<Real, Real> EigenSystem::get_eigenpair (unsigned int i)
261 {
262  // call the eigen_solver get_eigenpair method
263  return eigen_solver->get_eigenpair (i, *solution);
264 }
265 
266 } // namespace libMesh
267 
268 #endif // LIBMESH_HAVE_SLEPC

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo