From 9265c2db7f79983fe3b6d0999bd92a20b0034178 Mon Sep 17 00:00:00 2001
From: Olivier Parcollet <olivier.parcollet@cea.fr>
Date: Wed, 3 Sep 2014 23:18:16 +0200
Subject: [PATCH] Fix #119

- Cleaned of the eigensystems computations (worker is simpler, decision
  at runtime, etc..).
- Fix #119 : When the matrix is in C order, the fortran lapack
  sees in fact its conjugate, so we need to conjugate the eigenvectors at the end.
  NB : not true if the storage order of the matrix is already fortran of course.
---
 test/triqs/arrays/eigenelements.cpp   | 120 ++++++++----
 test/triqs/arrays/tridiag.cpp         |  11 +-
 triqs/arrays/linalg/eigenelements.hpp | 259 ++++++++++++--------------
 triqs/lattice/tight_binding.cpp       |   4 +-
 4 files changed, 209 insertions(+), 185 deletions(-)

diff --git a/test/triqs/arrays/eigenelements.cpp b/test/triqs/arrays/eigenelements.cpp
index f20ab538..14cdc03d 100644
--- a/test/triqs/arrays/eigenelements.cpp
+++ b/test/triqs/arrays/eigenelements.cpp
@@ -19,57 +19,109 @@
  *
  ******************************************************************************/
 #include "./common.hpp"
-#include <triqs/arrays/array.hpp>
-#include <triqs/arrays/vector.hpp>
-#include <triqs/arrays/matrix.hpp>
 #include <triqs/arrays/linalg/eigenelements.hpp>
+#include <triqs/utility/complex_ops.hpp>
 
 #include <iostream>
 
 using namespace triqs::arrays;
 using namespace triqs::arrays::linalg;
+using dcomplex = std::complex<double>;
+
+template <typename T> void check_eig(matrix<T> M, matrix<T> vectors, array<double, 1> values) {
+  auto _ = range();
+ for (auto i : range(0,first_dim(M)))  {
+  std::cerr << "check "<< i << std::endl;
+  std::cerr  << (M  -values(i))* vectors(i, _)<<std::endl;
+  assert_all_close(M * vectors(i, _), values(i) * vectors(i, _), 1.e-14);
+ }}
+
+template<typename M> void test(M A) { 
+  auto w = eigenelements(make_clone(A));
+  std::cerr << "A = " << A << std::endl;
+  std::cerr << " values = " <<w.first << std::endl;
+  std::cerr << " vectors = " << w.second << std::endl;
+  check_eig (A, w.second, w.first);
+}
 
 int main(int argc, char **argv) {
 
- 
+ {
+  matrix<double> A(3, 3);
 
- matrix<double> A(3,3);
+  for (int i = 0; i < 3; ++i)
+   for (int j = 0; j <= i; ++j) {
+    A(i, j) = (i > j ? i + 2 * j : i - j);
+    A(j, i) = A(i, j);
+   }
+  std::cerr << "A = " << A << std::endl;
+  
+  auto B = A;
+  auto w = eigenelements(B);
+  std::cout << "A = " << B << std::endl;
+  std::cout << " vectors = " << w.first << std::endl;
+  std::cout << " values = " << w.second << std::endl;
+  check_eig (A, w.second, w.first);
+ }
 
- for (int i =0; i<3; ++i)
-  for (int j=0; j<=i; ++j)
-  { A(i,j) = (i>j ? i+2*j : i-j); A(j,i) = A(i,j);}
+ {
+  matrix<double> A(3, 3);
+  A() = 0;
+  A(0, 1) = 1;
+  A(1, 0) = 1;
+  A(2, 2) = 8;
+  A(0, 2) = 2;
+  A(2, 0) = 2;
 
- std::cerr<<"A = "<<A<<std::endl;
- eigenelements_worker< matrix_view <double>, true> w (A());
- w.invoke();
- std::cout<<"A = "<<A<<std::endl;
- std::cout<<" vectors = "<< w.values()<<std::endl;
- std::cout<<" values = "<< w.vectors()<<std::endl;
+  auto B = A;
+  std::cout << "A = " << A << std::endl;
+  auto w =  eigenelements(B);
+  std::cout << " values = " <<w.first << std::endl;
+  std::cout << " vectors = " << w.second << std::endl;
+  check_eig (A, w.second, w.first);
+ }
 
- A() =0;
- A(0,1) = 1;
- A(1,0) = 1;
- A(2,2) = 8;
- A(0,2) = 2;
- A(2,0) = 2;
+ {
+  matrix<double> A(3, 3);
 
- std::cout<<"A = "<<A<<std::endl;
- std::cout<<" values = "<< eigenelements(A(),true).first<<std::endl;
- std::cout<<" vectors = "<< eigenelements(A(),true).second<<std::endl;
+  A() = 0;
+  A(0, 1) = 1;
+  A(1, 0) = 1;
+  A(2, 2) = 8;
 
- A() =0;
- A(0,1) = 1;
- A(1,0) = 1;
- A(2,2) = 8;
+  auto B = A;
+  std::cout << "A = " << A << std::endl;
+  auto w =  eigenelements(B);
+  std::cout << " vectors = " << w.second << std::endl;
+  std::cout << " values = " <<w.first << std::endl;
+  std::cout << "A = " << A << std::endl;
+  check_eig (A, w.second, w.first);
+ }
 
- std::cout<<"A = "<<A<<std::endl;
- std::cout<<" vectors = "<< eigenelements(A(),true).second<<std::endl;
- std::cout<<" values = "<< eigenelements(A(),true).first<<std::endl;
- std::cout<<"A = "<<A<<std::endl;
+ { // the complex case
+
+  matrix<dcomplex> M(2, 2);
+
+  M(0, 0) = 1;
+  M(0, 1) = 1.0_j;
+  M(1, 0) = -1.0_j;
+  M(1, 1) = 2;
+
+  test(M);
+ }
+
+{ // the complex case
+
+  matrix<dcomplex> M(2, 2, FORTRAN_LAYOUT);
+
+  M(0, 0) = 1;
+  M(0, 1) = 1.0_j;
+  M(1, 0) = -1.0_j;
+  M(1, 1) = 2;
+
+  test(M);
+ }
 
  return 0;
-
 }
 
