potrf.hpp
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Implements the default implementation of the POTRF algorithm
5  *
6  * \author O. Krause
7  * \date 2014
8  *
9  *
10  * \par Copyright 1995-2014 Shark Development Team
11  *
12  * <BR><HR>
13  * This file is part of Shark.
14  * <http://image.diku.dk/shark/>
15  *
16  * Shark is free software: you can redistribute it and/or modify
17  * it under the terms of the GNU Lesser General Public License as published
18  * by the Free Software Foundation, either version 3 of the License, or
19  * (at your option) any later version.
20  *
21  * Shark is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU Lesser General Public License for more details.
25  *
26  * You should have received a copy of the GNU Lesser General Public License
27  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
28  *
29  */
30 #ifndef SHARK_LINALG_BLAS_KERNELS_DEFAULT_POTRF_HPP
31 #define SHARK_LINALG_BLAS_KERNELS_DEFAULT_POTRF_HPP
32 
33 #include "../../matrix_proxy.hpp"
34 #include "../../vector_expression.hpp"
35 #include <boost/mpl/bool.hpp>
36 
37 namespace shark {
38 namespace blas {
39 namespace bindings {
40 
41 
42 //upper potrf(row-major)
43 template<class MatA>
44 std::size_t potrf_impl(
45  matrix_expression<MatA>& A,
46  row_major, lower
47 ) {
48  std::size_t m = A().size1();
49  for(size_t j = 0; j < m; j++) {
50  for(size_t i = j; i < m; i++) {
51  double s = A()(i, j);
52  for(size_t k = 0; k < j; k++) {
53  s -= A()(i, k) * A()(j, k);
54  }
55  if(i == j) {
56  if(s <= 0)
57  return i;
58  A()(i, j) = std::sqrt(s);
59  } else {
60  A()(i, j) = s / A()(j , j);
61  }
62  }
63  }
64  return 0;
65 }
66 
67 //lower potrf(row-major)
68 template<class MatA>
69 std::size_t potrf_impl(
70  matrix_expression<MatA>& A,
71  row_major, upper
72 ) {
73  std::size_t m = A().size1();
74  for(size_t i = 0; i < m; i++) {
75  double& Aii = A()(i, i);
76  if(Aii < 0)
77  return i;
78  using std::sqrt;
79  Aii = sqrt(Aii);
80  //update row
81  for(std::size_t j = i + 1; j < m; ++j) {
82  A()(i, j) /= Aii;
83  }
84  //rank-one update
85  for(size_t k = i + 1; k < m; k++) {
86  for(std::size_t j = k; j < m; ++j) {
87  A()(k, j) -= A()(i, k) * A()(i, j);
88  }
89  }
90  }
91  return 0;
92 }
93 
94 
95 //dispatcher for column major
96 template<class MatA, class Triangular>
97 std::size_t potrf_impl(
98  matrix_container<MatA>& A,
99  column_major, Triangular
100 ) {
101  blas::matrix_transpose<MatA> transA(A());
102  return potrf_impl(transA, row_major(), typename Triangular::transposed_orientation());
103 }
104 
105 //dispatcher
106 
107 template <class Triangular, typename MatA>
108 std::size_t potrf(
109  matrix_container<MatA>& A,
110  boost::mpl::false_//unoptimized
111 ) {
112  return potrf_impl(A, typename MatA::orientation(), Triangular());
113 }
114 
115 }
116 }
117 }
118 #endif