31 #ifndef SHARK_LINALG_BLAS_KERNELS_ATLAS_TRSM_HPP 32 #define SHARK_LINALG_BLAS_KERNELS_ATLAS_TRSM_HPP 34 #include "../../matrix_proxy.hpp" 35 #include "../../vector_expression.hpp" 36 #include <boost/mpl/bool.hpp> 38 namespace shark {
namespace blas {
namespace bindings {
43 template<
bool Unit,
class MatA,
class MatB>
45 matrix_expression<MatA>
const& A, matrix_expression<MatB>& B,
46 boost::mpl::false_, column_major
51 typedef typename MatA::value_type value_type;
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) {
60 B()(n, l) /= A()(n, n);
62 if (B()(n, l) != value_type()) {
63 matrix_column<MatB> columnMatrix =
column(B(),l);
70 template<
bool Unit,
class MatA,
class MatB>
72 matrix_expression<MatA>
const& A, matrix_expression<MatB>& B,
73 boost::mpl::false_, row_major
78 typedef typename MatA::value_type value_type;
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) {
87 row(B(),n)/=A()(n, n);
93 template<
bool Unit,
class MatA,
class MatB>
95 matrix_expression<MatA>
const& A, matrix_expression<MatB>& B,
96 boost::mpl::true_, column_major
101 typedef typename MatA::value_type value_type;
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);
110 row(B(),n) /= A()(n, n);
112 for (std::size_t l = 0; l < size2; ++ l) {
113 if (B()(n, l) != value_type()) {
114 matrix_column<MatB> columnMatrix =
column(B(),l);
122 template<
bool Unit,
class MatA,
class MatB>
124 matrix_expression<MatA>
const& A, matrix_expression<MatB>& B,
125 boost::mpl::true_, row_major
130 typedef typename MatA::value_type value_type;
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) {
140 row(B(),n)/=A()(n, n);
145 template <
bool Upper,
bool Unit,
typename TriangularA,
typename MatB>
147 matrix_expression<TriangularA>
const& A,
148 matrix_expression<MatB>& B,
153 boost::mpl::bool_<Upper>(),
154 typename TriangularA::orientation()