trsm.hpp
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief -
5  *
6  * \author O. Krause
7  * \date 2012
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_ATLAS_TRSM_HPP
32 #define SHARK_LINALG_BLAS_KERNELS_ATLAS_TRSM_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 //fixme: no handling of orientation of MatB. Might be solved by rewriting with block-algorithms.
41 
42 // Lower triangular(column major) - matrix
43 template<bool Unit, class MatA, class MatB>
44 void trsm_impl(
45  matrix_expression<MatA> const& A, matrix_expression<MatB>& B,
46  boost::mpl::false_, column_major
47 ) {
48  SIZE_CHECK(A().size1() == A().size2());
49  SIZE_CHECK(A().size2() == B().size1());
50 
51  typedef typename MatA::value_type value_type;
52 
53  std::size_t size1 = B().size1();
54  std::size_t size2 = B().size2();
55  for (std::size_t n = 0; n < size1; ++ n) {
56  matrix_column<MatA const> columnTriangular = column(A(),n);
57  for (std::size_t l = 0; l < size2; ++ l) {
58  if(!Unit){
59  RANGE_CHECK(A()(n, n) != value_type());//matrix is singular
60  B()(n, l) /= A()(n, n);
61  }
62  if (B()(n, l) != value_type/*zero*/()) {
63  matrix_column<MatB> columnMatrix = column(B(),l);
64  noalias(subrange(columnMatrix,n+1,size1)) -= B()(n,l) * subrange(columnTriangular,n+1,size1);
65  }
66  }
67  }
68 }
69 // Lower triangular(row major) - matrix
70 template<bool Unit, class MatA, class MatB>
71 void trsm_impl(
72  matrix_expression<MatA> const& A, matrix_expression<MatB>& B,
73  boost::mpl::false_, row_major
74 ) {
75  SIZE_CHECK(A().size1() == A().size2());
76  SIZE_CHECK(A().size2() == B().size1());
77 
78  typedef typename MatA::value_type value_type;
79 
80  std::size_t size1 = B().size1();
81  for (std::size_t n = 0; n < size1; ++ n) {
82  for (std::size_t m = 0; m < n; ++m) {
83  noalias(row(B(),n)) -= A()(n,m)*row(B(),m);
84  }
85  if(!Unit){
86  RANGE_CHECK(A()(n, n) != value_type());//matrix is singular
87  row(B(),n)/=A()(n, n);
88  }
89  }
90 }
91 
92 //Upper triangular(column major) - matrix
93 template<bool Unit, class MatA, class MatB>
94 void trsm_impl(
95  matrix_expression<MatA> const& A, matrix_expression<MatB>& B,
96  boost::mpl::true_, column_major
97 ) {
98  SIZE_CHECK(A().size1() == A().size2());
99  SIZE_CHECK(A().size2() == B().size1());
100 
101  typedef typename MatA::value_type value_type;
102 
103  std::size_t size1 = B().size1();
104  std::size_t size2 = B().size2();
105  for (std::size_t i = 0; i < size1; ++ i) {
106  std::size_t n = size1-i-1;
107  matrix_column<MatA const> columnTriangular = column(A(),n);
108  if(!Unit){
109  RANGE_CHECK(A()(n, n) != value_type());//matrix is singular
110  row(B(),n) /= A()(n, n);
111  }
112  for (std::size_t l = 0; l < size2; ++ l) {
113  if (B()(n, l) != value_type/*zero*/()) {
114  matrix_column<MatB> columnMatrix = column(B(),l);
115  noalias(subrange(columnMatrix,0,n)) -= B()(n,l) * subrange(columnTriangular,0,n);
116  }
117  }
118  }
119 }
120 
121 //Upper triangular(row major) - matrix
122 template<bool Unit, class MatA, class MatB>
123 void trsm_impl(
124  matrix_expression<MatA> const& A, matrix_expression<MatB>& B,
125  boost::mpl::true_, row_major
126 ) {
127  SIZE_CHECK(A().size1() == A().size2());
128  SIZE_CHECK(A().size2() == B().size1());
129 
130  typedef typename MatA::value_type value_type;
131 
132  std::size_t size1 = B().size1();
133  for (std::size_t i = 0; i < size1; ++ i) {
134  std::size_t n = size1-i-1;
135  for (std::size_t m = n+1; m < size1; ++m) {
136  noalias(row(B(),n)) -= A()(n,m)*row(B(),m);
137  }
138  if(!Unit){
139  RANGE_CHECK(A()(n, n) != value_type());//matrix is singular
140  row(B(),n)/=A()(n, n);
141  }
142  }
143 }
144 
145 template <bool Upper, bool Unit,typename TriangularA, typename MatB>
146 void trsm(
147  matrix_expression<TriangularA> const& A,
148  matrix_expression<MatB>& B,
149  boost::mpl::false_
150 ){
151  trsm_impl<Unit>(
152  A,B,
153  boost::mpl::bool_<Upper>(),
154  typename TriangularA::orientation()
155  );
156 }
157 
158 }}}
159 #endif