ViennaCL - The Vienna Computing Library  1.2.0
small_matrix.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_DETAIL_SPAI_SMALL_MATRIX_HPP
2 #define VIENNACL_LINALG_DETAIL_SPAI_SMALL_MATRIX_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 "boost/numeric/ublas/vector.hpp"
35 #include "boost/numeric/ublas/matrix.hpp"
36 #include "boost/numeric/ublas/matrix_proxy.hpp"
37 #include "boost/numeric/ublas/vector_proxy.hpp"
38 #include "boost/numeric/ublas/storage.hpp"
39 #include "boost/numeric/ublas/io.hpp"
40 #include "boost/numeric/ublas/lu.hpp"
41 #include "boost/numeric/ublas/triangular.hpp"
42 #include "boost/numeric/ublas/matrix_expression.hpp"
43 #include "boost/numeric/ublas/detail/matrix_assign.hpp"
44 
45 
46 
47 namespace viennacl
48 {
49  namespace linalg
50  {
51  namespace detail
52  {
53  namespace spai
54  {
55 
56  //
57  // Constructs an orthonormal sparse matrix M (with M^T M = Id). Is composed of elementary 2x2 rotation matrices with suitable renumbering.
58  //
59  template <typename MatrixType>
60  void make_rotation_matrix(MatrixType & mat, size_t new_size, size_t off_diagonal_distance = 4)
61  {
62  mat.resize(new_size, new_size, false);
63  mat.clear();
64 
65  double val = 1 / sqrt(2.0);
66 
67  for (size_t i=0; i<new_size; ++i)
68  mat(i,i) = val;
69 
70  for (size_t i=off_diagonal_distance; i<new_size; ++i)
71  {
72  mat(i-off_diagonal_distance, i) = val; mat(i, i-off_diagonal_distance) = -val;
73  }
74 
75  }
76 
77 
78  //calcualtes matrix determinant
79  template <typename MatrixType>
80  double determinant(boost::numeric::ublas::matrix_expression<MatrixType> const& mat_r)
81  {
82  double det = 1.0;
83 
84  MatrixType mLu(mat_r() );
85  boost::numeric::ublas::permutation_matrix<std::size_t> pivots(mat_r().size1() );
86 
87  int is_singular = static_cast<int>(lu_factorize(mLu, pivots));
88 
89  if (!is_singular)
90  {
91  for (std::size_t i=0; i < pivots.size(); ++i)
92  {
93  if (pivots(i) != i)
94  det *= -1.0;
95 
96  det *= mLu(i,i);
97  }
98  }
99  else
100  det = 0.0;
101 
102  return det;
103  }
104 
105  }
106  }
107  }
108 }
109 #endif