inline void Math::matrixEig(__attribute__((unused)) UInt n,
__attribute__((unused)) T * A,
__attribute__((unused)) T * d,
__attribute__((unused)) T * V) {
AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function! Or implement it for the type " << debug::demangle(typeid(T).name()));
}
-#ifdef AKANTU_USE_LAPACK
template<>
inline void Math::matrixEig<double>(UInt n, double * A, double * d, double * V) {
// Matrix A is row major, so the lapack function in fortran will process
// A^t. Asking for the left eigenvectors of A^t will give the transposed right
// eigenvectors of A so in the C++ code the right eigenvectors.
char jobvl;
if(V != NULL)
jobvl = 'V'; // compute left eigenvectors
else
jobvl = 'N'; // compute left eigenvectors
char jobvr('N'); // compute right eigenvectors
double * di = new double[n]; // imaginary part of the eigenvalues
inline void Math::inv(__attribute__((unused)) UInt n,
__attribute__((unused)) const T * A,
__attribute__((unused)) T * Ainv) {
AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function! Or implement it for the type " << debug::demangle(typeid(T).name()));
}
#ifdef AKANTU_USE_LAPACK
template<>
inline void Math::inv<double>(UInt n, const double * A, double * invA) {
int N = n;
int info;
int * ipiv = new int[N+1];
int lwork = N*N;
double * work = new double[lwork];
std::copy(A, A + n*n, invA);
dgetrf_(&N, &N, invA, &N, ipiv, &info);
if(info > 0) {
AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: "
<< info <<" )");
}
dgetri_(&N, invA, &N, ipiv, work, &lwork, &info);
if(info != 0) {
AKANTU_DEBUG_ERROR("Cannot invert the matrix (info: "<< info <<" )");