ug4
power_method.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2018: G-CSC, Goethe University Frankfurt
3  * Author: Martin Stepniewski
4  *
5  * This file is part of UG4.
6  *
7  * UG4 is free software: you can redistribute it and/or modify it under the
8  * terms of the GNU Lesser General Public License version 3 (as published by the
9  * Free Software Foundation) with the following additional attribution
10  * requirements (according to LGPL/GPL v3 §7):
11  *
12  * (1) The following notice must be displayed in the Appropriate Legal Notices
13  * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14  *
15  * (2) The following notice must be displayed at a prominent place in the
16  * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17  *
18  * (3) The following bibliography is recommended for citation and must be
19  * preserved in all covered files:
20  * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21  * parallel geometric multigrid solver on hierarchically distributed grids.
22  * Computing and visualization in science 16, 4 (2013), 151-164"
23  * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24  * flexible software system for simulating pde based models on high performance
25  * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26  *
27  * This program is distributed in the hope that it will be useful,
28  * but WITHOUT ANY WARRANTY; without even the implied warranty of
29  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  * GNU Lesser General Public License for more details.
31  */
32 
33 #ifndef __H__UG__LIB_ALGEBRA__POWER_METHOD_H__
34 #define __H__UG__LIB_ALGEBRA__POWER_METHOD_H__
35 
36 #include <vector>
37 #include <string>
38 #include "additional_math.h"
39 
44 
47 
48 
49 namespace ug{
50 
51 
59 template<typename TAlgebra>
61 {
62  public:
63  // Algebra type
64  typedef TAlgebra algebra_type;
65 
66  // Vector type
67  typedef typename TAlgebra::vector_type vector_type;
68  typedef typename TAlgebra::matrix_type matrix_type;
70 
71  private:
72  // Stiffness matrix
74 
75  // Mass matrix
80  //matrix_operator_type *m_pB;
81 
82  // Eigenvector and -value
87 
88  // Solver
90 
91  // Precision
93  double m_dPrecision;
94 
95  // Residual
97 
98  // Dirichlet node indexing
99  std::vector<bool> m_vbDirichlet;
101 
102  public:
103  // Number of iterations
104  size_t m_iteration;
105 
107  {
108  UG_LOG("Initializing PowerMethod." << std::endl);
109  m_spLinOpA = SPNULL;
110  m_spLinOpB = SPNULL;
111  m_spMatOpA = SPNULL;
112  m_spMatOpB = SPNULL;
113  m_maxIterations = 1000;
114  m_dPrecision = 1e-8;
115  m_iteration = 0;
116  m_dMaxEigenvalue = 0.0;
117  m_dMinEigenvalue = 0.0;
118  m_numDirichletRows = 0;
119  }
120 
122  {
123  m_spSolver = solver;
124  }
125 
127  {
128  m_spLinOpA = loA;
129 
130  // get dirichlet nodes
131  m_spMatOpA = m_spLinOpA.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
132  matrix_type& A = m_spMatOpA->get_matrix();
133  m_vbDirichlet.resize(A.num_rows());
134 
135  for(size_t i=0; i<A.num_rows(); i++)
136  {
137  m_vbDirichlet[i] = A.is_isolated(i);
138 
139  if(A.is_isolated(i) == 0)
141  }
142  }
143 
145  {
146  m_spLinOpB = loB;
147  m_spMatOpB = m_spLinOpB.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
148  }
149 
151  {
152  m_spEigenvector = vec;
153  }
154 
155  void set_max_iterations(size_t maxIterations)
156  {
157  m_maxIterations = maxIterations;
158  }
159 
160  void set_precision(double precision)
161  {
162  m_dPrecision = precision;
163  }
164 
166  {
167  PROFILE_FUNC_GROUP("PowerMethod");
168 
169  UG_COND_THROW(m_spSolver == SPNULL && m_spLinOpB != SPNULL, "PowerMethod::calculate_max_eigenvalue(): Solver not set, please specify.");
170  if(m_spLinOpB != SPNULL)
171  m_spSolver->init(m_spLinOpB);
172 
174 
176  {
178 
179  // residual = A v
180  UG_COND_THROW(m_spLinOpA == SPNULL, "PowerMethod::calculate_max_eigenvalue(): Linear operator A not set, please specify.");
182 
183  // v = B^-1 A v
184  if(m_spLinOpB != SPNULL)
186  // in case B is not set: v = A v
187  else
188  m_spEigenvector = m_spResidual->clone();
189 
190  // reset Dirichlet rows to 0
191  for(size_t i = 0; i < m_spEigenvector->size(); i++)
192  {
193  if(m_vbDirichlet[i])
194  {
195  (*m_spEigenvector)[i] = 0.0;
196  }
197  }
198 
199  // v = v / ||v||_B
201 
202  // in case B is not set, restore storage type of m_spEigenvector to PST_CONSISTENT
203  // changed by norm calculation in normalize_approximations()
204  if(m_spLinOpB == SPNULL)
205  {
206  #ifdef UG_PARALLEL
207  m_spEigenvector->change_storage_type(PST_CONSISTENT);
208  #endif
209  }
210 
211  // residual = v - v_old
213 
214  if(m_spResidual->norm() <= m_dPrecision)
215  {
216  UG_LOG("PowerMethod::calculate_max_eigenvalue() converged after " << m_iteration << " iterations." << std::endl);
217  break;
218  }
219 
221  UG_LOG("PowerMethod::calculate_max_eigenvalue() reached precision of " << m_spResidual->norm() << " after " << m_maxIterations << " iterations." << std::endl);
222  }
223 
224  // lambda = <v,Av>
227 
228  return true;
229  }
230 
232  {
233  PROFILE_FUNC_GROUP("PowerMethod");
234 
235  UG_COND_THROW(m_spLinOpA == SPNULL, "PowerMethod::calculate_min_eigenvalue(): Linear operator A not set, please specify.");
236  UG_COND_THROW(m_spSolver == SPNULL, "PowerMethod::calculate_min_eigenvalue(): Solver not set, please specify.");
237 
239 
240  m_spSolver->init(m_spLinOpA);
241 
243  {
245 
246  // residual = B v
247  if(m_spLinOpB != SPNULL)
249  // in case B is not set
250  else
251  {
252  m_spResidual = m_spEigenvector->clone();
253 
254  #ifdef UG_PARALLEL
255  m_spResidual->change_storage_type(PST_ADDITIVE);
256  #endif
257  }
258 
259  // v = A^-1 B v or v = A^-1 v
261 
262  // reset Dirichlet rows to 0
263  for(size_t i = 0; i < m_spEigenvector->size(); i++)
264  {
265  if(m_vbDirichlet[i])
266  {
267  (*m_spEigenvector)[i] = 0.0;
268  }
269  }
270 
271  // v = v / ||v||_B
273 
274  // in case B is not set, restore storage type of m_spEigenvector to PST_CONSISTENT
275  // changed by norm calculation in normalize_approximations()
276  if(m_spLinOpB == SPNULL)
277  {
278  #ifdef UG_PARALLEL
279  m_spEigenvector->change_storage_type(PST_CONSISTENT);
280  #endif
281  }
282 
283  // residual = v - v_old
285 
286  if(m_spResidual->norm() <= m_dPrecision)
287  {
288  UG_LOG("PowerMethod::calculate_min_eigenvalue() converged after " << m_iteration << " iterations." << std::endl);
289  break;
290  }
291 
293  UG_LOG("PowerMethod::calculate_min_eigenvalue() reached precision of " << m_spResidual->norm() << " after " << m_maxIterations << " iterations." << std::endl);
294  }
295 
296  // lambda = <v,Av>
299 
300  return true;
301  }
302 
304  {
305  return m_dMaxEigenvalue;
306  }
307 
309  {
310  return m_dMinEigenvalue;
311  }
312 
313  size_t get_iterations()
314  {
315  return m_iteration;
316  }
317 
319  {
320  PrintMatrix(m_spMatOpA->get_matrix(), "A");
321  }
322 
324  {
325  PrintMatrix(m_spMatOpB->get_matrix(), "B");
326  }
327 
329  {
330  m_spEigenvector->print();
331  }
332 
333  private:
334 
336  {
337  SmartPtr<vector_type> t = m_spEigenvector->clone_without_values();
338  t->set(0.0);
339  return t;
340  }
341 
342  double B_norm(vector_type &x)
343  {
344  if(m_spMatOpB != SPNULL)
345  return EnergyNorm(x, *m_spMatOpB);
346  else
347  return x.norm();
348  }
349 
351  {
353  }
354 };
355 
356 } // namespace ug
357 
358 
359 #endif // __H__UG__LIB_ALGEBRA__POWER_METHOD_H__
Definition: smart_pointer.h:108
describes a linear mapping X->Y
Definition: linear_operator.h:80
describes an inverse linear mapping X->Y
Definition: linear_operator_inverse.h:80
Definition: matrix_operator.h:49
Definition: power_method.h:61
size_t m_iteration
Definition: power_method.h:104
std::vector< bool > m_vbDirichlet
Definition: power_method.h:99
void print_matrix_B()
Definition: power_method.h:323
vector_type::value_type vector_value_type
Definition: power_method.h:69
TAlgebra algebra_type
Definition: power_method.h:64
SmartPtr< ILinearOperatorInverse< vector_type > > m_spSolver
Definition: power_method.h:89
void set_linear_operator_B(SmartPtr< ILinearOperator< vector_type > > loB)
Definition: power_method.h:144
void print_matrix_A()
Definition: power_method.h:318
double get_min_eigenvalue()
Definition: power_method.h:308
SmartPtr< matrix_operator_type > m_spMatOpB
Definition: power_method.h:79
SmartPtr< ILinearOperator< vector_type > > m_spLinOpA
Definition: power_method.h:73
int calculate_min_eigenvalue()
Definition: power_method.h:231
double B_norm(vector_type &x)
Definition: power_method.h:342
size_t get_iterations()
Definition: power_method.h:313
SmartPtr< vector_type > m_spEigenvector
Definition: power_method.h:83
TAlgebra::matrix_type matrix_type
Definition: power_method.h:68
int calculate_max_eigenvalue()
Definition: power_method.h:165
TAlgebra::vector_type vector_type
Definition: power_method.h:67
void set_linear_operator_A(SmartPtr< ILinearOperator< vector_type > > loA)
Definition: power_method.h:126
void print_eigenvector()
Definition: power_method.h:328
void set_solver(SmartPtr< ILinearOperatorInverse< vector_type > > solver)
Definition: power_method.h:121
int m_numDirichletRows
Definition: power_method.h:100
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Definition: power_method.h:76
SmartPtr< vector_type > m_spResidual
Definition: power_method.h:96
double m_dMinEigenvalue
Definition: power_method.h:86
size_t m_maxIterations
Definition: power_method.h:92
void set_max_iterations(size_t maxIterations)
Definition: power_method.h:155
SmartPtr< matrix_operator_type > m_spMatOpA
Definition: power_method.h:78
double m_dPrecision
Definition: power_method.h:93
PowerMethod()
Definition: power_method.h:106
SmartPtr< ILinearOperator< vector_type > > m_spLinOpB
Definition: power_method.h:77
double get_max_eigenvalue()
Definition: power_method.h:303
void set_precision(double precision)
Definition: power_method.h:160
void set_start_vector(SmartPtr< vector_type > vec)
Definition: power_method.h:150
SmartPtr< vector_type > create_approximation_vector()
Definition: power_method.h:335
double m_dMaxEigenvalue
Definition: power_method.h:85
void normalize_approximations()
Definition: power_method.h:350
SmartPtr< vector_type > m_spEigenvectorOld
Definition: power_method.h:84
@ PST_CONSISTENT
Definition: parallel_storage_type.h:68
@ PST_ADDITIVE
Definition: parallel_storage_type.h:69
const NullSmartPtr SPNULL
The equivalent to NULL for smart pointers.
Definition: smart_pointer.h:90
#define UG_LOG(msg)
Definition: log.h:367
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition: error.h:61
void VecScaleAdd(vector_t &vOut, typename vector_t::value_type s1, const vector_t &v1, typename vector_t::value_type s2, const vector_t &v2)
Scales two Vectors, adds them and returns the sum in a third vector.
Definition: math_vector_functions_common_impl.hpp:265
the ug namespace
double EnergyNorm(vector_type &x, matrix_type &A, vector_type &tmp)
Definition: additional_math.h:115
void PrintMatrix(const matrix_type &mat, const char *name)
Definition: additional_math.h:201
#define PROFILE_FUNC_GROUP(groups)
Definition: profiler.h:258
T value_type
Definition: sparsematrix_interface.h:2