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

@@ -19,33 +19,53 @@
  *
  ******************************************************************************/
 #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);
 
   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);}
-
+   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;
- 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;
+  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);
+ }
+
+ {
+  matrix<double> A(3, 3);
   A() = 0;
   A(0, 1) = 1;
   A(1, 0) = 1;
@@ -53,23 +73,55 @@ int main(int argc, char **argv) {
   A(0, 2) = 2;
   A(2, 0) = 2;
 
+  auto B = A;
   std::cout << "A = " << A << std::endl;
- std::cout<<" values = "<< eigenelements(A(),true).first<<std::endl;
- std::cout<<" vectors = "<< eigenelements(A(),true).second<<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);
+ }
+
+ {
+  matrix<double> A(3, 3);
 
   A() = 0;
   A(0, 1) = 1;
   A(1, 0) = 1;
   A(2, 2) = 8;
 
+  auto B = A;
   std::cout << "A = " << A << std::endl;
- std::cout<<" vectors = "<< eigenelements(A(),true).second<<std::endl;
- std::cout<<" values = "<< eigenelements(A(),true).first<<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;
-
- return 0;
-
+  check_eig (A, w.second, w.first);
  }
 
+ { // 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;
+}
 

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;
 }

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,76 +18,17 @@
  * 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;
-
- 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;
-    }
-
-   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'; }
-   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;
-    }
-  };
-
- //--------------------------------
   extern "C" {
   void TRIQS_FORTRAN_MANGLING(dsyev)(char *, char *, // JOBZ and UPLO
                                      int &,          // Matrix Size
@@ -112,33 +53,70 @@ namespace triqs { namespace arrays { namespace linalg {
                                      );
   }
 
- //--------------------------------
-
- 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;   
+  /**
+   * A worker to call lapack routine with the matrices. Handles both real and complex case.
+   */
+  template <typename T> class eigenelements_worker {
    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;
+   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 (!).
    }
   };
 
@@ -147,16 +125,12 @@ namespace triqs { namespace arrays { namespace linalg {
   /**
    * 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...
+   * @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 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());
+  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);
   }
 
   //--------------------------------
@@ -169,11 +143,10 @@ namespace triqs { namespace arrays { namespace linalg {
    *   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(); 
+  template <typename M> array<double, 1> eigenvalues(M &&m) {
+   return eigenelements_worker<typename std14::remove_reference_t<M>::value_type>().eigenvalues(m);
   }
-
-}}} // namespace triqs::arrays::linalg
-#endif
+ }
+}
+} // namespace triqs::arrays::linalg
 

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;
  }