1 #ifndef VIENNACL_LINALG_DETAIL_SPAI_SMALL_MATRIX_HPP
2 #define VIENNACL_LINALG_DETAIL_SPAI_SMALL_MATRIX_HPP
34 #include "boost/numeric/ublas/vector.hpp"
35 #include "boost/numeric/ublas/matrix.hpp"
36 #include "boost/numeric/ublas/matrix_proxy.hpp"
37 #include "boost/numeric/ublas/vector_proxy.hpp"
38 #include "boost/numeric/ublas/storage.hpp"
39 #include "boost/numeric/ublas/io.hpp"
40 #include "boost/numeric/ublas/lu.hpp"
41 #include "boost/numeric/ublas/triangular.hpp"
42 #include "boost/numeric/ublas/matrix_expression.hpp"
43 #include "boost/numeric/ublas/detail/matrix_assign.hpp"
59 template <
typename MatrixType>
62 mat.resize(new_size, new_size,
false);
65 double val = 1 / sqrt(2.0);
67 for (
size_t i=0; i<new_size; ++i)
70 for (
size_t i=off_diagonal_distance; i<new_size; ++i)
72 mat(i-off_diagonal_distance, i) = val; mat(i, i-off_diagonal_distance) = -val;
79 template <
typename MatrixType>
80 double determinant(boost::numeric::ublas::matrix_expression<MatrixType>
const& mat_r)
84 MatrixType mLu(mat_r() );
85 boost::numeric::ublas::permutation_matrix<std::size_t> pivots(mat_r().
size1() );
87 int is_singular =
static_cast<int>(
lu_factorize(mLu, pivots));
91 for (std::size_t i=0; i < pivots.size(); ++i)