ViennaCL - The Vienna Computing Library  1.5.2
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-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 "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 #include "boost/numeric/ublas/detail/matrix_assign.hpp"
45 
46 
47 
48 namespace viennacl
49 {
50  namespace linalg
51  {
52  namespace detail
53  {
54  namespace spai
55  {
56 
57  //
58  // Constructs an orthonormal sparse matrix M (with M^T M = Id). Is composed of elementary 2x2 rotation matrices with suitable renumbering.
59  //
60  template <typename MatrixType>
61  void make_rotation_matrix(MatrixType & mat, vcl_size_t new_size, vcl_size_t off_diagonal_distance = 4)
62  {
63  mat.resize(new_size, new_size, false);
64  mat.clear();
65 
66  double val = 1.0 / std::sqrt(2.0);
67 
68  for (vcl_size_t i=0; i<new_size; ++i)
69  mat(i,i) = val;
70 
71  for (vcl_size_t i=off_diagonal_distance; i<new_size; ++i)
72  {
73  mat(i-off_diagonal_distance, i) = val; mat(i, i-off_diagonal_distance) = -val;
74  }
75 
76  }
77 
78 
79  //calcualtes matrix determinant
80  template <typename MatrixType>
81  double determinant(boost::numeric::ublas::matrix_expression<MatrixType> const& mat_r)
82  {
83  double det = 1.0;
84 
85  MatrixType mLu(mat_r() );
86  boost::numeric::ublas::permutation_matrix<vcl_size_t> pivots(mat_r().size1() );
87 
88  int is_singular = static_cast<int>(lu_factorize(mLu, pivots));
89 
90  if (!is_singular)
91  {
92  for (vcl_size_t i=0; i < pivots.size(); ++i)
93  {
94  if (pivots(i) != i)
95  det *= -1.0;
96 
97  det *= mLu(i,i);
98  }
99  }
100  else
101  det = 0.0;
102 
103  return det;
104  }
105 
106  }
107  }
108  }
109 }
110 #endif
std::size_t vcl_size_t
Definition: forwards.h:58
void make_rotation_matrix(MatrixType &mat, vcl_size_t new_size, vcl_size_t off_diagonal_distance=4)
Definition: small_matrix.hpp:61
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:216
void lu_factorize(matrix< SCALARTYPE, viennacl::row_major > &A)
LU factorization of a row-major dense matrix.
Definition: lu.hpp:42
double determinant(boost::numeric::ublas::matrix_expression< MatrixType > const &mat_r)
Definition: small_matrix.hpp:81