ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
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-2015, 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 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 
71 template<typename SizeT>
72 bool isInIndexSet(std::vector<SizeT> const & J, SizeT ind)
73 {
74  return (std::find(J.begin(), J.end(), ind) != J.end());
75 }
76 
77 
78 
79 /********************************* STATIC SPAI FUNCTIONS******************************************/
80 
87 template<typename VectorT, typename SparseVectorT>
88 void fanOutVector(VectorT const & m_in, std::vector<unsigned int> const & J, SparseVectorT & m)
89 {
90  unsigned int cnt = 0;
91  for (vcl_size_t i = 0; i < J.size(); ++i)
92  m[J[i]] = m_in(cnt++);
93 }
94 
101 template<typename MatrixT, typename VectorT>
102 void backwardSolve(MatrixT const & R, VectorT const & y, VectorT & x)
103 {
104  for (long i2 = static_cast<long>(R.size2())-1; i2 >= 0; i2--)
105  {
106  vcl_size_t i = static_cast<vcl_size_t>(i2);
107  x(i) = y(i);
108  for (vcl_size_t j = static_cast<vcl_size_t>(i)+1; j < R.size2(); ++j)
109  x(i) -= R(i,j)*x(j);
110 
111  x(i) /= R(i,i);
112  }
113 }
114 
121 template<typename VectorT, typename NumericT>
122 void projectI(std::vector<unsigned int> const & I, VectorT & y, unsigned int ind)
123 {
124  for (vcl_size_t i = 0; i < I.size(); ++i)
125  {
126  //y.resize(y.size()+1);
127  if (I[i] == ind)
128  y(i) = NumericT(1.0);
129  else
130  y(i) = NumericT(0.0);
131  }
132 }
133 
139 template<typename SparseVectorT>
140 void buildColumnIndexSet(SparseVectorT const & v, std::vector<unsigned int> & J)
141 {
142  for (typename SparseVectorT::const_iterator vec_it = v.begin(); vec_it != v.end(); ++vec_it)
143  J.push_back(vec_it->first);
144 
145  std::sort(J.begin(), J.end());
146 }
147 
153 template<typename SparseMatrixT>
154 void initPreconditioner(SparseMatrixT const & A, SparseMatrixT & M)
155 {
156  typedef typename SparseMatrixT::value_type NumericType;
157 
158  M.resize(A.size1(), A.size2(), false);
159  for (typename SparseMatrixT::const_iterator1 row_it = A.begin1(); row_it!= A.end1(); ++row_it)
160  for (typename SparseMatrixT::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it)
161  M(col_it.index1(),col_it.index2()) = NumericType(1);
162 }
163 
170 template<typename SparseVectorT>
171 void projectRows(std::vector<SparseVectorT> const & A_v_c,
172  std::vector<unsigned int> const & J,
173  std::vector<unsigned int> & I)
174 {
175  for (vcl_size_t i = 0; i < J.size(); ++i)
176  {
177  for (typename SparseVectorT::const_iterator col_it = A_v_c[J[i]].begin(); col_it!=A_v_c[J[i]].end(); ++col_it)
178  {
179  if (!isInIndexSet(I, col_it->first))
180  I.push_back(col_it->first);
181  }
182  }
183  std::sort(I.begin(), I.end());
184 }
185 
186 
187 } //namespace spai
188 } //namespace detail
189 } //namespace linalg
190 } //namespace viennacl
191 
192 #endif
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.
Generic interface for the computation of inner products. See viennacl/linalg/vector_operations.hpp for implementations.
float NumericT
Definition: bisect.cpp:40
bool isInIndexSet(std::vector< SizeT > const &J, SizeT ind)
Determines if element ind is in set {J}.
Definition: spai-static.hpp:72
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:34
void projectRows(std::vector< SparseVectorT > const &A_v_c, std::vector< unsigned int > const &J, std::vector< unsigned int > &I)
Row projection for matrix A(:,J) -> A(I,J), building index set of non-zero rows.
void backwardSolve(MatrixT const &R, VectorT const &y, VectorT &x)
Solution of linear:R*x=y system by backward substitution.
Definition: blas3.hpp:36
Implementation of the compressed_matrix class.
Implementations of operations using sparse matrices.
void projectI(std::vector< unsigned int > const &I, VectorT &y, unsigned int ind)
Perform projection of set I on the unit-vector.
std::size_t vcl_size_t
Definition: forwards.h:75
The conjugate gradient method is implemented here.
void fanOutVector(VectorT const &m_in, std::vector< unsigned int > const &J, SparseVectorT &m)
Projects solution of LS problem onto original column m.
Definition: spai-static.hpp:88
void buildColumnIndexSet(SparseVectorT const &v, std::vector< unsigned int > &J)
Builds index set of projected columns for current column of preconditioner.
void initPreconditioner(SparseMatrixT const &A, SparseMatrixT &M)
Initialize preconditioner with sparcity pattern = p(A)
Implementation of the ViennaCL scalar class.