-
-
diff --git a/test/triqs/arrays/tridiag.cpp b/test/triqs/arrays/tridiag.cpp
index 7f51b2aa..d4ae8b59 100644
--- a/test/triqs/arrays/tridiag.cpp
+++ b/test/triqs/arrays/tridiag.cpp
@@ -43,13 +43,12 @@ int main(int argc, char **argv) {
  for (size_t i=0; i<4; ++i) { A(i,i+1) = DD(i); A(i+1,i) = DD(i);} 
 
  std::cerr<<"A = "<<A<<std::endl;
- linalg::eigenelements_worker< matrix_view <double>, true> w2 (A());
- w2.invoke();
- std::cerr<<" check values = "<< w2.values()<<std::endl;
- std::cerr<<" check vectors = "<< w2.vectors()<<std::endl;
+ auto eig = linalg::eigenelements(A());
+ std::cerr<<" check values = "<< eig.first<<std::endl;
+ std::cerr<<" check vectors = "<< eig.second<<std::endl;
 
- assert_all_close(w.vectors()(R,R), w2.vectors(), 1.e-10);
- assert_all_close(array_view<double,1>(w.values()), w2.values(), 1.e-10);
+ assert_all_close(w.vectors()(R,R), eig.second, 1.e-10);
+ assert_all_close(array_view<double,1>(w.values()), eig.first, 1.e-10);
 
  return 0;
 }
diff --git a/triqs/arrays/linalg/eigenelements.hpp b/triqs/arrays/linalg/eigenelements.hpp
index 8278cf6a..a08d4c03 100644
--- a/triqs/arrays/linalg/eigenelements.hpp
+++ b/triqs/arrays/linalg/eigenelements.hpp
@@ -2,7 +2,7 @@
  *
  * TRIQS: a Toolbox for Research in Interacting Quantum Systems
  *
- * Copyright (C) 2011 by O. Parcollet
+ * Copyright (C) 2011-2014 by O. Parcollet
  *
  * TRIQS is free software: you can redistribute it and/or modify it under the
  * terms of the GNU General Public License as published by the Free Software
@@ -18,162 +18,135 @@
  * TRIQS. If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef TRIQS_ARRAY_EIGENELEMENTS_H
-#define TRIQS_ARRAY_EIGENELEMENTS_H
+#pragma once
 #include <type_traits>
 #include "../array.hpp"
 #include "../matrix.hpp"
 #include "../vector.hpp"
 #include <triqs/utility/exceptions.hpp>
 
