ViennaCL - The Vienna Computing Library  1.2.0
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-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 
26 #include <utility>
27 #include <iostream>
28 #include <fstream>
29 #include <string>
30 #include <algorithm>
31 #include <vector>
32 #include <math.h>
33 #include <map>
34 //#include "spai-dynamic.hpp"
35 #include "boost/numeric/ublas/vector.hpp"
36 #include "boost/numeric/ublas/matrix.hpp"
37 #include "boost/numeric/ublas/matrix_proxy.hpp"
38 #include "boost/numeric/ublas/vector_proxy.hpp"
39 #include "boost/numeric/ublas/storage.hpp"
40 #include "boost/numeric/ublas/io.hpp"
41 #include "boost/numeric/ublas/lu.hpp"
42 #include "boost/numeric/ublas/triangular.hpp"
43 #include "boost/numeric/ublas/matrix_expression.hpp"
44 // ViennaCL includes
45 #include "viennacl/linalg/prod.hpp"
46 #include "viennacl/matrix.hpp"
50 #include "viennacl/scalar.hpp"
51 #include "viennacl/linalg/cg.hpp"
53 #include "viennacl/linalg/ilu.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 
66  /********************************* STATIC SPAI FUNCTIONS******************************************/
67 
73  template <typename VectorType, typename SparseVectorType>
74  void fanOutVector(const VectorType& m_in, const std::vector<unsigned int>& J, SparseVectorType& m){
75  unsigned int cnt = 0;
76  for (size_t i = 0; i < J.size(); ++i) {
77  m[J[i]] = m_in(cnt++);
78  }
79  }
85  template <typename MatrixType, typename VectorType>
86  void backwardSolve(const MatrixType& R, const VectorType& y, VectorType& x){
87  typedef typename MatrixType::value_type ScalarType;
88  for (long i = R.size2()-1; i >= 0 ; i--) {
89  x(i) = y(i);
90  for (size_t j = i+1; j < R.size2(); ++j) {
91  x(i) -= R(i,j)*x(j);
92  }
93  x(i) /= R(i,i);
94  }
95  }
101  template <typename VectorType, typename ScalarType>
102  void projectI(const std::vector<unsigned int>& I, VectorType& y, unsigned int ind){
103  for(size_t i = 0; i < I.size(); ++i){
104  //y.resize(y.size()+1);
105  if(I[i] == ind){
106  y(i) = static_cast<ScalarType>(1.0);
107  }
108  else{
109  y(i) = static_cast<ScalarType>(0.0);
110  }
111  }
112  }
113 
118  template <typename SparseVectorType>
119  void buildColumnIndexSet(const SparseVectorType& v, std::vector<unsigned int>& J){
120  //typedef typename VectorType::value_type ScalarType;
121  unsigned int tmp_v;
122  for(typename SparseVectorType::const_iterator vec_it = v.begin(); vec_it != v.end(); ++vec_it){
123  tmp_v = vec_it->first;
124  J.push_back(vec_it->first);
125  }
126  std::sort(J.begin(), J.end());
127  }
128 
133  template <typename SparseMatrixType>
134  void initPreconditioner(const SparseMatrixType& A, SparseMatrixType& M){
135  typedef typename SparseMatrixType::value_type ScalarType;
136  M.resize(A.size1(), A.size2(), false);
137  for(typename SparseMatrixType::const_iterator1 row_it = A.begin1(); row_it!= A.end1(); ++row_it){
138  //
139  for(typename SparseMatrixType::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it){
140  M(col_it.index1(),col_it.index2()) = static_cast<ScalarType>(1);
141  }
142  }
143  }
144 
150  template <typename SparseVectorType>
151  void projectRows(const std::vector<SparseVectorType>& A_v_c, const std::vector<unsigned int>& J, std::vector<unsigned int>& I){
152  for(size_t i = 0; i < J.size(); ++i){
153  for(typename SparseVectorType::const_iterator col_it = A_v_c[J[i]].begin(); col_it!=A_v_c[J[i]].end(); ++col_it){
154  if(!isInIndexSet(I, col_it->first)){
155  I.push_back(col_it->first);
156  }
157  }
158  }
159  std::sort(I.begin(), I.end());
160  }
161  }
162  }
163  }
164 }
165 
166 #endif