tpmv.hpp
Go to the documentation of this file.
1 //===========================================================================
2 /*!
3  *
4  * \brief Implements the triangular packed matrix-vector multiplication
5  *
6  * \author O. Krause
7  * \date 2010
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 //===========================================================================
31 #ifndef SHARK_LINALG_BLAS_KERNELS_DEFAULT_tpmv_HPP
32 #define SHARK_LINALG_BLAS_KERNELS_DEFAULT_tpmv_HPP
33 
34 #include "../../matrix_proxy.hpp"
35 #include "../../vector_expression.hpp"
36 #include <boost/mpl/bool.hpp>
37 
38 namespace shark{ namespace blas{ namespace bindings{
39 
40 //Lower triangular(row-major) - vector product
41 // computes the row-wise inner product of the matrix
42 // starting with row 0 and stores the result in b(i)
43 //this does not interfere with the next products as
44 // the element b(i) is not needed for iterations j > i
45 template<class TriangularA, class V>
46 void tpmv_impl(
47  matrix_expression<TriangularA> const& A,
48  vector_expression<V>& b,
49  row_major,
50  upper
51 ){
52  typedef typename V::value_type value_type;
53  typedef typename V::size_type size_type;
54  typedef typename TriangularA::const_row_iterator row_iterator;
55  size_type size = A().size1();
56 
57  for(size_type i = 0; i != size; ++i){
58  value_type sum(0);
59  row_iterator end = A().row_end(i);
60  for(row_iterator pos = A().row_begin(i); pos != end; ++pos){
61  sum += *pos*b()(pos.index());
62  }
63  b()(i) = sum;
64  }
65 }
66 
67 //Upper triangular(row-major) - vector product
68 // computes the row-wise inner product of the matrix
69 // starting with the last row and stores the result in b(i)
70 // this does not interfere with the next products as
71 // the element b(i) is not needed for row products j < i
72 template<class TriangularA, class V>
73 void tpmv_impl(
74  matrix_expression<TriangularA> const& A,
75  vector_expression<V>& b,
76  row_major,
77  lower
78 ){
79  typedef typename V::value_type value_type;
80  typedef typename V::size_type size_type;
81  typedef typename TriangularA::const_row_iterator row_iterator;
82  size_type size = A().size1();
83 
84  for(size_type irev = size; irev != 0; --irev){
85  size_type i= irev-1;
86  value_type sum(0);
87  row_iterator end = A().row_end(i);
88  for(row_iterator pos = A().row_begin(i); pos != end; ++pos){
89  sum += *pos*b()(pos.index());
90  }
91  b()(i) = sum;
92  }
93 }
94 
95 //Upper triangular(column-major) - vector product
96 //computes the product as a series of vector-additions
97 //on b starting with the last column of A.
98 template<class TriangularA, class V>
99 void tpmv_impl(
100  matrix_expression<TriangularA> const& A,
101  vector_expression<V>& b,
102  column_major,
103  upper
104 ){
105  typedef typename TriangularA::const_column_iterator column_iterator;
106  typedef typename V::size_type size_type;
107  typedef typename V::value_type value_type;
108  size_type size = A().size1();
109  for(size_type i = 0; i != size; ++i){
110  value_type bi = b()(i);
111  b()(i)= value_type/*zero*/();
112  column_iterator end = A().column_end(i);
113  for(column_iterator pos = A().column_begin(i); pos != end; ++pos){
114  b()(pos.index()) += *pos*bi;
115  }
116  }
117 
118 }
119 
120 //Lower triangular(column-major) - vector product
121 // computes the product as a series of vector-additions
122 // on b starting with the first column of A.
123 template<class TriangularA, class V>
124 void tpmv_impl(
125  matrix_expression<TriangularA> const& A,
126  vector_expression<V>& b,
127  column_major,
128  lower
129 ){
130  typedef typename TriangularA::const_column_iterator column_iterator;
131  typedef typename V::size_type size_type;
132  typedef typename V::value_type value_type;
133  size_type size = A().size1();
134 
135  for(size_type irev = size; irev != 0; --irev){
136  size_type i= irev-1;
137  value_type bi = b()(i);
138  b()(i)= value_type/*zero*/();
139  column_iterator end = A().column_end(i);
140  for(column_iterator pos = A().column_begin(i); pos != end; ++pos){
141  b()(pos.index()) += *pos*bi;
142  }
143  }
144 }
145 
146 //dispatcher
147 template <typename TriangularA, typename V>
148 void tpmv(
149  matrix_expression<TriangularA> const& A,
150  vector_expression<V>& b,
151  boost::mpl::false_//unoptimized
152 ){
153  SIZE_CHECK(A().size1() == A().size2());
154  SIZE_CHECK(A().size2() == b().size());
155  tpmv_impl(
156  A, b,
157  typename TriangularA::orientation::orientation(),
158  typename TriangularA::orientation::triangular_type()
159  );
160 }
161 
162 }}}
163 #endif