-namespace triqs { namespace arrays { namespace linalg { 
+namespace triqs {
+namespace arrays {
+ namespace linalg {
 
- /**
-  * A worker to call lapack routine with the matrices.
-  *
-  * Handles both real and complex case.
-  *
-  * Usage : 
-  *  - construct from a VIEW of a matrix, that MUST be contiguous.
-  *  - call invoke()
-  *  - read the eigenvalues/vectors in values and vectors resp.
-  *  NB : the content of the matrix is destroyed by the computation (it is .vectors() in fact, if Compute_Eigenvectors is true).
-  *  For a one shot usage, prefer eigenelements, eigenvalues functions.
-  */
- template<typename MatrixViewType, bool Compute_Eigenvectors > struct eigenelements_worker;
- template<typename T, ull_t Opt, bool Compute_Eigenvectors > struct eigenelements_worker_base;
+  extern "C" {
+  void TRIQS_FORTRAN_MANGLING(dsyev)(char *, char *, // JOBZ and UPLO
+                                     int &,          // Matrix Size
+                                     double[],       // matrix
+                                     int &,          // LDA of the matrix
+                                     double[],       // Eigenvalues array
+                                     double[],       // WORK
+                                     int &,          // LWORK
+                                     int &           // INFO
+                                     );
 
- template<typename T, ull_t Opt >
-  struct eigenelements_worker_base <T,Opt,false> { 
-   private:
-    void operator = ( eigenelements_worker_base const & x);
-   protected:
-    matrix_view <T,Opt> mat;
-    triqs::arrays::vector<double> ev; 
-    triqs::arrays::vector<T> work;
-    int dim,info,lwork;
-    char uplo,compz;
-    bool has_run;
 
-    eigenelements_worker_base ( matrix_view <T,Opt> the_matrix) : mat(the_matrix) { 
-     if (mat.is_empty())   TRIQS_RUNTIME_ERROR<<"eigenelements_worker : the matrix is empty : matrix =  "<<mat<<"  ";
-     if (!mat.is_square())   TRIQS_RUNTIME_ERROR<<"eigenelements_worker : the matrix "<<mat<<" is not square ";
-     if (!mat.indexmap().is_contiguous())   TRIQS_RUNTIME_ERROR<<"eigenelements_worker : the matrix "<<mat<<" is not contiguous in memory";
-     dim = first_dim(mat);
-     ev.resize(dim);
-     lwork = 64*dim;
-     work.resize(lwork);
-     uplo='U';compz='N' ;
-     has_run = false;
-    }
+  void TRIQS_FORTRAN_MANGLING(zheev)(char *, char *,          // JOBZ and UPLO
+                                     int &,                   // Matrix Size
+                                     std::complex<double> [], // matrix
+                                     int &,                   // LDA of the matrix
+                                     double[],                // Eigenvalues array
+                                     std::complex<double> [], // WORK
+                                     int &,                   // LWORK
+                                     double[],                // WORK2
+                                     int &                    // INFO
+                                     );
+  }
 
-   public :
-    array<double,1> values() const { 
-     if (!has_run)  TRIQS_RUNTIME_ERROR<<"eigenelements_worker has not been invoked !";
-     return ev;
-    }
-  };
-
- //--------------------------------
-
- template<typename T, ull_t Opt>
-  struct eigenelements_worker_base <T,Opt,true> : eigenelements_worker_base <T,Opt,false>  {
-   protected:
-    eigenelements_worker_base ( matrix_view <T,Opt> the_matrix) :  eigenelements_worker_base <T,Opt,false>  (the_matrix) {this->compz='V'; }
+  /**
+   * A worker to call lapack routine with the matrices. Handles both real and complex case.
+   */
+  template <typename T> class eigenelements_worker {
    public:
-    typename matrix_view<T,Opt>::regular_type vectors() const { 
-     if (!this->has_run)  TRIQS_RUNTIME_ERROR<<"eigenelements_worker has not been invoked !";
-     return this->mat;
-    }
+   eigenelements_worker() = default;
+
+   /// The eigenvalues
+   template <typename M> array<double, 1> eigenvalues(M &mat) const {
+    _prepare(mat);
+    _invoke(is_complex<T>(), 'N', mat);
+    return ev;
+   }
+
+   /// The eigensystems
+   template <typename M> std::pair<array<double, 1>, matrix<T>> eigenelements(M &mat) const {
+    _prepare(mat);
+    _invoke(is_complex<T>(), 'V', mat);
+    return {ev, _conj(mat, is_complex<T>())};
+   }
+
+   private:
+   mutable array<double, 1> ev, work2; // work2 only used for T complex
+   mutable array<T, 1> work;
+   mutable int dim, lwork, info;
+
+   // dispatch the implementation of invoke for T = double or complex
+   void _invoke(std::false_type, char compz, matrix_view<double> mat) const { // the case double (is_complex = false)
+    char uplo = 'U';
+    TRIQS_FORTRAN_MANGLING(dsyev)(&compz, &uplo, dim, mat.data_start(), dim, ev.data_start(), work.data_start(), lwork, info);
+    if (info) TRIQS_RUNTIME_ERROR << "eigenelements_worker :error code dsyev : " << info << " for matrix " << mat;
+   }
+
+   void _invoke(std::true_type, char compz, matrix_view<std::complex<double>> mat) const { // the case complex (is_complex = true)
+    char uplo = 'U';
+    TRIQS_FORTRAN_MANGLING(zheev)(&compz, &uplo, dim, mat.data_start(), dim, ev.data_start(), work.data_start(), lwork,
+                                  work2.data_start(), info);
+    if (info) TRIQS_RUNTIME_ERROR << "eigenelements_worker :error code zheev : " << info << " for matrix " << mat;
+   }
+
+   template <typename M> void _prepare(M const &mat) const {
+    if (mat.is_empty()) TRIQS_RUNTIME_ERROR << "eigenelements_worker : the matrix is empty : matrix =  " << mat << "  ";
+    if (!mat.is_square()) TRIQS_RUNTIME_ERROR << "eigenelements_worker : the matrix " << mat << " is not square ";
+    if (!mat.indexmap().is_contiguous())
+     TRIQS_RUNTIME_ERROR << "eigenelements_worker : the matrix " << mat << " is not contiguous in memory";
+    dim = first_dim(mat);
+    ev.resize(dim);
+    lwork = 64 * dim;
+    work.resize(lwork);
+    if (is_complex<T>::value) work2.resize(lwork);
+   }
+
+   template <typename M> matrix<double> _conj(M const &m, std::false_type) const { return m; }
+
+   // impl : since we call fortran lapack, if the order is C (!), the matrix is transposed, or conjugated, so we obtain
+   // the conjugate of the eigenvectors... Fix #119.
+   // Do nothing if the order is fortran already...
+   template <typename M> matrix<std::complex<double>> _conj(M const &m, std::true_type) const {
+    if (m.memory_layout_is_c())
+     return conj(m);
+    else
+     return m.transpose(); // the matrix mat is understood as a fortran matrix. After the lapack, in memory, it contains the
+                           // correct answer.
+    // but since it is a fortran matrix, the C will see its transpose. We need to compensate this transpose (!).
+   }
   };
 
- //--------------------------------
- extern "C" { 
- void TRIQS_FORTRAN_MANGLING(dsyev)(char*,char*,        // JOBZ and UPLO
-   int &,              // Matrix Size
-   double[],            // matrix
-   int&,               // LDA of the matrix
-   double[],           // Eigenvalues array
-   double[],            // WORK
-   int&,               // LWORK
-   int &               // INFO
-   );
+  //--------------------------------
 
+  /**
+   * Simple diagonalization call, return all eigenelements.
+   * Handles both real and complex case.
+   * @param M : the matrix or view. MUST be contiguous. It is modified by the call.
+   *            If you wish not to modify it, call eigenelements(make_clone(A))
+   */
+  template <typename M>
+  std::pair<array<double, 1>, matrix<typename std14::remove_reference_t<M>::value_type>> eigenelements(M &&m) {
+   return eigenelements_worker<typename std14::remove_reference_t<M>::value_type>().eigenelements(m);
+  }
 
- void TRIQS_FORTRAN_MANGLING(zheev)(char*,char*,        // JOBZ and UPLO
-   int &,              // Matrix Size
-   std::complex<double> [],            // matrix
-   int&,               // LDA of the matrix
-   double[],           // Eigenvalues array
-   std::complex<double>[],            // WORK
-   int &,               // LWORK
-   double[],  // WORK2
-   int &               // INFO
-   );
+  //--------------------------------
+
+  /**
+   * Simple diagonalization call, returning only the eigenvalues.
+   * Handles both real and complex case.
+   * @param M : the matrix VIEW : it MUST be contiguous
+   * @param take_copy : makes a copy of the matrix before calling lapack, so that the original is preserved.
+   *   if false : no copy is made and the content of the matrix M is destroyed.
+   *   if true : a copy is made, M is preserved, but of course it is slower...
+   */
+  template <typename M> array<double, 1> eigenvalues(M &&m) {
+   return eigenelements_worker<typename std14::remove_reference_t<M>::value_type>().eigenvalues(m);
+  }
  }
-
- //--------------------------------
-
- template<ull_t Opt, bool Compute_Eigenvectors >
-  struct eigenelements_worker< matrix_view<double,Opt> ,Compute_Eigenvectors > :eigenelements_worker_base<double,Opt,Compute_Eigenvectors> { 
-   eigenelements_worker ( matrix_view <double,Opt> the_matrix) : eigenelements_worker_base<double,Opt,Compute_Eigenvectors> (the_matrix) {}
-   void invoke() {
-    int info;
-    //fortran_int_t info;
-    TRIQS_FORTRAN_MANGLING(dsyev) (&this->compz,&this->uplo,this->dim,this->mat.data_start(),this->dim,this->ev.data_start(),this->work.data_start(),this->lwork,info);
-    if (info)  TRIQS_RUNTIME_ERROR<<"eigenelements_worker :error code dsyev : "<<info<<" for matrix "<<this->mat;
-    this->has_run = true;
-   }
-  };
-
- //--------------------------------
-
- template<ull_t Opt, bool Compute_Eigenvectors >
-  struct eigenelements_worker< matrix_view<std::complex<double>, Opt>,Compute_Eigenvectors > :eigenelements_worker_base<std::complex<double>,Opt,Compute_Eigenvectors> { 
-   triqs::arrays::vector <double> work2;   
-   public :
-   eigenelements_worker ( matrix_view <std::complex<double>,Opt> the_matrix) : eigenelements_worker_base<std::complex<double>,Opt,Compute_Eigenvectors> (the_matrix) {  work2.resize(this->lwork);}
-   void invoke() {
-    int info;
-    TRIQS_FORTRAN_MANGLING(zheev) (&this->compz,&this->uplo,this->dim,this->mat.data_start(),
-      this->dim,this->ev.data_start(),this->work.data_start(),this->lwork,this->work2.data_start(),info);
-    if (info)  TRIQS_RUNTIME_ERROR<<"eigenelements_worker :error code zheev : "<<info<<" for matrix "<<this->mat;
-    this->has_run = true;
-   }
-  };
-
- //--------------------------------
-
- /**
-  * Simple diagonalization call, return all eigenelements.
-  * Handles both real and complex case.
-  * @param M : the matrix VIEW : it MUST be contiguous
-  * @param take_copy : makes a copy of the matrix before calling lapack, so that the original is preserved.
-  *   if false : no copy is made and the content of the matrix M is destroyed (it is equal to vectors()).
-  *   if true : a copy is made, M is preserved, but of course it is slower...
-  */
- template<typename MatrixViewType >   
-  std::pair<array<double,1>, typename MatrixViewType::regular_type> eigenelements( MatrixViewType const & M, bool take_copy =false) { 
-   eigenelements_worker<typename MatrixViewType::view_type, true> W(take_copy ? MatrixViewType(make_clone(M)) : M); 
-   W.invoke(); 
-   return std::make_pair(W.values(),W.vectors());
-  }
-
- //--------------------------------
-
- /**
-  * Simple diagonalization call, returning only the eigenvalues.
-  * Handles both real and complex case.
-  * @param M : the matrix VIEW : it MUST be contiguous
-  * @param take_copy : makes a copy of the matrix before calling lapack, so that the original is preserved.
-  *   if false : no copy is made and the content of the matrix M is destroyed.
-  *   if true : a copy is made, M is preserved, but of course it is slower...
-  */
- template<typename MatrixViewType >   
-  triqs::arrays::vector_view <double> eigenvalues( MatrixViewType const & M, bool take_copy = false) { 
-   eigenelements_worker<MatrixViewType,false> W(take_copy ? MatrixViewType(make_clone(M)) : M); W.invoke(); return W.values(); 
-  }
-
-}}} // namespace triqs::arrays::linalg
-#endif
+}
+} // namespace triqs::arrays::linalg
 
diff --git a/triqs/lattice/tight_binding.cpp b/triqs/lattice/tight_binding.cpp
index 5187529d..9f38042e 100644
--- a/triqs/lattice/tight_binding.cpp
+++ b/triqs/lattice/tight_binding.cpp
@@ -59,7 +59,7 @@ namespace lattice {
   array<double, 2> eval(norb, n_pts);
   k_t dk = (K2 - K1) / double(n_pts), k = K1;
   for (int i = 0; i < n_pts; ++i, k += dk) {
-   eval(range(), i) = linalg::eigenvalues(TK(k(range(0, ndim)))(), false);
+   eval(range(), i) = linalg::eigenvalues(TK(k(range(0, ndim)))());
   }
   return eval;
  }
@@ -86,7 +86,7 @@ namespace lattice {
   grid_generator grid(ndim, n_pts);
   array<double, 2> eval(norb, grid.size());
   for (; grid; ++grid) {
-   eval(range(), grid.index()) = linalg::eigenvalues(TK((*grid)(range(0, ndim)))(), false);
+   eval(range(), grid.index()) = linalg::eigenvalues(TK((*grid)(range(0, ndim)))());
   }
   return eval;
  }