30 #ifndef SHARK_LINALG_BLAS_KERNELS_DEFAULT_POTRF_HPP 31 #define SHARK_LINALG_BLAS_KERNELS_DEFAULT_POTRF_HPP 33 #include "../../matrix_proxy.hpp" 34 #include "../../vector_expression.hpp" 35 #include <boost/mpl/bool.hpp> 44 std::size_t potrf_impl(
45 matrix_expression<MatA>& A,
48 std::size_t m = A().size1();
49 for(
size_t j = 0; j < m; j++) {
50 for(
size_t i = j; i < m; i++) {
52 for(
size_t k = 0; k < j; k++) {
53 s -= A()(i, k) * A()(j, k);
58 A()(i, j) = std::sqrt(s);
60 A()(i, j) = s / A()(j , j);
69 std::size_t potrf_impl(
70 matrix_expression<MatA>& A,
73 std::size_t m = A().size1();
74 for(
size_t i = 0; i < m; i++) {
75 double& Aii = A()(i, i);
81 for(std::size_t j = i + 1; j < m; ++j) {
85 for(
size_t k = i + 1; k < m; k++) {
86 for(std::size_t j = k; j < m; ++j) {
87 A()(k, j) -= A()(i, k) * A()(i, j);
96 template<
class MatA,
class Triangular>
97 std::size_t potrf_impl(
98 matrix_container<MatA>& A,
99 column_major, Triangular
101 blas::matrix_transpose<MatA> transA(A());
102 return potrf_impl(transA, row_major(),
typename Triangular::transposed_orientation());
107 template <
class Triangular,
typename MatA>
109 matrix_container<MatA>& A,
112 return potrf_impl(A,
typename MatA::orientation(), Triangular());