ViennaCL - The Vienna Computing Library  1.2.0
spai.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_SPAI_HPP
2 #define VIENNACL_LINALG_SPAI_HPP
3 
4 /* =========================================================================
5  Copyright (c) 2010-2011, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8 
9  -----------------
10  ViennaCL - The Vienna Computing Library
11  -----------------
12 
13  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
14 
15  (A list of authors and contributors can be found in the PDF manual)
16 
17  License: MIT (X11), see file LICENSE in the base directory
18 ============================================================================= */
19 
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 
36 //local includes
38 #include "viennacl/linalg/qr.hpp"
46 
47 //boost includes
48 #include "boost/numeric/ublas/vector.hpp"
49 #include "boost/numeric/ublas/matrix.hpp"
50 #include "boost/numeric/ublas/matrix_proxy.hpp"
51 #include "boost/numeric/ublas/vector_proxy.hpp"
52 #include "boost/numeric/ublas/storage.hpp"
53 #include "boost/numeric/ublas/io.hpp"
54 #include "boost/numeric/ublas/lu.hpp"
55 #include "boost/numeric/ublas/triangular.hpp"
56 #include "boost/numeric/ublas/matrix_expression.hpp"
57 
58 // ViennaCL includes
59 #include "viennacl/linalg/prod.hpp"
60 #include "viennacl/matrix.hpp"
64 #include "viennacl/scalar.hpp"
66 #include "viennacl/linalg/ilu.hpp"
67 #include "viennacl/ocl/backend.hpp"
70 
71 
72 namespace viennacl
73 {
74  namespace linalg
75  {
76 
79 
84  //UBLAS version
85  template <typename MatrixType>
87  {
88  public:
89  typedef typename MatrixType::value_type ScalarType;
90  typedef typename boost::numeric::ublas::vector<ScalarType> VectorType;
95  spai_precond(const MatrixType& A,
96  const spai_tag& tag): _tag(tag){
97 
98  //VCLMatrixType vcl_Ap((unsigned int)A.size2(), (unsigned int)A.size1()), vcl_A((unsigned int)A.size1(), (unsigned int)A.size2()),
99  //vcl_At((unsigned int)A.size1(), (unsigned int)A.size2());
100  //UBLASDenseMatrixType dA = A;
101  MatrixType pA(A.size1(), A.size2());
102  MatrixType At;
103  //std::cout<<A<<std::endl;
104  if(!_tag.getIsRight()){
106  }else{
107  At = A;
108  }
109  pA = At;
112  //(At, pA, _tag.getIsRight(), _tag.getIsStatic(), (ScalarType)_tag.getResidualNormThreshold(), (unsigned int)_tag.getIterationLimit(),
113  //_spai_m);
114 
115  }
119  void apply(VectorType& vec) const {
120  vec = viennacl::linalg::prod(_spai_m, vec);
121  }
122  private:
123  // variables
124  spai_tag _tag;
125  // result of SPAI
126  MatrixType _spai_m;
127  };
128 
129  //VIENNACL version
130  template <typename ScalarType, unsigned int MAT_ALIGNMENT>
131  class spai_precond< viennacl::compressed_matrix<ScalarType, MAT_ALIGNMENT> >
132  {
134  typedef boost::numeric::ublas::compressed_matrix<ScalarType> UBLASSparseMatrixType;
137 
138  typedef boost::numeric::ublas::vector<ScalarType> UBLASVectorType;
139  public:
140 
146  const spai_tag& tag): _tag(tag)
147  {
148  viennacl::linalg::kernels::spai<ScalarType, 1>::init();
149 
150  MatrixType At(A.size1(), A.size2());
151  UBLASSparseMatrixType ubls_A, ubls_spai_m;
152  UBLASSparseMatrixType ubls_At;
153  viennacl::copy(A, ubls_A);;
154  if(!_tag.getIsRight()){
156  }
157  else{
158  ubls_At = ubls_A;
159  }
160  //current pattern is A
161  //pA = ubls_At;
162  //execute SPAI with ublas matrix types
164  viennacl::copy(ubls_At, At);
165  viennacl::linalg::detail::spai::computeSPAI(At, ubls_At, ubls_spai_m, _spai_m, _tag);
166  //viennacl::copy(ubls_spai_m, _spai_m);
167 
168  }
172  void apply(VectorType& vec) const {
173  vec = viennacl::linalg::prod(_spai_m, vec);
174  }
175  private:
176  // variables
177  spai_tag _tag;
178  // result of SPAI
179  MatrixType _spai_m;
180  };
181 
182 
183  //
184  // FSPAI
185  //
186 
191  //UBLAS version
192  template <typename MatrixType>
194  {
195  typedef typename MatrixType::value_type ScalarType;
196  typedef typename boost::numeric::ublas::vector<ScalarType> VectorType;
197  typedef typename boost::numeric::ublas::matrix<ScalarType> UBLASDenseMatrixType;
199  public:
200 
205  fspai_precond(const MatrixType& A,
206  const fspai_tag& tag): tag_(tag)
207  {
208  MatrixType pA = A;
209  viennacl::linalg::detail::spai::computeFSPAI(A, pA, L, L_trans, tag_);
210  }
211 
215  void apply(VectorType& vec) const
216  {
217  VectorType temp = viennacl::linalg::prod(L_trans, vec);
218  vec = viennacl::linalg::prod(L, temp);
219  }
220 
221  private:
222  // variables
223  const fspai_tag & tag_;
224  // result of SPAI
225  MatrixType L;
226  MatrixType L_trans;
227  };
228 
229 
230 
231 
232 
233  //
234  // ViennaCL version
235  //
236  template <typename ScalarType, unsigned int MAT_ALIGNMENT>
237  class fspai_precond< viennacl::compressed_matrix<ScalarType, MAT_ALIGNMENT> >
238  {
242  typedef boost::numeric::ublas::compressed_matrix<ScalarType> UBLASSparseMatrixType;
243  typedef boost::numeric::ublas::vector<ScalarType> UBLASVectorType;
244  public:
245 
251  const fspai_tag & tag): tag_(tag){
252  //UBLASSparseMatrixType ubls_A;
253  UBLASSparseMatrixType ublas_A(A.size1(), A.size2());
254  UBLASSparseMatrixType pA(A.size1(), A.size2());
255  UBLASSparseMatrixType ublas_L(A.size1(), A.size2());
256  UBLASSparseMatrixType ublas_L_trans(A.size1(), A.size2());
257  viennacl::copy(A, ublas_A);
258  //viennacl::copy(ubls_A, vcl_A);
259  //vcl_At = viennacl::linalg::prod(vcl_A, vcl_A);
260  //vcl_pA = viennacl::linalg::prod(vcl_A, vcl_At);
261  //viennacl::copy(vcl_pA, pA);
262  pA = ublas_A;
263  //execute SPAI with ublas matrix types
264  viennacl::linalg::detail::spai::computeFSPAI(ublas_A, pA, ublas_L, ublas_L_trans, tag_);
265  //copy back to GPU
266  viennacl::copy(ublas_L, L);
267  viennacl::copy(ublas_L_trans, L_trans);
268  }
269 
270 
274  void apply(VectorType& vec) const
275  {
276  VectorType temp(vec.size());
277  temp = viennacl::linalg::prod(L_trans, vec);
278  vec = viennacl::linalg::prod(L, temp);
279  }
280 
281  private:
282  // variables
283  const fspai_tag & tag_;
284  MatrixType L;
285  MatrixType L_trans;
286  };
287 
288 
289  }
290 }
291 #endif