fml  0.1-0
Fused Matrix Library
linalg_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_MPI_LINALG_LINALG_EIGEN_H
6 #define FML_MPI_LINALG_LINALG_EIGEN_H
7 #pragma once
8 
9 
10 #include <stdexcept>
11 
12 #include "../../_internals/linalgutils.hh"
13 #include "../../cpu/cpuvec.hh"
14 
15 #include "../internals/bcutils.hh"
16 #include "../internals/mpi_utils.hh"
17 
18 #include "../copy.hh"
19 #include "../mpimat.hh"
20 
21 #include "linalg_blas.hh"
22 #include "linalg_err.hh"
23 #include "scalapack.hh"
24 
25 
26 namespace fml
27 {
28 namespace linalg
29 {
30  namespace
31  {
32  template <typename REAL>
33  int eig_sym_internals(const bool only_values, mpimat<REAL> &x,
34  cpuvec<REAL> &values, mpimat<REAL> &vectors)
35  {
36  if (!x.is_square())
37  throw std::runtime_error("'x' must be a square matrix");
38 
39  int info = 0;
40  int val_found, vec_found;
41  char jobz;
42 
43  len_t n = x.nrows();
44  values.resize(n);
45 
46  if (only_values)
47  jobz = 'N';
48  else
49  {
50  jobz = 'V';
51  vectors.resize(n, n, x.bf_rows(), x.bf_cols());
52  }
53 
54  REAL worksize;
55  int lwork, liwork;
56 
57  fml::scalapack::syevr(jobz, 'A', 'L', n, x.data_ptr(), x.desc_ptr(),
58  (REAL) 0.f, (REAL) 0.f, 0, 0, &val_found, &vec_found,
59  values.data_ptr(), vectors.data_ptr(), vectors.desc_ptr(),
60  &worksize, -1, &liwork, -1, &info);
61 
62  lwork = (int) worksize;
63  cpuvec<REAL> work(lwork);
64  cpuvec<int> iwork(liwork);
65 
66  fml::scalapack::syevr(jobz, 'A', 'L', n, x.data_ptr(), x.desc_ptr(),
67  (REAL) 0.f, (REAL) 0.f, 0, 0, &val_found, &vec_found,
68  values.data_ptr(), vectors.data_ptr(), vectors.desc_ptr(),
69  work.data_ptr(), lwork, iwork.data_ptr(), liwork, &info);
70 
71  return info;
72  }
73  }
74 
97  template <typename REAL>
99  {
100  mpimat<REAL> ignored(x.get_grid());
101 
102  int info = eig_sym_internals(true, x, values, ignored);
103  fml::linalgutils::check_info(info, "syevr");
104  }
105 
107  template <typename REAL>
108  void eigen_sym(mpimat<REAL> &x, cpuvec<REAL> &values, mpimat<REAL> &vectors)
109  {
110  err::check_grid(x, vectors);
111 
112  int info = eig_sym_internals(false, x, values, vectors);
113  fml::linalgutils::check_info(info, "syevr");
114  }
115 }
116 }
117 
118 
119 #endif
fml::mpimat
Matrix class for data distributed over MPI in the 2-d block cyclic format.
Definition: mpimat.hh:40
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: linalg_eigen.hh:93