ViennaCL - The Vienna Computing Library  1.5.1
spai-static.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_DETAIL_SPAI_SPAI_STATIC_HPP
2 #define VIENNACL_LINALG_DETAIL_SPAI_SPAI_STATIC_HPP
3 
4 /* =========================================================================
5  Copyright (c) 2010-2014, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8  Portions of this software are copyright by UChicago Argonne, LLC.
9 
10  -----------------
11  ViennaCL - The Vienna Computing Library
12  -----------------
13 
14  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
15 
16  (A list of authors and contributors can be found in the PDF manual)
17 
18  License: MIT (X11), see file LICENSE in the base directory
19 ============================================================================= */
20 
27 #include <utility>
28 #include <iostream>
29 #include <fstream>
30 #include <string>
31 #include <algorithm>
32 #include <vector>
33 #include <math.h>
34 #include <map>
35 //#include "spai-dynamic.hpp"
36 #include "boost/numeric/ublas/vector.hpp"
37 #include "boost/numeric/ublas/matrix.hpp"
38 #include "boost/numeric/ublas/matrix_proxy.hpp"
39 #include "boost/numeric/ublas/vector_proxy.hpp"
40 #include "boost/numeric/ublas/storage.hpp"
41 #include "boost/numeric/ublas/io.hpp"
42 #include "boost/numeric/ublas/lu.hpp"
43 #include "boost/numeric/ublas/triangular.hpp"
44 #include "boost/numeric/ublas/matrix_expression.hpp"
45 // ViennaCL includes
46 #include "viennacl/linalg/prod.hpp"
47 #include "viennacl/matrix.hpp"
51 #include "viennacl/scalar.hpp"
52 #include "viennacl/linalg/cg.hpp"
54 
55 //#include "boost/numeric/ublas/detail/matrix_assign.hpp"
56 
57 namespace viennacl
58 {
59  namespace linalg
60  {
61  namespace detail
62  {
63  namespace spai
64  {
65 
70  template <typename SizeType>
71  bool isInIndexSet(const std::vector<SizeType>& J, SizeType ind)
72  {
73  return (std::find(J.begin(), J.end(), ind) != J.end());
74  }
75 
76 
77 
78  /********************************* STATIC SPAI FUNCTIONS******************************************/
79 
85  template <typename VectorType, typename SparseVectorType>
86  void fanOutVector(const VectorType& m_in, const std::vector<unsigned int>& J, SparseVectorType& m)
87  {
88  unsigned int cnt = 0;
89  for (vcl_size_t i = 0; i < J.size(); ++i)
90  m[J[i]] = m_in(cnt++);
91  }
97  template <typename MatrixType, typename VectorType>
98  void backwardSolve(const MatrixType& R, const VectorType& y, VectorType& x)
99  {
100  for (long i = static_cast<long>(R.size2())-1; i >= 0 ; i--)
101  {
102  x(i) = y(i);
103  for (vcl_size_t j = i+1; j < R.size2(); ++j)
104  x(i) -= R(i,j)*x(j);
105 
106  x(i) /= R(i,i);
107  }
108  }
114  template <typename VectorType, typename ScalarType>
115  void projectI(const std::vector<unsigned int>& I, VectorType& y, unsigned int ind)
116  {
117  for(vcl_size_t i = 0; i < I.size(); ++i)
118  {
119  //y.resize(y.size()+1);
120  if(I[i] == ind)
121  y(i) = static_cast<ScalarType>(1.0);
122  else
123  y(i) = static_cast<ScalarType>(0.0);
124  }
125  }
126 
131  template <typename SparseVectorType>
132  void buildColumnIndexSet(const SparseVectorType& v, std::vector<unsigned int>& J)
133  {
134  //typedef typename VectorType::value_type ScalarType;
135  //unsigned int tmp_v;
136  for(typename SparseVectorType::const_iterator vec_it = v.begin(); vec_it != v.end(); ++vec_it)
137  {
138  //tmp_v = vec_it->first;
139  J.push_back(vec_it->first);
140  }
141  std::sort(J.begin(), J.end());
142  }
143 
148  template <typename SparseMatrixType>
149  void initPreconditioner(const SparseMatrixType& A, SparseMatrixType& M)
150  {
151  typedef typename SparseMatrixType::value_type ScalarType;
152  M.resize(A.size1(), A.size2(), false);
153  for(typename SparseMatrixType::const_iterator1 row_it = A.begin1(); row_it!= A.end1(); ++row_it)
154  {
155  //
156  for(typename SparseMatrixType::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it)
157  {
158  M(col_it.index1(),col_it.index2()) = static_cast<ScalarType>(1);
159  }
160  }
161  }
162 
168  template <typename SparseVectorType>
169  void projectRows(const std::vector<SparseVectorType>& A_v_c, const std::vector<unsigned int>& J, std::vector<unsigned int>& I)
170  {
171  for(vcl_size_t i = 0; i < J.size(); ++i)
172  {
173  for(typename SparseVectorType::const_iterator col_it = A_v_c[J[i]].begin(); col_it!=A_v_c[J[i]].end(); ++col_it)
174  {
175  if(!isInIndexSet(I, col_it->first))
176  I.push_back(col_it->first);
177  }
178  }
179  std::sort(I.begin(), I.end());
180  }
181 
182 
183  } //namespace spai
184  } //namespace detail
185  } //namespace linalg
186 } //namespace viennacl
187 
188 #endif
std::size_t vcl_size_t
Definition: forwards.h:58
Implementations of dense matrix related operations including matrix-vector products.
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
Implementation of the dense matrix class.
void backwardSolve(const MatrixType &R, const VectorType &y, VectorType &x)
Solution of linear:R*x=y system by backward substitution.
Definition: spai-static.hpp:98
bool isInIndexSet(const std::vector< SizeType > &J, SizeType ind)
Determines if element ind is in set {J}.
Definition: spai-static.hpp:71
void projectRows(const std::vector< SparseVectorType > &A_v_c, const std::vector< unsigned int > &J, std::vector< unsigned int > &I)
Row projection for matrix A(:,J) -> A(I,J), building index set of non-zero rows.
Definition: spai-static.hpp:169
Generic interface for the computation of inner products. See viennacl/linalg/vector_operations.hpp for implementations.
Implementation of the compressed_matrix class.
Implementations of operations using sparse matrices.
void fanOutVector(const VectorType &m_in, const std::vector< unsigned int > &J, SparseVectorType &m)
Projects solution of LS problem onto original column m.
Definition: spai-static.hpp:86
The conjugate gradient method is implemented here.
void buildColumnIndexSet(const SparseVectorType &v, std::vector< unsigned int > &J)
Builds index set of projected columns for current column of preconditioner.
Definition: spai-static.hpp:132
void initPreconditioner(const SparseMatrixType &A, SparseMatrixType &M)
Initialize preconditioner with sparcity pattern = p(A)
Definition: spai-static.hpp:149
void projectI(const std::vector< unsigned int > &I, VectorType &y, unsigned int ind)
Perform projection of set I on the unit-vector.
Definition: spai-static.hpp:115
Implementation of the ViennaCL scalar class.