Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Rupp
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__LIB_ALGEBRA__LAPACK_LU_OPERATOR__
34 : #define __H__LIB_ALGEBRA__LAPACK_LU_OPERATOR__
35 : #include <iostream>
36 : #include <sstream>
37 :
38 : #include "common/common.h"
39 : #include "lib_algebra/operator/interface/matrix_operator_inverse.h"
40 :
41 : #ifdef UG_PARALLEL
42 : #include "lib_algebra/parallelization/parallelization.h"
43 : #endif
44 : #include "../preconditioner/ilut_scalar.h"
45 : #include "../interface/preconditioned_linear_operator_inverse.h"
46 : #include "linear_solver.h"
47 :
48 : #include "lib_algebra/cpu_algebra_types.h"
49 :
50 : #include "../operator_util.h"
51 : namespace ug{
52 :
53 : template <typename TAlgebra>
54 : class LU
55 : : public IMatrixOperatorInverse<typename TAlgebra::matrix_type,
56 : typename TAlgebra::vector_type>
57 : {
58 : public:
59 : /// Algebra type
60 : typedef TAlgebra algebra_type;
61 :
62 : /// Vector type
63 : typedef typename TAlgebra::vector_type vector_type;
64 :
65 : /// Matrix type
66 : typedef typename TAlgebra::matrix_type matrix_type;
67 :
68 : /// Base type
69 : typedef IMatrixOperatorInverse<matrix_type,vector_type> base_type;
70 :
71 : using base_type::init;
72 :
73 : protected:
74 : using base_type::convergence_check;
75 :
76 : public:
77 : /// constructor
78 0 : LU() : m_spOperator(NULL), m_mat(), m_bSortSparse(true), m_bInfo(false), m_bShowProgress(true)
79 : {
80 : #ifdef LAPACK_AVAILABLE
81 : m_iMinimumForSparse = 4000;
82 : #else
83 0 : m_iMinimumForSparse = 1000;
84 : #endif
85 0 : };
86 :
87 : /// returns if parallel solving is supported
88 0 : virtual bool supports_parallel() const {return false;}
89 :
90 : ///
91 0 : void set_minimum_for_sparse(size_t N)
92 : {
93 0 : m_iMinimumForSparse=N;
94 0 : }
95 :
96 0 : void set_sort_sparse(bool b)
97 : {
98 0 : m_bSortSparse = b;
99 0 : }
100 :
101 0 : void set_info(bool b)
102 : {
103 0 : m_bInfo = b;
104 0 : }
105 :
106 0 : void set_show_progress(bool b)
107 : {
108 0 : m_bShowProgress = b;
109 0 : }
110 :
111 0 : virtual const char* name() const {return "LU";}
112 :
113 : private:
114 :
115 0 : void print_info(const matrix_type &A)
116 : {
117 : size_t blockSize = block_traits<typename matrix_type::value_type>::static_num_rows;
118 : UG_LOG("Matrix " << A.num_rows() << " x " << A.num_rows() << " with " << blockSize << " x " << blockSize << " blocks");
119 0 : }
120 :
121 : /// returns name of solver
122 0 : bool init_dense(const matrix_type &A)
123 : {
124 : PROFILE_FUNC();
125 0 : m_bDense = true;
126 :
127 0 : if(m_bInfo)
128 : {
129 : UG_LOG("LU, using DenseLU on ");
130 0 : print_info(A);
131 0 : UG_LOG("\n DenseLU needs " << GetBytesSizeString(m_size*m_size*sizeof(double)) << " of memory.\n");
132 : }
133 :
134 0 : m_size = GetDenseDoubleFromSparse(m_mat, A);
135 :
136 0 : if(m_mat.invert() == false)
137 : {
138 0 : UG_THROW("ERROR in Matrix is singular");
139 : return false;
140 : }
141 : else
142 : return true;
143 : }
144 :
145 : void init_var(const matrix_type &A)
146 : {
147 : UG_ASSERT(0, "not yet tested");
148 : /*m_size = 0;
149 : std::vector<size_t> blockbegin(A.num_rows()+1);
150 :
151 : for(size_t i=0; i<A.num_rows(); i++)
152 : {
153 : bool bFound;
154 : typename matrix_type::const_row_iterator it = A.get_connection(i,i, bFound);
155 : UG_ASSERT(bFound, "Matrix has to have entry A(" << i << ", " << i << ")");
156 : size_t s = GetRows(it.value());
157 : UG_ASSERT(s == GetCols(it.value()), "diagonal elements have to be square");
158 : if(i == 0)
159 : blockbegin[i] = m_size;
160 : m_size += s;
161 : }
162 : blockbegin[A.num_rows()] = m_size;
163 :
164 : m_mat.resize(m_size);
165 :
166 : for(size_t r=0; r<A.num_rows(); r++)
167 : for(typename matrix_type::const_row_iterator it = A.begin_row(r); it != A.end_row(r); ++it)
168 : {
169 : size_t c = it.index();
170 : const typename matrix_type::value_type &val = it.value();
171 : UG_ASSERT(blockbegin[r]+GetRows(val) == blockbegin[r+1], "blocksizes in matrix inconsistent");
172 : UG_ASSERT(blockbegin[c]+GetCols(val) == blockbegin[c+1], "blocksizes in matrix inconsistent");
173 : for(size_t r2=0; r2 < GetRows(val); r2++)
174 : for(size_t c2=0; c2 < GetCols(val); c2++)
175 : m_mat(blockbegin[r] + r2, blockbegin[c]+c2) = BlockRef(val, r2, c2);
176 : }
177 : */
178 : }
179 :
180 :
181 0 : bool init_sparse(const matrix_type &A)
182 : {
183 : try{
184 : PROFILE_FUNC();
185 0 : m_bDense = false;
186 :
187 0 : if(m_bInfo)
188 : {
189 : UG_LOG("LU using Sparse LU on ");
190 0 : print_info(A);
191 : UG_LOG("\n");
192 : }
193 0 : ilut_scalar = make_sp(new ILUTScalarPreconditioner<algebra_type>(0.0));
194 0 : ilut_scalar->set_sort(m_bSortSparse);
195 0 : ilut_scalar->set_info(m_bInfo);
196 0 : ilut_scalar->set_show_progress(m_bShowProgress);
197 0 : ilut_scalar->preprocess(A);
198 :
199 0 : }UG_CATCH_THROW("LU::" << __FUNCTION__ << " failed")
200 0 : return true;
201 : }
202 :
203 0 : bool solve_dense(vector_type &x, const vector_type &b)
204 : {
205 : try{
206 : PROFILE_FUNC();
207 0 : x = b;
208 0 : m_tmp.resize(m_size);
209 0 : for(size_t i=0, k=0; i<b.size(); i++)
210 : {
211 0 : for(size_t j=0; j<GetSize(b[i]); j++)
212 0 : m_tmp[k++] = BlockRef(b[i],j);
213 : }
214 0 : m_mat.apply(m_tmp);
215 :
216 :
217 0 : for(size_t i=0, k=0; i<b.size(); i++)
218 : {
219 0 : for(size_t j=0; j<GetSize(b[i]); j++)
220 0 : BlockRef(x[i],j) = m_tmp[k++];
221 : }
222 :
223 0 : }UG_CATCH_THROW("LU::" << __FUNCTION__ << " failed")
224 0 : return true;
225 : }
226 :
227 : bool solve_sparse(vector_type &x, const vector_type &b)
228 : {
229 : PROFILE_FUNC();
230 0 : ilut_scalar->solve(x, b);
231 : return true;
232 : }
233 :
234 : public:
235 : /// initializes the solver for a matrix A
236 0 : bool init_lu(const matrix_type *pA)
237 : {
238 : try{
239 : // get matrix of Operator
240 0 : m_pMatrix = pA;
241 0 : if(m_pMatrix->num_rows() == 0) return true;
242 :
243 : // check that matrix exist
244 : if(m_pMatrix == NULL)
245 : {
246 : UG_LOG("ERROR in 'LU::init': No Matrix given.\n");
247 : return false;
248 : }
249 :
250 : const matrix_type &A = *pA;
251 :
252 : PROFILE_FUNC();
253 : PROFILE_BEGIN_GROUP(LU_init_lu, "algebra lu");
254 : if(block_traits<typename vector_type::value_type>::is_static)
255 : {
256 : const size_t nrOfRows = block_traits<typename matrix_type::value_type>::static_num_rows;
257 : UG_ASSERT(nrOfRows == block_traits<typename matrix_type::value_type>::static_num_cols, "only square matrices supported");
258 0 : m_size = A.num_rows() * nrOfRows;
259 :
260 0 : if(m_size > m_iMinimumForSparse)
261 0 : init_sparse(A);
262 : else
263 0 : init_dense(A);
264 : }
265 : else
266 : init_var(A);
267 0 : }UG_CATCH_THROW("LU::" << __FUNCTION__ << " failed")
268 : return true;
269 : }
270 :
271 0 : bool apply_lu(vector_type &x, const vector_type &b)
272 : {
273 0 : if(m_pMatrix->num_rows() == 0) return true;
274 : try{
275 : PROFILE_BEGIN_GROUP(LU_apply_lu, "algebra lu");
276 : #ifndef NDEBUG
277 : if(block_traits<typename vector_type::value_type>::is_static)
278 : {
279 : const size_t static_size = block_traits<typename vector_type::value_type>::static_size;
280 : UG_ASSERT(m_size == b.size() * static_size && m_size == x.size() * static_size,
281 : " wrong size! has to be " << m_size << ", but is " << b << " and " << x);
282 : }
283 : else
284 : {
285 : size_t b_size = 0;
286 : for(size_t i=0; i<b.size(); i++)
287 : {
288 : UG_ASSERT(GetSize(b[i]) == GetSize(x[i]), "wrong size! Sizes of b and x must be the same, but is "
289 : << GetSize(b[i]) << " and " << GetSize(x[i]) << "!");
290 : b_size += GetSize(b[i]);
291 : }
292 : UG_ASSERT(m_size == b_size, " wrong size! has to be " << m_size << ", but is " << b_size << "!");
293 : }
294 : #endif
295 :
296 0 : if(m_bDense)
297 0 : return solve_dense(x, b);
298 : else
299 : return solve_sparse(x, b);
300 :
301 0 : }UG_CATCH_THROW("LU::" << __FUNCTION__ << " failed")
302 : }
303 :
304 : /// set operator L, that will be inverted
305 0 : virtual bool init(SmartPtr<MatrixOperator<matrix_type, vector_type> > Op)
306 : {
307 : // remember operator
308 0 : m_spOperator = Op;
309 :
310 : // init LU operator
311 0 : if(!init_lu(&m_spOperator->get_matrix()))
312 : {
313 : UG_LOG("ERROR in 'LU::init': Cannot init LU Decomposition.\n");
314 0 : return false;
315 : }
316 :
317 : // we're done
318 : return true;
319 : }
320 :
321 : /// Compute u = L^{-1} * f
322 0 : virtual bool apply(vector_type& u, const vector_type& f)
323 : {
324 : PROFILE_FUNC();
325 0 : convergence_check()->set_symbol('%');
326 0 : convergence_check()->set_name("LU Solver");
327 :
328 : #ifdef UG_PARALLEL
329 : if(!f.has_storage_type(PST_ADDITIVE))
330 : {
331 : UG_LOG("ERROR: In 'LU::apply': "
332 : "Inadequate storage format of Vector f.\n");
333 : return false;
334 : }
335 : if(!u.has_storage_type(PST_CONSISTENT))
336 : {
337 : UG_LOG("ERROR: In 'LU::apply': "
338 : "Inadequate storage format of Vector u.\n");
339 : return false;
340 : }
341 : #endif
342 : UG_ASSERT(f.size() == m_pMatrix->num_rows(), "Vector ["<<f.size()<<"] and Row "<<m_pMatrix->num_rows()<<" size mismatch");
343 : UG_ASSERT(u.size() == m_pMatrix->num_cols(), "Vector ["<<u.size()<<"] and Col "<<m_pMatrix->num_cols()<<" size mismatch");
344 : UG_ASSERT(f.size() == u.size(), "Vector sizes have to match!");
345 :
346 0 : if(!apply_lu(u, f))
347 : {
348 : UG_LOG("ERROR in 'LU::apply': "
349 : "Cannot apply LU decomposition.\n");
350 0 : return false;
351 : }
352 :
353 : #ifdef UG_PARALLEL
354 : // todo: we set solution to consistent here, but that is only true for
355 : // serial case. Handle parallel case.
356 : u.set_storage_type(PST_CONSISTENT);
357 : #endif
358 :
359 : // we're done
360 : return true;
361 : }
362 :
363 : /// Compute u = L^{-1} * f AND return defect f := f - L*u
364 0 : virtual bool apply_return_defect(vector_type& u, vector_type& f)
365 : {
366 : PROFILE_BEGIN_GROUP(LU_apply_return_defect, "algebra lu");
367 : // solve u
368 0 : if(!apply(u, f)) return false;
369 :
370 : // update defect
371 0 : if(!m_pMatrix->matmul_minus(f, u))
372 : {
373 : UG_LOG("ERROR in 'LU::apply_return_defect': "
374 : "Cannot apply matmul_minus.\n");
375 : return false;
376 : }
377 :
378 : // we're done
379 0 : return true;
380 : }
381 :
382 0 : virtual std::string config_string() const
383 : {
384 0 : std::stringstream ss;
385 0 : ss << "LU Decomposition: Direct Solver for Linear Equation Systems.\n";
386 0 : ss << " Minimum Entries for Sparse LU: " << m_iMinimumForSparse;
387 0 : if(m_iMinimumForSparse==0)
388 0 : ss << " (= always Sparse LU)";
389 0 : return ss.str();
390 0 : }
391 :
392 :
393 : /// Destructor
394 0 : virtual ~LU() {};
395 :
396 : protected:
397 : /// Operator to invert
398 : SmartPtr<MatrixOperator<matrix_type, vector_type> > m_spOperator;
399 :
400 : /// matrix to invert
401 : const matrix_type* m_pMatrix;
402 :
403 : /// inverse
404 : DenseMatrixInverse<DenseMatrix<VariableArray2<double> > > m_mat;
405 : DenseVector<VariableArray1<double> > m_tmp;
406 : CPUAlgebra::vector_type m_u;
407 : CPUAlgebra::vector_type m_b;
408 : size_t m_size;
409 :
410 : bool m_bDense;
411 : SmartPtr<ILUTScalarPreconditioner<algebra_type> > ilut_scalar;
412 : size_t m_iMinimumForSparse;
413 : bool m_bSortSparse, m_bInfo, m_bShowProgress;
414 : };
415 :
416 : } // end namespace ug
417 :
418 : #endif /* __H__LIB_ALGEBRA__LAPACK_LU_OPERATOR__ */
|
