33#ifndef __H__UG__CPU_ALGEBRA__LU_DECOMP_H__
34#define __H__UG__CPU_ALGEBRA__LU_DECOMP_H__
36#include "../small_matrix/densematrix.h"
37#include "../small_matrix/densevector.h"
38#include "../../common/operations.h"
42#ifndef LAPACK_AVAILABLE
47template<
typename matrix_t>
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;
54 size_t n = A.num_rows();
59 for(
size_t k=0; k<n; k++)
62 for(
size_t j=k+1; j<n; j++)
63 if(
dabs(A(biggest, k)) <
dabs(A(j,k)))
66 for(
size_t j=0; j<n; j++)
67 std::swap(A(k, j), A(biggest, j));
69 pInterchange[k] = biggest;
70 if(
dabs(A(k,k)) < 1e-10)
73 for(
size_t i=k+1; i<n; i++)
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);
83 for(
size_t k=0; k<n; k++)
85 if(
dabs(A(k,k)) < 1e-10)
88 for(
size_t i=k+1; i<n; i++)
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);
99template<
typename matrix_t>
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;
106 size_t n = A.num_rows();
111 for(
size_t i=0; i < n; i++)
113 UG_LOG(
"i=" << i <<
": \n")
115 UG_LOG(
"A(i,i) = " << A(i,i) <<
"\n");
116 for(
size_t j=i+1; j<n; j++)
118 UG_LOG(
"A(" << j <<
", " << i <<
")=" << A(j,i) <<
"\n");
119 if(
dabs(A(biggest, i)) <
dabs(A(j,i))) biggest = j;
121 UG_LOG(biggest <<
" is biggest.");
124 for(
size_t j=0; j<n; j++)
125 std::swap(A(i, j), A(biggest, j));
127 pInterchange[i] = biggest;
128 if(
dabs(A(i,i)) < 1e-10)
132 for(
size_t k=0; k<i; k++)
137 double a_ik = (A(i,k) /= A(k,k));
139 for(
size_t j=k+1; j<n; j++)
140 A(i,j) -= A(k,j) * a_ik;
148 for(
size_t i=0; i < n; i++)
150 if(
dabs(A(i,i)) < 1e-10)
153 for(
size_t k=0; k<i; k++)
158 double a_ik = (A(i,k) /= A(k,k));
160 for(
size_t j=k+1; j<n; j++)
161 A(i,j) -= A(k,j) * a_ik;
168template<
typename matrix_t,
typename vector_t>
171 size_t n=A.num_rows();
174 for(
size_t i=0; i<n; i++)
175 if(i < pInterchange[i])
176 std::swap(x[i], x[pInterchange[i]]);
181 for(
size_t i=0; i<n; i++)
184 for(
size_t k=0; k<i; k++)
191 for(
size_t i=n-1; ; i--)
194 for(
size_t k=i+1; k<n; k++)
209template<
typename TStorage>
210class DenseMatrixInverse
250 UG_ASSERT(mat.num_rows() == mat.num_cols(),
"only for square matrices");
264 UG_LOG(
"ERROR in 'DenseMatrixInverse::invert': Matrix is singular, "
265 "cannot calculate Inverse.\n");
274 template<
typename vector_t>
290std::ostream &
operator << (std::ostream &out,
const DenseMatrixInverse<T> &mat)
293 for(
size_t r=0; r<mat.num_rows(); ++r)
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 <<
"| ";
300 out <<
" (DenseMatrixInverse " << mat.num_rows() <<
"x" << mat.num_cols() <<
", " << ((T::ordering ==
ColMajor) ?
"ColMajor)" :
"RowMajor)");
307struct matrix_algebra_type_traits<DenseMatrixInverse<T> >
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
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