mirror of
https://github.com/triqs/dft_tools
synced 2024-11-01 19:53:45 +01:00
072b45ac1c
Not in the concept, not needed, just an annoyance. replaced by free functions : first_dim(A), second_dim(A), get_shape(A) and so on...
180 lines
7.4 KiB
C++
180 lines
7.4 KiB
C++
/*******************************************************************************
|
|
*
|
|
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
|
*
|
|
* Copyright (C) 2011 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
|
|
* Foundation, either version 3 of the License, or (at your option) any later
|
|
* version.
|
|
*
|
|
* TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
|
|
* details.
|
|
*
|
|
* You should have received a copy of the GNU General Public License along with
|
|
* TRIQS. If not, see <http://www.gnu.org/licenses/>.
|
|
*
|
|
******************************************************************************/
|
|
#ifndef TRIQS_ARRAY_EIGENELEMENTS_H
|
|
#define TRIQS_ARRAY_EIGENELEMENTS_H
|
|
#include <type_traits>
|
|
#include "../array.hpp"
|
|
#include "../matrix.hpp"
|
|
#include "../vector.hpp"
|
|
#include <triqs/utility/exceptions.hpp>
|
|
|
|
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
|
|
double[], // matrix
|
|
int&, // LDA of the matrix
|
|
double[], // Eigenvalues array
|
|
double[], // WORK
|
|
int&, // LWORK
|
|
int & // INFO
|
|
);
|
|
|
|
|
|
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
|
|
);
|
|
}
|
|
|
|
//--------------------------------
|
|
|
|
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
|
|
|