lu.hpp
Go to the documentation of this file.
1 /*!
2  * \brief Implements the lu decomposiion algorithm
3  *
4  * \author O. Krause
5  * \date 2013
6  *
7  *
8  * \par Copyright 1995-2015 Shark Development Team
9  *
10  * <BR><HR>
11  * This file is part of Shark.
12  * <http://image.diku.dk/shark/>
13  *
14  * Shark is free software: you can redistribute it and/or modify
15  * it under the terms of the GNU Lesser General Public License as published
16  * by the Free Software Foundation, either version 3 of the License, or
17  * (at your option) any later version.
18  *
19  * Shark is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU Lesser General Public License for more details.
23  *
24  * You should have received a copy of the GNU Lesser General Public License
25  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
26  *
27  */
28 #ifndef SHARK_LINALG_BLAS_LU_HPP
29 #define SHARK_LINALG_BLAS_LU_HPP
30 
31 #include "permutation.hpp"
32 #include "vector_proxy.hpp"
33 #include "matrix_proxy.hpp"
34 
35 namespace shark {
36 namespace blas {
37 
38 // LU factorization without pivoting
39 template<class M>
40 typename M::size_type lu_factorize(M &m) {
41  typedef typename M::size_type size_type;
42  typedef typename M::value_type value_type;
43 
44  size_type singular = 0;
45  size_type size1 = m.size1();
46  size_type size2 = m.size2();
47  size_type size = (std::min)(size1, size2);
48  for (size_type i = 0; i < size; ++ i) {
49  matrix_column<M> mci(column(m, i));
50  matrix_row<M> mri(row(m, i));
51  if (m(i, i) != value_type/*zero*/()) {
52  value_type m_inv = value_type(1) / m(i, i);
53  subrange(mci, i + 1, size1) *= m_inv;
54  } else if (singular == 0) {
55  singular = i + 1;
56  }
57  noalias(subrange(m, i + 1, size1, i + 1, size2))-=(
58  outer_prod(subrange(mci, i + 1, size1),
59  subrange(mri, i + 1, size2)));
60  }
61  return singular;
62 }
63 
64 // LU factorization with partial pivoting
65 template<class M, class PM>
66 typename M::size_type lu_factorize(M &m, PM &pm) {
67  typedef typename M::size_type size_type;
68  typedef typename M::value_type value_type;
69 
70  size_type singular = 0;
71  size_type size1 = m.size1();
72  size_type size2 = m.size2();
73  size_type size = (std::min)(size1, size2);
74  for (size_type i = 0; i < size; ++ i) {
75  matrix_column<M> mci(column(m, i));
76  matrix_row<M> mri(row(m, i));
77  size_type i_norm_inf = i + index_norm_inf(subrange(mci, i, size1));
78  SIZE_CHECK(i_norm_inf < size1);
79  if (m(i_norm_inf, i) != value_type/*zero*/()) {
80  if (i_norm_inf != i) {
81  pm(i) = i_norm_inf;
82  swap_rows(m,i_norm_inf,i);
83  } else {
84  SIZE_CHECK(pm(i) == i_norm_inf);
85  }
86  value_type m_inv = value_type(1) / m(i, i);
87  subrange(mci, i + 1, size1) *= m_inv;
88  } else if (singular == 0) {
89  singular = i + 1;
90  }
91  noalias(subrange(m,i + 1, size1, i + 1, size2))-=(
92  outer_prod(subrange(mci, i + 1, size1),
93  subrange(mri, i + 1, size2)));
94  }
95  return singular;
96 }
97 }
98 }
99 
100 #endif