Line data Source code
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 :
40 : #include "lib_algebra/operator/interface/linear_operator.h"
41 : #include "lib_algebra/operator/interface/linear_operator_inverse.h"
42 : #include "lib_algebra/operator/interface/preconditioner.h"
43 : #include "lib_algebra/operator/interface/matrix_operator.h"
44 :
45 : #include "lib_algebra/algebra_common/vector_util.h"
46 : #include "common/profiler/profiler.h"
47 :
48 :
49 : namespace ug{
50 :
51 :
52 : /**
53 : * Power Method Eigensolver
54 : *
55 : * The Power Method solves generalized eigenvalue problems of the form A v = lambda B v
56 : * calculating the largest/smallest (in terms of abs(lambda)) eigenvalues of the problem
57 : * for sparse matrices A and B which e.g. emerge from FE or FV discretizations.
58 : */
59 : template<typename TAlgebra>
60 : class PowerMethod
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;
69 : typedef typename vector_type::value_type vector_value_type;
70 :
71 : private:
72 : // Stiffness matrix
73 : SmartPtr<ILinearOperator<vector_type> > m_spLinOpA;
74 :
75 : // Mass matrix
76 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
77 : SmartPtr<ILinearOperator<vector_type> > m_spLinOpB;
78 : SmartPtr<matrix_operator_type> m_spMatOpA;
79 : SmartPtr<matrix_operator_type> m_spMatOpB;
80 : //matrix_operator_type *m_pB;
81 :
82 : // Eigenvector and -value
83 : SmartPtr<vector_type> m_spEigenvector;
84 : SmartPtr<vector_type> m_spEigenvectorOld;
85 : double m_dMaxEigenvalue;
86 : double m_dMinEigenvalue;
87 :
88 : // Solver
89 : SmartPtr<ILinearOperatorInverse<vector_type> > m_spSolver;
90 :
91 : // Precision
92 : size_t m_maxIterations;
93 : double m_dPrecision;
94 :
95 : // Residual
96 : SmartPtr<vector_type> m_spResidual;
97 :
98 : // Dirichlet node indexing
99 : std::vector<bool> m_vbDirichlet;
100 : int m_numDirichletRows;
101 :
102 : public:
103 : // Number of iterations
104 : size_t m_iteration;
105 :
106 0 : PowerMethod()
107 : {
108 : UG_LOG("Initializing PowerMethod." << std::endl);
109 0 : m_spLinOpA = SPNULL;
110 0 : m_spLinOpB = SPNULL;
111 0 : m_spMatOpA = SPNULL;
112 0 : m_spMatOpB = SPNULL;
113 0 : m_maxIterations = 1000;
114 0 : m_dPrecision = 1e-8;
115 0 : m_iteration = 0;
116 0 : m_dMaxEigenvalue = 0.0;
117 0 : m_dMinEigenvalue = 0.0;
118 0 : m_numDirichletRows = 0;
119 0 : }
120 :
121 0 : void set_solver(SmartPtr<ILinearOperatorInverse<vector_type> > solver)
122 : {
123 0 : m_spSolver = solver;
124 0 : }
125 :
126 0 : void set_linear_operator_A(SmartPtr<ILinearOperator<vector_type> > loA)
127 : {
128 0 : m_spLinOpA = loA;
129 :
130 : // get dirichlet nodes
131 0 : m_spMatOpA = m_spLinOpA.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
132 0 : matrix_type& A = m_spMatOpA->get_matrix();
133 0 : m_vbDirichlet.resize(A.num_rows());
134 :
135 0 : for(size_t i=0; i<A.num_rows(); i++)
136 : {
137 0 : m_vbDirichlet[i] = A.is_isolated(i);
138 :
139 0 : if(A.is_isolated(i) == 0)
140 0 : m_numDirichletRows++;
141 : }
142 0 : }
143 :
144 0 : void set_linear_operator_B(SmartPtr<ILinearOperator<vector_type> > loB)
145 : {
146 0 : m_spLinOpB = loB;
147 0 : m_spMatOpB = m_spLinOpB.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
148 0 : }
149 :
150 0 : void set_start_vector(SmartPtr<vector_type> vec)
151 : {
152 0 : m_spEigenvector = vec;
153 0 : }
154 :
155 0 : void set_max_iterations(size_t maxIterations)
156 : {
157 0 : m_maxIterations = maxIterations;
158 0 : }
159 :
160 0 : void set_precision(double precision)
161 : {
162 0 : m_dPrecision = precision;
163 0 : }
164 :
165 0 : int calculate_max_eigenvalue()
166 : {
167 : PROFILE_FUNC_GROUP("PowerMethod");
168 :
169 0 : UG_COND_THROW(m_spSolver == SPNULL && m_spLinOpB != SPNULL, "PowerMethod::calculate_max_eigenvalue(): Solver not set, please specify.");
170 0 : if(m_spLinOpB != SPNULL)
171 0 : m_spSolver->init(m_spLinOpB);
172 :
173 0 : m_spResidual = create_approximation_vector();
174 :
175 0 : for(m_iteration = 0; m_iteration < m_maxIterations; ++m_iteration)
176 : {
177 0 : m_spEigenvectorOld = m_spEigenvector->clone();
178 :
179 : // residual = A v
180 0 : UG_COND_THROW(m_spLinOpA == SPNULL, "PowerMethod::calculate_max_eigenvalue(): Linear operator A not set, please specify.");
181 0 : m_spLinOpA->apply(*m_spResidual, *m_spEigenvector);
182 :
183 : // v = B^-1 A v
184 0 : if(m_spLinOpB != SPNULL)
185 0 : m_spSolver->apply(*m_spEigenvector, *m_spResidual);
186 : // in case B is not set: v = A v
187 : else
188 0 : m_spEigenvector = m_spResidual->clone();
189 :
190 : // reset Dirichlet rows to 0
191 0 : for(size_t i = 0; i < m_spEigenvector->size(); i++)
192 : {
193 0 : if(m_vbDirichlet[i])
194 : {
195 0 : (*m_spEigenvector)[i] = 0.0;
196 : }
197 : }
198 :
199 : // v = v / ||v||_B
200 0 : normalize_approximations();
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
212 0 : VecScaleAdd(*m_spResidual, 1.0, *m_spEigenvectorOld, -1.0, *m_spEigenvector);
213 :
214 0 : if(m_spResidual->norm() <= m_dPrecision)
215 : {
216 0 : UG_LOG("PowerMethod::calculate_max_eigenvalue() converged after " << m_iteration << " iterations." << std::endl);
217 : break;
218 : }
219 :
220 0 : if(m_iteration == m_maxIterations-1)
221 0 : UG_LOG("PowerMethod::calculate_max_eigenvalue() reached precision of " << m_spResidual->norm() << " after " << m_maxIterations << " iterations." << std::endl);
222 : }
223 :
224 : // lambda = <v,Av>
225 0 : m_spLinOpA->apply(*m_spResidual, *m_spEigenvector);
226 0 : m_dMaxEigenvalue = m_spEigenvector->dotprod(*m_spResidual);
227 :
228 0 : return true;
229 : }
230 :
231 0 : int calculate_min_eigenvalue()
232 : {
233 : PROFILE_FUNC_GROUP("PowerMethod");
234 :
235 0 : UG_COND_THROW(m_spLinOpA == SPNULL, "PowerMethod::calculate_min_eigenvalue(): Linear operator A not set, please specify.");
236 0 : UG_COND_THROW(m_spSolver == SPNULL, "PowerMethod::calculate_min_eigenvalue(): Solver not set, please specify.");
237 :
238 0 : m_spResidual = create_approximation_vector();
239 :
240 0 : m_spSolver->init(m_spLinOpA);
241 :
242 0 : for(m_iteration = 0; m_iteration < m_maxIterations; ++m_iteration)
243 : {
244 0 : m_spEigenvectorOld = m_spEigenvector->clone();
245 :
246 : // residual = B v
247 0 : if(m_spLinOpB != SPNULL)
248 0 : m_spLinOpB->apply(*m_spResidual, *m_spEigenvector);
249 : // in case B is not set
250 : else
251 : {
252 0 : 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
260 0 : m_spSolver->apply(*m_spEigenvector, *m_spResidual);
261 :
262 : // reset Dirichlet rows to 0
263 0 : for(size_t i = 0; i < m_spEigenvector->size(); i++)
264 : {
265 0 : if(m_vbDirichlet[i])
266 : {
267 0 : (*m_spEigenvector)[i] = 0.0;
268 : }
269 : }
270 :
271 : // v = v / ||v||_B
272 0 : normalize_approximations();
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
284 0 : VecScaleAdd(*m_spResidual, 1.0, *m_spEigenvectorOld, -1.0, *m_spEigenvector);
285 :
286 0 : if(m_spResidual->norm() <= m_dPrecision)
287 : {
288 0 : UG_LOG("PowerMethod::calculate_min_eigenvalue() converged after " << m_iteration << " iterations." << std::endl);
289 : break;
290 : }
291 :
292 0 : if(m_iteration == m_maxIterations-1)
293 0 : UG_LOG("PowerMethod::calculate_min_eigenvalue() reached precision of " << m_spResidual->norm() << " after " << m_maxIterations << " iterations." << std::endl);
294 : }
295 :
296 : // lambda = <v,Av>
297 0 : m_spLinOpA->apply(*m_spResidual, *m_spEigenvector);
298 0 : m_dMinEigenvalue = m_spEigenvector->dotprod(*m_spResidual);
299 :
300 0 : return true;
301 : }
302 :
303 0 : double get_max_eigenvalue()
304 : {
305 0 : return m_dMaxEigenvalue;
306 : }
307 :
308 0 : double get_min_eigenvalue()
309 : {
310 0 : return m_dMinEigenvalue;
311 : }
312 :
313 0 : size_t get_iterations()
314 : {
315 0 : return m_iteration;
316 : }
317 :
318 0 : void print_matrix_A()
319 : {
320 0 : PrintMatrix(m_spMatOpA->get_matrix(), "A");
321 0 : }
322 :
323 0 : void print_matrix_B()
324 : {
325 0 : PrintMatrix(m_spMatOpB->get_matrix(), "B");
326 0 : }
327 :
328 0 : void print_eigenvector()
329 : {
330 0 : m_spEigenvector->print();
331 0 : }
332 :
333 : private:
334 :
335 0 : SmartPtr<vector_type> create_approximation_vector()
336 : {
337 0 : SmartPtr<vector_type> t = m_spEigenvector->clone_without_values();
338 : t->set(0.0);
339 0 : return t;
340 : }
341 :
342 0 : double B_norm(vector_type &x)
343 : {
344 0 : if(m_spMatOpB != SPNULL)
345 0 : return EnergyNorm(x, *m_spMatOpB);
346 : else
347 0 : return x.norm();
348 : }
349 :
350 0 : void normalize_approximations()
351 : {
352 0 : *m_spEigenvector *= 1/ B_norm(*m_spEigenvector);
353 0 : }
354 : };
355 :
356 : } // namespace ug
357 :
358 :
359 : #endif // __H__UG__LIB_ALGEBRA__POWER_METHOD_H__
|
