Commit ed1c7654 authored by Alvarez, Gonzalo's avatar Alvarez, Gonzalo
Browse files

inverse now in Matrix.h

parent 64f652d7

src/Matrix.cpp

+0 −52
Original line number Diff line number Diff line
@@ -210,57 +210,5 @@ void geev(char jobvl, 
#endif

}



template<typename T>

typename std::enable_if<IsComplexNumber<T>::True, void>::type

inverse(Matrix<T> &m)

{

#ifdef NO_LAPACK

	throw RuntimeError("inverse: NO LAPACK!\n");

#else

	int n = m.rows();

	int info = 0;

	Vector<int>::Type ipiv(n,0);

	psimag::LAPACK::GETRF(n, n, &(m(0,0)), n, &(ipiv[0]), info);

	int lwork = -1;

	typename Vector<T>::Type work(2);

	psimag::LAPACK::GETRI(n, &(m(0,0)), n, &(ipiv[0]), &(work[0]), lwork, info);

	lwork = static_cast<int>(PsimagLite::real(work[0]));

	work.resize(lwork+2);

	psimag::LAPACK::GETRI(n, &(m(0,0)), n, &(ipiv[0]), &(work[0]), lwork, info);

	String s = "[cz]getri_ failed\n";

	if (info!=0) throw RuntimeError(s.c_str());

#endif

}



template<typename T>

typename std::enable_if<Loki::TypeTraits<T>::isArith, void>::type

inverse(Matrix<T> &m)

{

#ifdef NO_LAPACK

	throw RuntimeError("inverse: NO LAPACK!\n");

#else

	int n = m.rows();

	int info = 0;

	Vector<int>::Type ipiv(n,0);

	psimag::LAPACK::GETRF( n, n, &(m(0,0)), n, &(ipiv[0]), info);

	int lwork = -1;

	typename Vector<T>::Type work(2);

	psimag::LAPACK::GETRI(n, &(m(0,0)), n, &(ipiv[0]), &(work[0]), lwork, info);

	lwork = static_cast<int>(work[0]);

	work.resize(lwork+2);

	psimag::LAPACK::GETRI(n, &(m(0,0)), n, &(ipiv[0]), &(work[0]), lwork, info);

	String s = "[sd]getri_ failed\n";

	if (info!=0) throw RuntimeError(s.c_str());

#endif

}



template void inverse<double>(Matrix<double>&);



template void inverse<float>(Matrix<float>&);



template void inverse<std::complex<double> >(Matrix<std::complex<double> >&);



template void inverse<std::complex<float> >(Matrix<std::complex<float> >&);



} // namespace PsimagLite

src/Matrix.h

+42 −5
Original line number Diff line number Diff line
@@ -528,13 +528,50 @@ void diag(Matrix<float> &m,Vector<float> ::Type& eigs,char option); 


void diag(Matrix<std::complex<float> > &m,Vector<float> ::Type& eigs,char option);



void inverse(Matrix<std::complex<double> >& m);



void inverse(Matrix<double>& m);

template<typename T>

typename std::enable_if<IsComplexNumber<T>::True, void>::type

inverse(Matrix<T> &m)

{

#ifdef NO_LAPACK

	throw RuntimeError("inverse: NO LAPACK!\n");

#else

	int n = m.rows();

	int info = 0;

	Vector<int>::Type ipiv(n,0);

	psimag::LAPACK::GETRF(n, n, &(m(0,0)), n, &(ipiv[0]), info);

	int lwork = -1;

	typename Vector<T>::Type work(2);

	psimag::LAPACK::GETRI(n, &(m(0,0)), n, &(ipiv[0]), &(work[0]), lwork, info);

	lwork = static_cast<int>(PsimagLite::real(work[0]));

	work.resize(lwork+2);

	psimag::LAPACK::GETRI(n, &(m(0,0)), n, &(ipiv[0]), &(work[0]), lwork, info);

	String s = "[cz]getri_ failed\n";

	if (info!=0) throw RuntimeError(s.c_str());

#endif

}



void inverse(Matrix<std::complex<float> >& m);

template<typename T>

typename std::enable_if<Loki::TypeTraits<T>::isArith, void>::type

inverse(Matrix<T> &m)

{

#ifdef NO_LAPACK

	throw RuntimeError("inverse: NO LAPACK!\n");

#else

	int n = m.rows();

	int info = 0;

	Vector<int>::Type ipiv(n,0);

	psimag::LAPACK::GETRF( n, n, &(m(0,0)), n, &(ipiv[0]), info);

	int lwork = -1;

	typename Vector<T>::Type work(2);

	psimag::LAPACK::GETRI(n, &(m(0,0)), n, &(ipiv[0]), &(work[0]), lwork, info);

	lwork = static_cast<int>(work[0]);

	work.resize(lwork+2);

	psimag::LAPACK::GETRI(n, &(m(0,0)), n, &(ipiv[0]), &(work[0]), lwork, info);

	String s = "[sd]getri_ failed\n";

	if (info!=0) throw RuntimeError(s.c_str());

#endif

}



void inverse(Matrix<float>& m);

// end in Matrix.cpp



template<typename T>