fml  0.1-0
Fused Matrix Library
eigen.hh
1 // This file is part of fml which is released under the Boost Software
2 // License, Version 1.0. See accompanying file LICENSE or copy at
3 // https://www.boost.org/LICENSE_1_0.txt
4 
5 #ifndef FML_CPU_LINALG_LINALG_EIGEN_H
6 #define FML_CPU_LINALG_LINALG_EIGEN_H
7 #pragma once
8 
9 
10 #include <cmath>
11 #include <stdexcept>
12 
13 #include "../../_internals/linalgutils.hh"
14 #include "../../_internals/omp.hh"
15 
16 #include "../cpumat.hh"
17 #include "../cpuvec.hh"
18 
19 #include "internals/lapack.hh"
20 
21 
22 namespace fml
23 {
24 namespace linalg
25 {
26  namespace
27  {
28  template <typename REAL>
29  int eig_sym_internals(const bool only_values, cpumat<REAL> &x,
30  cpuvec<REAL> &values, cpumat<REAL> &vectors)
31  {
32  if (!x.is_square())
33  throw std::runtime_error("'x' must be a square matrix");
34 
35  int info = 0;
36  int nfound;
37  char jobz;
38 
39  len_t n = x.nrows();
40  values.resize(n);
41  cpuvec<int> support;
42 
43  if (only_values)
44  jobz = 'N';
45  else
46  {
47  jobz = 'V';
48  vectors.resize(n, n);
49  support.resize(2*n);
50  }
51 
52  REAL worksize;
53  int lwork, liwork;
54  fml::lapack::syevr(jobz, 'A', 'L', n, x.data_ptr(), n, (REAL) 0.f, (REAL) 0.f,
55  0, 0, (REAL) 0.f, &nfound, values.data_ptr(), vectors.data_ptr(), n,
56  support.data_ptr(), &worksize, -1, &liwork, -1,
57  &info);
58 
59  lwork = (int) worksize;
60  cpuvec<REAL> work(lwork);
61  cpuvec<int> iwork(liwork);
62 
63  fml::lapack::syevr(jobz, 'A', 'L', n, x.data_ptr(), n, (REAL) 0.f, (REAL) 0.f,
64  0, 0, (REAL) 0.f, &nfound, values.data_ptr(), vectors.data_ptr(), n,
65  support.data_ptr(), work.data_ptr(), lwork, iwork.data_ptr(), liwork,
66  &info);
67 
68  return info;
69  }
70  }
71 
72 
73 
94  template <typename REAL>
96  {
97  cpumat<REAL> ignored;
98 
99  int info = eig_sym_internals(true, x, values, ignored);
100  fml::linalgutils::check_info(info, "syevr");
101  }
102 
103 
104 
106  template <typename REAL>
107  void eigen_sym(cpumat<REAL> &x, cpuvec<REAL> &values, cpumat<REAL> &vectors)
108  {
109  int info = eig_sym_internals(false, x, values, vectors);
110  fml::linalgutils::check_info(info, "syevr");
111  }
112 }
113 }
114 
115 
116 #endif
fml::cpumat
Matrix class for data held on a single CPU.
Definition: cpumat.hh:36
fml::cpuvec
Vector class for data held on a single CPU.
Definition: cpuvec.hh:31
fml
Core namespace.
Definition: dimops.hh:10
fml::linalg::eigen_sym
void eigen_sym(cpumat< REAL > &x, cpuvec< REAL > &values)
Compute the eigenvalues and optionally the eigenvectors for a symmetric matrix.
Definition: eigen.hh:95