syev.hpp
Go to the documentation of this file.
1 //===========================================================================
2 /*!
3  *
4  *
5  * \brief Contains the lapack bindings for the symmetric eigenvalue problem syev.
6  *
7  * \author O. Krause
8  * \date 2010
9  *
10  *
11  * \par Copyright 1995-2015 Shark Development Team
12  *
13  * <BR><HR>
14  * This file is part of Shark.
15  * <http://image.diku.dk/shark/>
16  *
17  * Shark is free software: you can redistribute it and/or modify
18  * it under the terms of the GNU Lesser General Public License as published
19  * by the Free Software Foundation, either version 3 of the License, or
20  * (at your option) any later version.
21  *
22  * Shark is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU Lesser General Public License for more details.
26  *
27  * You should have received a copy of the GNU Lesser General Public License
28  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
29  *
30  */
31 //===========================================================================
32 #ifndef SHARK_LINALG_BLAS_KERNELS_LAPACK_SYEV_HPP
33 #define SHARK_LINALG_BLAS_KERNELS_LAPACK_SYEV_HPP
34 
35 #include "fortran.hpp"
36 #include "../traits.hpp"
37 
38 #define SHARK_LAPACK_DSYEV FORTRAN_ID(dsyev)
39 
40 extern "C"{
41 void SHARK_LAPACK_DSYEV(
42  const char* jobz, const char* uplo, const int *n,
43  double* a, const int * lda, double* w,
44  double* work, const int * lwork, int* info
45 );
46 }
47 
48 
49 
50 namespace shark { namespace blas { namespace bindings {
51 
52 inline void syev(
53  int n, bool upper,
54  double* A, int lda,
55  double* eigenvalues
56 ){
57  if(n == 0) return;
58  int lwork = std::min<int>(130,4*n)*n;
59  double* work = new double[lwork];
60  int info;
61  char job = 'V';
62  char uplo = upper?'U':'L';
63  SHARK_LAPACK_DSYEV(&job, &uplo, &n, A, &lda,eigenvalues,work,&lwork,&info);
64  delete[] work;
65 
66 }
67 
68 
69 template <typename MatrA, typename VectorB>
70 void syev(
71  matrix_expression<MatrA>& matA,
72  vector_expression<VectorB>& eigenValues
73 ) {
74  SIZE_CHECK(matA().size1() == matA().size2());
75  SIZE_CHECK(matA().size1() == eigenValues().size());
76 
77  std::size_t n = matA().size1();
78  bool upper = false;
79  //lapack is column major storage.
80  if(boost::is_same<typename MatrA::orientation, blas::row_major>::value){
81  upper = !upper;
82  }
83  syev(
84  n, upper,
85  traits::storage(matA()),
86  traits::leading_dimension(matA()),
87  traits::storage(eigenValues())
88  );
89 
90  matA() = trans(matA);
91 
92  //reverse eigenvectors and eigenvalues
93  for (int i = 0; i < (int)n-i-1; i++)
94  {
95  int l = n-i-1;
96  std::swap(eigenValues()( l ),eigenValues()( i ));
97  }
98  for (int j = 0; j < (int)n; j++) {
99  for (int i = 0; i < (int)n-i-1; i++)
100  {
101  int l = n-i-1;
102  std::swap(matA()( j , l ), matA()( j , i ));
103  }
104  }
105 }
106 
107 }}}
108 
109 #undef SHARK_LAPACK_DSYEV
110 
111 #endif