Line data Source code
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 :
44 : namespace ug
45 : {
46 :
47 : template<typename matrix_t>
48 0 : bool 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 0 : if(A.num_rows() != A.num_cols()) return false;
53 :
54 : size_t n = A.num_rows();
55 :
56 0 : if(pInterchange)
57 : {
58 0 : pInterchange[0] = 0;
59 0 : for(size_t k=0; k<n; k++)
60 : {
61 : size_t biggest = k;
62 0 : for(size_t j=k+1; j<n; j++)
63 0 : if(dabs(A(biggest, k)) < dabs(A(j,k)))
64 : biggest = j; // costly.
65 0 : if(biggest != k)
66 0 : for(size_t j=0; j<n; j++)
67 : std::swap(A(k, j), A(biggest, j));
68 :
69 0 : pInterchange[k] = biggest;
70 0 : if(dabs(A(k,k)) < 1e-10)
71 : return false;
72 :
73 0 : for(size_t i=k+1; i<n; i++)
74 : {
75 0 : A(i,k) = A(i,k)/A(k,k);
76 0 : for(size_t j=k+1; j<n; j++)
77 0 : A(i,j) = A(i,j) - A(i,k)*A(k,j);
78 : }
79 : }
80 : }
81 : else
82 : {
83 0 : for(size_t k=0; k<n; k++)
84 : {
85 0 : if(dabs(A(k,k)) < 1e-10)
86 : return false;
87 :
88 0 : for(size_t i=k+1; i<n; i++)
89 : {
90 0 : A(i,k) = A(i,k)/A(k,k);
91 0 : for(size_t j=k+1; j<n; j++)
92 0 : A(i,j) = A(i,j) - A(i,k)*A(k,j);
93 : }
94 : }
95 : }
96 : return true;
97 : }
98 :
99 : template<typename matrix_t>
100 : bool 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 : }
168 : template<typename matrix_t, typename vector_t>
169 0 : bool SolveLU(const DenseMatrix<matrix_t> &A, DenseVector<vector_t> &x, const size_t *pInterchange)
170 : {
171 : size_t n=A.num_rows();
172 :
173 0 : if(pInterchange)
174 0 : for(size_t i=0; i<n; i++)
175 0 : 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 0 : for(size_t i=0; i<n; i++)
182 : {
183 0 : s = x[i];
184 0 : for(size_t k=0; k<i; k++)
185 0 : s -= A(i, k)*x[k];
186 0 : x[i] = s;
187 : }
188 :
189 : // -> x = U^{-1} (L^{-1} b)
190 : // solve upper right
191 0 : for(size_t i=n-1; ; i--)
192 : {
193 0 : s = x[i];
194 0 : for(size_t k=i+1; k<n; k++)
195 0 : s -= A(i, k)*x[k];
196 0 : x[i] = s/A(i,i);
197 0 : if(i==0) break;
198 : }
199 :
200 0 : return true;
201 : }
202 :
203 :
204 : ///////////////////////////////////////////////////////////////////////////////////////
205 : /**
206 : * smallInverse<size_t n>
207 : * A class to hold a inverse of a smallMatrix<n>
208 : */
209 : template<typename TStorage>
210 : class DenseMatrixInverse
211 : {
212 : private: // storage
213 : DenseMatrix<TStorage> densemat;
214 : std::vector<size_t> interchange;
215 :
216 : public:
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 0 : inline void resize(size_t k, size_t k2)
228 : {
229 0 : UG_COND_THROW(k!=k2, "only square matrices supported");
230 0 : resize(k);
231 0 : }
232 0 : inline void resize(size_t k)
233 : {
234 0 : densemat.resize(k,k);
235 0 : densemat = 0.0;
236 0 : }
237 :
238 : double &operator()(int r, int c)
239 : {
240 0 : return densemat(r,c);
241 : }
242 : const double &operator()(int r, int c) const
243 : {
244 : return densemat(r,c);
245 : }
246 : public:
247 : //! initializes this object as inverse of mat
248 : bool set_as_inverse_of(const DenseMatrix<TStorage> &mat)
249 : {
250 : UG_ASSERT(mat.num_rows() == mat.num_cols(), "only for square matrices");
251 0 : densemat = mat;
252 0 : return invert();
253 : }
254 :
255 0 : bool invert()
256 : {
257 0 : interchange.resize(densemat.num_rows());
258 :
259 0 : if(!interchange.empty()){
260 0 : bool bLUDecomp = LUDecomp(densemat, &interchange[0]);
261 :
262 0 : if(bLUDecomp!=true)
263 : {
264 : UG_LOG("ERROR in 'DenseMatrixInverse::invert': Matrix is singular, "
265 : "cannot calculate Inverse.\n");
266 : }
267 :
268 0 : return bLUDecomp;
269 : }
270 : else
271 : return true;
272 : }
273 :
274 : template<typename vector_t>
275 : void apply(DenseVector<vector_t> &vec) const
276 : {
277 0 : if(!interchange.empty())
278 0 : 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 :
289 : template<typename T>
290 : std::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 :
306 : template<typename T>
307 : struct 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__
|
