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 REMORA_KERNELS_LAPACK_SYEV_HPP
33 #define REMORA_KERNELS_LAPACK_SYEV_HPP
34 
35 #include "fortran.hpp"
36 #include "../../detail/traits.hpp"
37 
38 #define REMORA_LAPACK_DSYEV FORTRAN_ID(dsyev)
39 
40 extern "C"{
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 remora {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  REMORA_LAPACK_DSYEV(&job, &uplo, &n, A, &lda,eigenvalues,work,&lwork,&info);
64  delete[] work;
65 
66 }
67 
68 
69 template <typename MatA, typename VectorB>
70 void syev(
71  matrix_expression<MatA, cpu_tag>& A,
72  vector_expression<VectorB, cpu_tag>& eigenValues
73 ) {
74  REMORA_SIZE_CHECK(A().size1() == A().size2());
75  REMORA_SIZE_CHECK(A().size1() == eigenValues().size());
76 
77  std::size_t n = A().size1();
78  bool upper = false;
79  //lapack is column major storage.
80  if(std::is_same<typename MatA::orientation, row_major>::value){
81  upper = !upper;
82  }
83  auto storageA = A().raw_storage();
84  auto storageEig = eigenValues().raw_storage();
85  syev(
86  n, upper,
87  storageA.values,
88  storageA.leading_dimension,
89  storageEig.values
90  );
91 
92  A() = trans(A);
93 
94  //reverse eigenvectors and eigenvalues
95  for (int i = 0; i < (int)n-i-1; i++)
96  {
97  int l = n-i-1;
98  std::swap(eigenValues()( l ),eigenValues()( i ));
99  }
100  for (int j = 0; j < (int)n; j++) {
101  for (int i = 0; i < (int)n-i-1; i++)
102  {
103  int l = n-i-1;
104  std::swap(A()( j , l ), A()( j , i ));
105  }
106  }
107 }
108 
109 }}
110 
111 #undef REMORA_LAPACK_DSYEV
112 
113 #endif