Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
lu_decomp.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2010-2014: 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__UG__CPU_ALGEBRA__LU_DECOMP_H__
34#define __H__UG__CPU_ALGEBRA__LU_DECOMP_H__
35
36#include "../small_matrix/densematrix.h"
37#include "../small_matrix/densevector.h"
38#include "../../common/operations.h"
39
40#include <vector>
41
42#ifndef LAPACK_AVAILABLE
43
44namespace ug
45{
46
47template<typename matrix_t>
48bool LUDecomp(DenseMatrix<matrix_t> &A, size_t *pInterchange)
49{
50 // LU Decomposition, IKJ Variant
51 UG_ASSERT(A.num_rows() == A.num_cols(), "LU decomposition only for square matrices");
52 if(A.num_rows() != A.num_cols()) return false;
53
54 size_t n = A.num_rows();
55
56 if(pInterchange)
57 {
58 pInterchange[0] = 0;
59 for(size_t k=0; k<n; k++)
60 {
61 size_t biggest = k;
62 for(size_t j=k+1; j<n; j++)
63 if(dabs(A(biggest, k)) < dabs(A(j,k)))
64 biggest = j; // costly.
65 if(biggest != k)
66 for(size_t j=0; j<n; j++)
67 std::swap(A(k, j), A(biggest, j));
68
69 pInterchange[k] = biggest;
70 if(dabs(A(k,k)) < 1e-10)
71 return false;
72
73 for(size_t i=k+1; i<n; i++)
74 {
75 A(i,k) = A(i,k)/A(k,k);
76 for(size_t j=k+1; j<n; j++)
77 A(i,j) = A(i,j) - A(i,k)*A(k,j);
78 }
79 }
80 }
81 else
82 {
83 for(size_t k=0; k<n; k++)
84 {
85 if(dabs(A(k,k)) < 1e-10)
86 return false;
87
88 for(size_t i=k+1; i<n; i++)
89 {
90 A(i,k) = A(i,k)/A(k,k);
91 for(size_t j=k+1; j<n; j++)
92 A(i,j) = A(i,j) - A(i,k)*A(k,j);
93 }
94 }
95 }
96 return true;
97}
98
99template<typename matrix_t>
100bool LUDecompIKJ(DenseMatrix<matrix_t> &A, size_t *pInterchange)
101{
102 // LU Decomposition, IKJ Variant
103 UG_ASSERT(A.num_rows() == A.num_cols(), "LU decomposition only for square matrices");
104 if(A.num_rows() != A.num_cols()) return false;
105
106 size_t n = A.num_rows();
107
108 if(pInterchange)
109 {
110 pInterchange[0] = 0;
111 for(size_t i=0; i < n; i++)
112 {
113 UG_LOG("i=" << i << ": \n")
114 size_t biggest = i;
115 UG_LOG("A(i,i) = " << A(i,i) << "\n");
116 for(size_t j=i+1; j<n; j++)
117 {
118 UG_LOG("A(" << j << ", " << i << ")=" << A(j,i) << "\n");
119 if(dabs(A(biggest, i)) < dabs(A(j,i))) biggest = j; // costly.
120 }
121 UG_LOG(biggest << " is biggest.");
122
123 if(biggest != i)
124 for(size_t j=0; j<n; j++)
125 std::swap(A(i, j), A(biggest, j));
126
127 pInterchange[i] = biggest;
128 if(dabs(A(i,i)) < 1e-10)
129 return false;
130
131 // eliminate all entries A(i, k) with k<i with rows A(k, .) and k<i
132 for(size_t k=0; k<i; k++)
133 {
134 // add row k to row i by A(i, .) -= A(k,.) A(i,k) / A(k,k)
135 // so that A(i,k) is zero (elements A(i, <k) are already "zero")
136 // safe A(i,k)/A(k,k) in A(i,k)
137 double a_ik = (A(i,k) /= A(k,k));
138
139 for(size_t j=k+1; j<n; j++)
140 A(i,j) -= A(k,j) * a_ik;
141 }
142 }
143 // P A = L R
144 }
145 else
146 {
147 // for all rows
148 for(size_t i=0; i < n; i++)
149 {
150 if(dabs(A(i,i)) < 1e-10)
151 return false;
152 // eliminate all entries A(i, k) with k<i with rows A(k, .) and k<i
153 for(size_t k=0; k<i; k++)
154 {
155 // add row k to row i by A(i, .) -= A(k,.) A(i,k) / A(k,k)
156 // so that A(i,k) is zero (elements A(i, <k) are already "zero")
157 // safe A(i,k)/A(k,k) in A(i,k)
158 double a_ik = (A(i,k) /= A(k,k));
159
160 for(size_t j=k+1; j<n; j++)
161 A(i,j) -= A(k,j) * a_ik;
162 }
163 }
164 }
165
166 return true;
167}
168template<typename matrix_t, typename vector_t>
169bool SolveLU(const DenseMatrix<matrix_t> &A, DenseVector<vector_t> &x, const size_t *pInterchange)
170{
171 size_t n=A.num_rows();
172
173 if(pInterchange)
174 for(size_t i=0; i<n; i++)
175 if(i < pInterchange[i])
176 std::swap(x[i], x[pInterchange[i]]);
177
178 // LU x = b, -> U x = L^{-1} b
179 // solve lower left
180 double s;
181 for(size_t i=0; i<n; i++)
182 {
183 s = x[i];
184 for(size_t k=0; k<i; k++)
185 s -= A(i, k)*x[k];
186 x[i] = s;
187 }
188
189 // -> x = U^{-1} (L^{-1} b)
190 // solve upper right
191 for(size_t i=n-1; ; i--)
192 {
193 s = x[i];
194 for(size_t k=i+1; k<n; k++)
195 s -= A(i, k)*x[k];
196 x[i] = s/A(i,i);
197 if(i==0) break;
198 }
199
200 return true;
201}
202
203
205
209template<typename TStorage>
210class DenseMatrixInverse
211{
212private: // storage
213 DenseMatrix<TStorage> densemat;
214 std::vector<size_t> interchange;
215
216public:
217 inline size_t num_cols() const
218 {
219 return densemat.num_cols();
220 }
221
222 inline size_t num_rows() const
223 {
224 return densemat.num_rows();
225 }
226
227 inline void resize(size_t k, size_t k2)
228 {
229 UG_COND_THROW(k!=k2, "only square matrices supported");
230 resize(k);
231 }
232 inline void resize(size_t k)
233 {
234 densemat.resize(k,k);
235 densemat = 0.0;
236 }
237
238 double &operator()(int r, int c)
239 {
240 return densemat(r,c);
241 }
242 const double &operator()(int r, int c) const
243 {
244 return densemat(r,c);
245 }
246public:
249 {
250 UG_ASSERT(mat.num_rows() == mat.num_cols(), "only for square matrices");
251 densemat = mat;
252 return invert();
253 }
254
255 bool invert()
256 {
257 interchange.resize(densemat.num_rows());
258
259 if(!interchange.empty()){
260 bool bLUDecomp = LUDecomp(densemat, &interchange[0]);
261
262 if(bLUDecomp!=true)
263 {
264 UG_LOG("ERROR in 'DenseMatrixInverse::invert': Matrix is singular, "
265 "cannot calculate Inverse.\n");
266 }
267
268 return bLUDecomp;
269 }
270 else
271 return true;
272 }
273
274 template<typename vector_t>
276 {
277 if(!interchange.empty())
278 SolveLU(densemat, vec, &interchange[0]);
279 }
280
281 // todo: implem
282
283 // todo: implement operator *=
284
285 template<typename T> friend std::ostream &operator << (std::ostream &out, const DenseMatrixInverse<T> &mat);
286};
287
288
289template<typename T>
290std::ostream &operator << (std::ostream &out, const DenseMatrixInverse<T> &mat)
291{
292 out << "[ ";
293 for(size_t r=0; r<mat.num_rows(); ++r)
294 {
295 for(size_t c=0; c<mat.num_cols(); ++c)
296 out << mat.densemat(r, c) << " ";
297 if(r != mat.num_rows()-1) out << "| ";
298 }
299 out << "]";
300 out << " (DenseMatrixInverse " << mat.num_rows() << "x" << mat.num_cols() << ", " << ((T::ordering == ColMajor) ? "ColMajor)" : "RowMajor)");
301
302 return out;
303}
304
305
306template<typename T>
307struct matrix_algebra_type_traits<DenseMatrixInverse<T> >
308{
309 static const int type = MATRIX_USE_GLOBAL_FUNCTIONS;
310};
311
312} // namespace ug
313
314#endif // not LAPACK_AVAILABLE
315
316#endif // __H__UG__CPU_ALGEBRA__LU_DECOMP_H__
parameterString s
Definition densematrix.h:57
Definition lapack_densematrix_inverse.h:51
void resize(size_t k, size_t k2)
Definition lu_decomp.h:227
std::vector< size_t > interchange
Definition lu_decomp.h:214
void resize(size_t k)
Definition lapack_densematrix_inverse.h:67
bool invert()
Definition lapack_densematrix_inverse.h:97
std::vector< lapack_int > interchange
Definition lapack_densematrix_inverse.h:54
void apply(DenseVector< vector_t > &vec) const
Definition lu_decomp.h:275
size_t num_cols() const
Definition lu_decomp.h:217
friend std::ostream & operator<<(std::ostream &out, const DenseMatrixInverse< T > &mat)
Definition lapack_densematrix_inverse.h:130
const double & operator()(int r, int c) const
Definition lu_decomp.h:242
DenseMatrix< TStorage > densemat
Definition lapack_densematrix_inverse.h:53
size_t num_rows() const
Definition lu_decomp.h:222
double & operator()(int r, int c)
Definition lu_decomp.h:238
bool set_as_inverse_of(const DenseMatrix< TStorage > &mat)
initializes this object as inverse of mat
Definition lu_decomp.h:248
Definition densevector.h:101
std::ostream & operator<<(std::ostream &outStream, const ug::MathMatrix< 2, 2 > &m)
Definition math_matrix.cpp:38
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#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
the ug namespace
double dabs(double a)
Definition blocks.h:41
@ ColMajor
Definition storage.h:49
@ MATRIX_USE_GLOBAL_FUNCTIONS
Definition matrix_algebra_types.h:69
bool LUDecomp(DenseMatrix< matrix_t > &A, size_t *pInterchange)
Definition lu_decomp.h:48
bool LUDecompIKJ(DenseMatrix< matrix_t > &A, size_t *pInterchange)
Definition lu_decomp.h:100
bool SolveLU(const DenseMatrix< matrix_t > &A, DenseVector< vector_t > &x, const size_t *pInterchange)
Definition lu_decomp.h:169
static const int type
Definition matrix_algebra_types.h:80