trsv.hpp
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief -
5  *
6  * \author O. Krause
7  * \date 2011
8  *
9  *
10  * \par Copyright 1995-2015 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_TRSV_HPP
31 #define SHARK_LINALG_BLAS_KERNELS_DEFAULT_TRSV_HPP
32 
33 #include "../../matrix_proxy.hpp"
34 #include "../../vector_expression.hpp"
35 #include <boost/mpl/bool.hpp>
36 
37 namespace shark {namespace blas {namespace bindings {
38 
39 //tag encoding
40 // Upper matrix -> boost::mpl::true_
41 // Lower matrix -> boost::mpl::false_
42 
43 //Lower triangular(row-major) - vector
44 template<bool Unit, class TriangularA, class V>
45 void trsv_impl(
46  matrix_expression<TriangularA> const& A,
47  vector_expression<V> &b,
48  boost::mpl::false_, column_major
49 ) {
50  SIZE_CHECK(A().size1() == A().size2());
51  SIZE_CHECK(A().size2() == b().size());
52 
53  typedef typename TriangularA::value_type value_type;
54 
55  std::size_t size = b().size();
56  for (std::size_t n = 0; n != size; ++ n) {
57  if(!Unit){
58  RANGE_CHECK(A()(n, n) != value_type());//matrix is singular
59  b()(n) /= A()(n, n);
60  }
61  if (b()(n) != value_type/*zero*/()){
62  matrix_column<TriangularA const> col = column(A(),n);
63  noalias(subrange(b(),n+1,size)) -= b()(n) * subrange(col,n+1,size);
64  }
65  }
66 }
67 //Lower triangular(column-major) - vector
68 template<bool Unit, class TriangularA, class V>
69 void trsv_impl(
70  matrix_expression<TriangularA> const& A,
71  vector_expression<V> &b,
72  boost::mpl::false_, row_major
73 ) {
74  SIZE_CHECK(A().size1() == A().size2());
75  SIZE_CHECK(A().size2() == b().size());
76 
77  typedef typename TriangularA::value_type value_type;
78 
79  std::size_t size = b().size();
80  for (std::size_t n = 0; n < size; ++ n) {
81  matrix_row<TriangularA const> matRow = row(A(),n);
82  b()(n) -= inner_prod(subrange(matRow,0,n),subrange(b(),0,n));
83 
84  if(!Unit){
85  RANGE_CHECK(A()(n, n) != value_type());//matrix is singular
86  b()(n) /= A()(n, n);
87  }
88  }
89 }
90 
91 //upper triangular(column-major)-vector
92 template<bool Unit, class TriangularA, class V>
93 void trsv_impl(
94  matrix_expression<TriangularA> const& A,
95  vector_expression<V> &b,
96  boost::mpl::true_, column_major
97 ) {
98  SIZE_CHECK(A().size1() == A().size2());
99  SIZE_CHECK(A().size2() == b().size());
100 
101  typedef typename TriangularA::value_type value_type;
102 
103  std::size_t size = b().size();
104  for (std::size_t i = 0; i < size; ++ i) {
105  std::size_t n = size-i-1;
106  if(!Unit){
107  RANGE_CHECK(A()(n, n) != value_type());//matrix is singular
108  b()(n) /= A()(n, n);
109  }
110  if (b()(n) != value_type/*zero*/()) {
111  matrix_column<TriangularA const> col = column(A(),n);
112  noalias(subrange(b(),0,n)) -= b()(n) * subrange(col,0,n);
113  }
114  }
115 }
116 //upper triangular(row-major)-vector
117 template<bool Unit, class TriangularA, class V>
118 void trsv_impl(
119  matrix_expression<TriangularA> const& A,
120  vector_expression<V> &b,
121  boost::mpl::true_, row_major
122 ) {
123  SIZE_CHECK(A().size1() == A().size2());
124  SIZE_CHECK(A().size2() == b().size());
125 
126  typedef typename TriangularA::value_type value_type;
127 
128  std::size_t size = A().size1();
129  for (std::size_t i = 0; i < size; ++ i) {
130  std::size_t n = size-i-1;
131  matrix_row<TriangularA const> matRow = row(A(),n);
132  b()(n) -= inner_prod(subrange(matRow,n+1,size),subrange(b(),n+1,size));
133  if(!Unit){
134  RANGE_CHECK(A()(n, n) != value_type());//matrix is singular
135  b()(n) /= A()(n, n);
136  }
137  }
138 }
139 
140 //dispatcher
141 
142 template <bool Upper,bool Unit,typename TriangularA, typename V>
143 void trsv(
144  matrix_expression<TriangularA> const& A,
145  vector_expression<V> & b,
146  boost::mpl::false_//unoptimized
147 ){
148  trsv_impl<Unit>(A, b, boost::mpl::bool_<Upper>(), typename TriangularA::orientation());
149 }
150 
151 }}}
152 #endif