1 #ifndef VIENNACL_LINALG_DETAIL_SPAI_FSPAI_HPP
2 #define VIENNACL_LINALG_DETAIL_SPAI_FSPAI_HPP
31 #include "boost/numeric/ublas/vector.hpp"
32 #include "boost/numeric/ublas/matrix.hpp"
33 #include "boost/numeric/ublas/matrix_proxy.hpp"
34 #include "boost/numeric/ublas/vector_proxy.hpp"
35 #include "boost/numeric/ublas/storage.hpp"
36 #include "boost/numeric/ublas/io.hpp"
37 #include "boost/numeric/ublas/lu.hpp"
38 #include "boost/numeric/ublas/triangular.hpp"
39 #include "boost/numeric/ublas/matrix_expression.hpp"
79 double residual_norm_threshold = 1e-3,
80 unsigned int iteration_limit = 5,
81 bool is_static =
false,
82 bool is_right =
false) :
83 residual_norm_threshold_(residual_norm_threshold),
84 iteration_limit_(iteration_limit),
85 is_static_(is_static),
86 is_right_(is_right) {}
89 {
return residual_norm_threshold_; }
91 {
return iteration_limit_; }
93 {
return is_static_; }
97 if(residual_norm_threshold > 0)
98 residual_norm_threshold_ = residual_norm_threshold;
101 if(iteration_limit > 0)
102 iteration_limit_ = iteration_limit;
105 is_right_ = is_right;
108 is_static_ = is_static;
112 double residual_norm_threshold_;
113 unsigned long iteration_limit_;
123 template <
typename MatrixType,
typename ScalarType>
126 STL_A.resize(A.size1());
127 for (
typename MatrixType::const_iterator1 row_it = A.begin1();
131 for (
typename MatrixType::const_iterator2 col_it = row_it.begin();
132 col_it != row_it.end();
135 if (col_it.index1() >= col_it.index2())
136 STL_A[col_it.index1()][
static_cast<unsigned int>(col_it.index2())] = *col_it;
147 template <
typename MatrixType>
148 void generateJ(MatrixType
const & A, std::vector<std::vector<vcl_size_t> > & J)
150 for (
typename MatrixType::const_iterator1 row_it = A.begin1();
154 for (
typename MatrixType::const_iterator2 col_it = row_it.begin();
155 col_it != row_it.end();
158 if (col_it.index1() > col_it.index2())
160 J[col_it.index2()].push_back(col_it.index1());
161 J[col_it.index1()].push_back(col_it.index2());
174 template <
typename ScalarType,
typename MatrixType,
typename VectorType>
175 void fill_blocks(std::vector< std::map<unsigned int, ScalarType> > & A,
176 std::vector<MatrixType> & blocks,
177 std::vector<std::vector<vcl_size_t> >
const & J,
178 std::vector<VectorType> & Y)
182 std::vector<vcl_size_t>
const & Jk = J[k];
183 VectorType & yk = Y[k];
184 MatrixType & block_k = blocks[k];
186 yk.resize(Jk.size());
187 block_k.resize(Jk.size(), Jk.size());
193 std::map<unsigned int, ScalarType> & A_row = A[row_index];
196 yk[i] = A_row[
static_cast<unsigned int>(k)];
201 if (col_index <= row_index && A_row.find(static_cast<unsigned int>(col_index)) != A_row.end())
202 block_k(i, j) = A_row[
static_cast<unsigned int>(col_index)];
212 template <
typename MatrixType>
219 std::cout <<
"k: " << k << std::endl;
220 std::cout <<
"A(k,k): " << A(k,k) << std::endl;
225 A(k,k) = std::sqrt(A(k,k));
231 A(i,j) -= A(i,k) * A(j,k);
240 template <
typename MatrixType,
typename VectorType>
247 b[i] -= L(i,j) * b[j];
255 b[i] -= L(k,i) * b[k];
268 template <
typename MatrixType,
typename VectorType1>
271 MatrixType & L_trans,
272 std::vector<VectorType1> & Y,
273 std::vector<std::vector<vcl_size_t> > & J)
275 typedef typename VectorType1::value_type ScalarType;
276 typedef std::vector<std::map<unsigned int, ScalarType> > STLSparseMatrixType;
278 STLSparseMatrixType L_temp(A.size1());
282 std::vector<vcl_size_t>
const & Jk = J[k];
283 VectorType1
const & yk = Y[k];
286 ScalarType Lkk = A(k,k);
288 Lkk -= A(Jk[i],k) * yk[i];
290 Lkk = ScalarType(1) / std::sqrt(Lkk);
291 L_temp[k][
static_cast<unsigned int>(k)] = Lkk;
297 L_temp[Jk[i]][
static_cast<unsigned int>(k)] = -Lkk * yk[i];
298 L_trans(k, Jk[i]) = -Lkk * yk[i];
305 for (
typename std::map<unsigned int, ScalarType>::const_iterator it = L_temp[i].begin();
306 it != L_temp[i].end();
308 L(i, it->first) = it->second;
315 template <
typename MatrixType>
317 MatrixType
const & PatternA,
319 MatrixType & L_trans,
322 typedef typename MatrixType::value_type ScalarType;
323 typedef boost::numeric::ublas::matrix<ScalarType> DenseMatrixType;
324 typedef std::vector<std::map<unsigned int, ScalarType> > SparseMatrixType;
330 std::vector<std::vector<ScalarType> > y_k(A.size1());
331 SparseMatrixType STL_A(A.size1());
339 std::vector<std::vector<vcl_size_t> > J(A.size1());
346 std::vector<DenseMatrixType> subblocks_A(A.size1());
354 for (
vcl_size_t i=0; i<subblocks_A.size(); ++i)
373 if (subblocks_A[i].
size1() > 0)
389 L.resize(A.size1(), A.size2(),
false);
390 L.reserve(A.nnz(),
false);
391 L_trans.resize(A.size1(), A.size2(),
false);
392 L_trans.reserve(A.nnz(),
false);
void setIsRight(bool is_right)
Definition: fspai.hpp:104
unsigned long getIterationLimit() const
Definition: fspai.hpp:90
std::size_t vcl_size_t
Definition: forwards.h:58
Implementations of dense matrix related operations including matrix-vector products.
fspai_tag(double residual_norm_threshold=1e-3, unsigned int iteration_limit=5, bool is_static=false, bool is_right=false)
Constructor.
Definition: fspai.hpp:78
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
Implementation of the dense matrix class.
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:216
void setResidualNormThreshold(double residual_norm_threshold)
Definition: fspai.hpp:96
void generateJ(MatrixType const &A, std::vector< std::vector< vcl_size_t > > &J)
Definition: fspai.hpp:148
Generic interface for the computation of inner products. See viennacl/linalg/vector_operations.hpp for implementations.
void sym_sparse_matrix_to_stl(MatrixType const &A, std::vector< std::map< unsigned int, ScalarType > > &STL_A)
Definition: fspai.hpp:124
void computeFSPAI(MatrixType const &A, MatrixType const &PatternA, MatrixType &L, MatrixType &L_trans, fspai_tag)
Definition: fspai.hpp:316
A tag for FSPAI. Experimental. Contains values for the algorithm. Must be passed to spai_precond cons...
Definition: fspai.hpp:70
double getResidualNormThreshold() const
Definition: fspai.hpp:88
Implementations of incomplete factorization preconditioners. Convenience header file.
Implementation of the compressed_matrix class.
bool getIsRight() const
Definition: fspai.hpp:94
Implementations of operations using sparse matrices.
void cholesky_decompose(MatrixType &A)
Definition: fspai.hpp:213
void cholesky_solve(MatrixType const &L, VectorType &b)
Definition: fspai.hpp:241
The conjugate gradient method is implemented here.
void computeL(MatrixType const &A, MatrixType &L, MatrixType &L_trans, std::vector< VectorType1 > &Y, std::vector< std::vector< vcl_size_t > > &J)
Definition: fspai.hpp:269
void fill_blocks(std::vector< std::map< unsigned int, ScalarType > > &A, std::vector< MatrixType > &blocks, std::vector< std::vector< vcl_size_t > > const &J, std::vector< VectorType > &Y)
Definition: fspai.hpp:175
void setIterationLimit(unsigned long iteration_limit)
Definition: fspai.hpp:100
bool getIsStatic() const
Definition: fspai.hpp:92
void setIsStatic(bool is_static)
Definition: fspai.hpp:107
Implementation of the ViennaCL scalar class.