trmv.hpp
Go to the documentation of this file.
1 //===========================================================================
2 /*!
3  *
4  * \brief -
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_TRMV_HPP
32 #define SHARK_LINALG_BLAS_KERNELS_DEFAULT_TRMV_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
41 template<bool Unit, class TriangularA, class V>
42 void trmv_impl(
43  matrix_expression<TriangularA> const& A,
44  vector_expression<V> &b,
45  boost::mpl::false_, row_major
46 ){
47  typedef typename TriangularA::value_type value_typeA;
48  typedef typename V::value_type value_typeV;
49  std::size_t size = A().size1();
50  std::size_t const blockSize = 128;
51  std::size_t numBlocks = size/blockSize;
52  if(numBlocks*blockSize < size) ++numBlocks;
53 
54  //this implementation partitions A into
55  //a set of panels, where a Panel is a set
56  // of columns. We start with the last panel
57  //and compute the product of it with the part of the vector
58  // and than just add the previous panel on top etc.
59 
60  //tmporary storage for subblocks of b
61  value_typeV valueStorage[blockSize];
62 
63  for(std::size_t bi = 1; bi <= numBlocks; ++bi){
64  std::size_t startbi = blockSize*(numBlocks-bi);
65  std::size_t sizebi = std::min(blockSize,size-startbi);
66  dense_vector_adaptor<value_typeA> values(valueStorage,sizebi);
67 
68  //store and save the values of b we ar enow changing
69  noalias(values) = subrange(b,startbi,startbi+sizebi);
70 
71  //multiply with triangular element
72  for (std::size_t i = 0; i != sizebi; ++i) {
73  std::size_t posi = startbi+i;
74  b()(posi) = 0;
75  for(std::size_t j = 0; j < i; ++j){
76  b()(posi) += A()(posi,startbi+j)*values(j);
77  }
78  b()(posi) += values(i)*(Unit? value_typeA(1):A()(posi,posi));
79  }
80  //now compute the remaining inner products
81  for(std::size_t posi = startbi+sizebi; posi != size; ++posi){
82  b()(posi) += inner_prod(values,subrange(row(A,posi),startbi,startbi+sizebi));
83  }
84  }
85 }
86 
87 //upper triangular(row-major)-vector
88 template<bool Unit, class TriangularA, class V>
89 void trmv_impl(
90  matrix_expression<TriangularA> const& A,
91  vector_expression<V>& b,
92  boost::mpl::true_, row_major
93 ){
94  std::size_t size = A().size1();
95  for (std::size_t i = 0; i < size; ++ i) {
96  if(!Unit){
97  b()(i) *= A()(i,i);
98  }
99  b()(i) += inner_prod(subrange(row(A,i),i+1,size),subrange(b,i+1,size));
100  }
101 }
102 
103 //Lower triangular(column-major) - vector
104 template<bool Unit, class TriangularA, class V>
105 void trmv_impl(
106  matrix_expression<TriangularA> const& A,
107  vector_expression<V> &b,
108  boost::mpl::false_, column_major
109 ){
110 
111  std::size_t size = A().size1();
112  for (std::size_t n = 1; n <= size; ++n) {
113  std::size_t i = size-n;
114  double bi = b()(i);
115  if(!Unit){
116  b()(i) *= A()(i,i);
117  }
118  noalias(subrange(b,i+1,size))+= bi * subrange(column(A,i),i+1,size);
119  }
120 }
121 
122 //upper triangular(column-major)-vector
123 template<bool Unit, class TriangularA, class V>
124 void trmv_impl(
125  matrix_expression<TriangularA> const& A,
126  vector_expression<V>& b,
127  boost::mpl::true_, column_major
128 ){
129  typedef typename TriangularA::value_type value_typeA;
130  typedef typename V::value_type value_typeV;
131  std::size_t size = A().size1();
132  std::size_t const blockSize = 128;
133  std::size_t numBlocks = size/blockSize;
134  if(numBlocks*blockSize < size) ++numBlocks;
135 
136  //this implementation partitions A into
137  //a set of panels, where a Panel is a set
138  // of rows. We start with the first panel
139  //and compute the product of it with the part of the vector
140  // and than just add the next panel on top etc.
141 
142  //tmporary storage for subblocks of b
143  value_typeV valueStorage[blockSize];
144 
145  for(std::size_t bj = 0; bj != numBlocks; ++bj){
146  std::size_t startbj = blockSize*bj;
147  std::size_t sizebj = std::min(blockSize,size-startbj);
148  dense_vector_adaptor<value_typeA> values(valueStorage,sizebj);
149 
150  //store and save the values of b we ar enow changing
151  noalias(values) = subrange(b,startbj,startbj+sizebj);
152  subrange(b,startbj,startbj+sizebj).clear();
153  //multiply with triangular element
154  for (std::size_t j = 0; j != sizebj; ++j) {
155  std::size_t posj = startbj+j;
156  for(std::size_t i = 0; i < j; ++i){
157  b()(startbj+i) += A()(startbj+i,posj)*values(j);
158  }
159  b()(posj) += values(j)*(Unit? 1.0:A()(posj,posj));
160  }
161  //now compute the remaining inner products
162  for(std::size_t posj = startbj+sizebj; posj != size; ++posj){
163  noalias(subrange(b,startbj,startbj+sizebj)) += b()(posj)*subrange(column(A,posj),startbj,startbj+sizebj);
164  }
165  }
166 }
167 
168 //dispatcher
169 template <bool Upper,bool Unit,typename TriangularA, typename V>
170 void trmv(
171  matrix_expression<TriangularA> const& A,
172  vector_expression<V> & b,
173  boost::mpl::false_//unoptimized
174 ){
175  SIZE_CHECK(A().size1() == A().size2());
176  SIZE_CHECK(A().size2() == b().size());
177  trmv_impl<Unit>(A, b, boost::mpl::bool_<Upper>(), typename TriangularA::orientation());
178 }
179 
180 }}}
181 #endif