1 #ifndef VIENNACL_LINALG_DETAIL_SPAI_SPAI_STATIC_HPP
2 #define VIENNACL_LINALG_DETAIL_SPAI_SPAI_STATIC_HPP
35 #include "boost/numeric/ublas/vector.hpp"
36 #include "boost/numeric/ublas/matrix.hpp"
37 #include "boost/numeric/ublas/matrix_proxy.hpp"
38 #include "boost/numeric/ublas/vector_proxy.hpp"
39 #include "boost/numeric/ublas/storage.hpp"
40 #include "boost/numeric/ublas/io.hpp"
41 #include "boost/numeric/ublas/lu.hpp"
42 #include "boost/numeric/ublas/triangular.hpp"
43 #include "boost/numeric/ublas/matrix_expression.hpp"
73 template <
typename VectorType,
typename SparseVectorType>
74 void fanOutVector(
const VectorType& m_in,
const std::vector<unsigned int>& J, SparseVectorType& m){
76 for (
size_t i = 0; i < J.size(); ++i) {
77 m[J[i]] = m_in(cnt++);
85 template <
typename MatrixType,
typename VectorType>
86 void backwardSolve(
const MatrixType& R,
const VectorType& y, VectorType& x){
87 typedef typename MatrixType::value_type ScalarType;
88 for (
long i = R.size2()-1; i >= 0 ; i--) {
90 for (
size_t j = i+1; j < R.size2(); ++j) {
101 template <
typename VectorType,
typename ScalarType>
102 void projectI(
const std::vector<unsigned int>& I, VectorType& y,
unsigned int ind){
103 for(
size_t i = 0; i < I.size(); ++i){
106 y(i) =
static_cast<ScalarType
>(1.0);
109 y(i) =
static_cast<ScalarType
>(0.0);
118 template <
typename SparseVectorType>
122 for(
typename SparseVectorType::const_iterator vec_it = v.begin(); vec_it != v.end(); ++vec_it){
123 tmp_v = vec_it->first;
124 J.push_back(vec_it->first);
126 std::sort(J.begin(), J.end());
133 template <
typename SparseMatrixType>
135 typedef typename SparseMatrixType::value_type ScalarType;
136 M.resize(A.size1(), A.size2(),
false);
137 for(
typename SparseMatrixType::const_iterator1 row_it = A.begin1(); row_it!= A.end1(); ++row_it){
139 for(
typename SparseMatrixType::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it){
140 M(col_it.index1(),col_it.index2()) = static_cast<ScalarType>(1);
150 template <
typename SparseVectorType>
151 void projectRows(
const std::vector<SparseVectorType>& A_v_c,
const std::vector<unsigned int>& J, std::vector<unsigned int>& I){
152 for(
size_t i = 0; i < J.size(); ++i){
153 for(
typename SparseVectorType::const_iterator col_it = A_v_c[J[i]].begin(); col_it!=A_v_c[J[i]].end(); ++col_it){
155 I.push_back(col_it->first);
159 std::sort(I.begin(), I.end());