1 #ifndef VIENNACL_LINALG_DETAIL_SPAI_FSPAI_HPP
2 #define VIENNACL_LINALG_DETAIL_SPAI_FSPAI_HPP
30 #include "boost/numeric/ublas/vector.hpp"
31 #include "boost/numeric/ublas/matrix.hpp"
32 #include "boost/numeric/ublas/matrix_proxy.hpp"
33 #include "boost/numeric/ublas/vector_proxy.hpp"
34 #include "boost/numeric/ublas/storage.hpp"
35 #include "boost/numeric/ublas/io.hpp"
36 #include "boost/numeric/ublas/lu.hpp"
37 #include "boost/numeric/ublas/triangular.hpp"
38 #include "boost/numeric/ublas/matrix_expression.hpp"
78 double residual_norm_threshold = 1e-3,
79 unsigned int iteration_limit = 5,
80 bool is_static =
false,
81 bool is_right =
false) :
82 _residual_norm_threshold(residual_norm_threshold),
83 _iteration_limit(iteration_limit),
84 _is_static(is_static),
85 _is_right(is_right){};
88 {
return _residual_norm_threshold; }
90 {
return _iteration_limit; }
92 {
return _is_static; }
96 if(residual_norm_threshold > 0)
97 _residual_norm_threshold = residual_norm_threshold;
100 if(iteration_limit > 0)
101 _iteration_limit = iteration_limit;
104 _is_right = is_right;
107 _is_static = is_static;
111 double _residual_norm_threshold;
112 unsigned long _iteration_limit;
122 template <
typename MatrixType,
typename ScalarType>
125 STL_A.resize(A.size1());
126 for (
typename MatrixType::const_iterator1 row_it = A.begin1();
130 for (
typename MatrixType::const_iterator2 col_it = row_it.begin();
131 col_it != row_it.end();
134 if (col_it.index1() >= col_it.index2())
135 STL_A[col_it.index1()][col_it.index2()] = *col_it;
146 template <
typename MatrixType>
147 void generateJ(MatrixType
const & A, std::vector<std::vector<size_t> > & J)
149 for (
typename MatrixType::const_iterator1 row_it = A.begin1();
153 for (
typename MatrixType::const_iterator2 col_it = row_it.begin();
154 col_it != row_it.end();
157 if (col_it.index1() > col_it.index2())
159 J[col_it.index2()].push_back(col_it.index1());
160 J[col_it.index1()].push_back(col_it.index2());
173 template <
typename ScalarType,
typename MatrixType,
typename VectorType>
174 void fill_blocks(std::vector< std::map<unsigned int, ScalarType> > & A,
175 std::vector<MatrixType> & blocks,
176 std::vector<std::vector<size_t> >
const & J,
177 std::vector<VectorType> & Y)
179 for (
size_t k=0; k<A.size(); ++k)
181 std::vector<size_t>
const & Jk = J[k];
182 VectorType & yk = Y[k];
183 MatrixType & block_k = blocks[k];
185 yk.resize(Jk.size());
186 block_k.resize(Jk.size(), Jk.size());
189 for (
size_t i=0; i<Jk.size(); ++i)
191 size_t row_index = Jk[i];
192 std::map<unsigned int, ScalarType> & A_row = A[row_index];
197 for (
size_t j=0; j<Jk.size(); ++j)
199 size_t col_index = Jk[j];
200 if (col_index <= row_index && A_row.find(col_index) != A_row.end())
201 block_k(i, j) = A_row[col_index];
211 template <
typename MatrixType>
214 for (
size_t k=0; k<A.size2(); ++k)
218 std::cout <<
"k: " << k << std::endl;
219 std::cout <<
"A(k,k): " << A(k,k) << std::endl;
224 A(k,k) = std::sqrt(A(k,k));
226 for (
size_t i=k+1; i<A.size1(); ++i)
229 for (
size_t j=k+1; j<=i; ++j)
230 A(i,j) -= A(i,k) * A(j,k);
239 template <
typename MatrixType,
typename VectorType>
242 typedef typename VectorType::value_type ScalarType;
245 for (
size_t i=0; i<L.size1(); ++i)
247 for (
size_t j=0; j<i; ++j)
248 b[i] -= L(i,j) * b[j];
253 for (
size_t i=L.size1()-1; ; --i)
255 for (
size_t k=i+1; k<L.size1(); ++k)
256 b[i] -= L(k,i) * b[k];
269 template <
typename MatrixType,
typename VectorType1>
272 MatrixType & L_trans,
273 std::vector<VectorType1> & Y,
274 std::vector<std::vector<size_t> > & J)
276 typedef typename VectorType1::value_type ScalarType;
277 typedef std::vector<std::map<unsigned int, ScalarType> > STLSparseMatrixType;
279 STLSparseMatrixType L_temp(A.size1());
281 for (
size_t k=0; k<A.size1(); ++k)
283 std::vector<size_t>
const & Jk = J[k];
284 VectorType1
const & yk = Y[k];
287 ScalarType Lkk = A(k,k);
288 for (
size_t i=0; i<Jk.size(); ++i)
289 Lkk -= A(Jk[i],k) * yk[i];
291 Lkk = 1.0 / sqrt(Lkk);
296 for (
size_t i=0; i<Jk.size(); ++i)
298 L_temp[Jk[i]][k] = -Lkk * yk[i];
299 L_trans(k, Jk[i]) = -Lkk * yk[i];
305 for (
size_t i=0; i<L_temp.size(); ++i)
306 for (
typename std::map<unsigned int, ScalarType>::const_iterator it = L_temp[i].begin();
307 it != L_temp[i].end();
309 L(i, it->first) = it->second;
316 template <
typename MatrixType>
318 MatrixType
const & PatternA,
320 MatrixType & L_trans,
323 typedef typename MatrixType::value_type ScalarType;
324 typedef boost::numeric::ublas::matrix<ScalarType> DenseMatrixType;
325 typedef std::vector<std::map<unsigned int, ScalarType> > SparseMatrixType;
331 std::vector<std::vector<ScalarType> > y_k(A.size1());
332 SparseMatrixType STL_A(A.size1());
340 std::vector<std::vector<size_t> > J(A.size1());
347 std::vector<DenseMatrixType> subblocks_A(A.size1());
355 for (
size_t i=0; i<subblocks_A.size(); ++i)
372 for (
size_t i=0; i<y_k.size(); ++i)
374 if (subblocks_A[i].
size1() > 0)
390 L.resize(A.size1(), A.size2(),
false);
391 L.reserve(A.nnz(),
false);
392 L_trans.resize(A.size1(), A.size2(),
false);
393 L_trans.reserve(A.nnz(),
false);