samerafach/Polymath
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' of https://git.afach.de/samerafach/Polymath

Browse Source
This commit is contained in:
Samer Afach 2017-09-02 23:19:31 +02:00
commit 67e7bb707b
2 changed files with 2181 additions and 2066 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
include/internal/LapackAdapters.h
View File

@ -12,26 +12,33 @@ typedef int lapack_int;
typedef std::complex<float> lapack_complex_float;
typedef std::complex<double> lapack_complex_double;
extern "C" void sgetrf_( lapack_int* m, lapack_int* n, float* a, lapack_int* lda, lapack_int* ipiv, lapack_int *info );
extern "C" void dgetrf_( lapack_int* m, lapack_int* n, double* a, lapack_int* lda, lapack_int* ipiv, lapack_int *info );
extern "C" void cgetrf_( lapack_int* m, lapack_int* n, lapack_complex_float* a, lapack_int* lda, lapack_int* ipiv, lapack_int *info );
extern "C" void zgetrf_( lapack_int* m, lapack_int* n, lapack_complex_double* a, lapack_int* lda, lapack_int* ipiv, lapack_int *info );
extern "C" void sgetri_( lapack_int* n, float* a, lapack_int* lda, lapack_int* ipiv, float* work, lapack_int* lwork, lapack_int *info );
extern "C" void dgetri_( lapack_int* n, double* a, lapack_int* lda, lapack_int* ipiv, double* work, lapack_int* lwork, lapack_int *info );
extern "C" void cgetri_( lapack_int* n, lapack_complex_float* a, lapack_int* lda, lapack_int* ipiv, lapack_complex_float* work, lapack_int* lwork, lapack_int *info );
extern "C" void zgetri_( lapack_int* n, lapack_complex_double* a, lapack_int* lda, lapack_int* ipiv, lapack_complex_double* work, lapack_int* lwork, lapack_int *info );
extern "C" void sspsv_( char* uplo, lapack_int* n, lapack_int* nrhs, float* ap, lapack_int* ipiv, float* b, lapack_int* ldb, lapack_int *info );
extern "C" void dspsv_( char* uplo, lapack_int* n, lapack_int* nrhs, double* ap, lapack_int* ipiv, double* b, lapack_int* ldb, lapack_int *info );
extern "C" void cspsv_( char* uplo, lapack_int* n, lapack_int* nrhs, lapack_complex_float* ap, lapack_int* ipiv, lapack_complex_float* b, lapack_int* ldb, lapack_int *info );
extern "C" void zspsv_( char* uplo, lapack_int* n, lapack_int* nrhs, lapack_complex_double* ap, lapack_int* ipiv, lapack_complex_double* b, lapack_int* ldb, lapack_int *info );
extern "C" void cheev_( char* jobz, char* uplo, lapack_int* n, lapack_complex_float* a, lapack_int* lda, float* w, lapack_complex_float* work, lapack_int* lwork, float* rwork, lapack_int *info );
extern "C" void zheev_( char* jobz, char* uplo, lapack_int* n, lapack_complex_double* a, lapack_int* lda, double* w, lapack_complex_double* work, lapack_int* lwork, double* rwork, lapack_int *info );
extern "C" void cheevd_( char* jobz, char* uplo, lapack_int* n, lapack_complex_float* a, lapack_int* lda, float* w, lapack_complex_float* work, lapack_int* lwork, float* rwork, lapack_int* lrwork, lapack_int* iwork, lapack_int* liwork, lapack_int *info );
extern "C" void zheevd_( char* jobz, char* uplo, lapack_int* n, lapack_complex_double* a, lapack_int* lda, double* w, lapack_complex_double* work, lapack_int* lwork, double* rwork, lapack_int* lrwork, lapack_int* iwork, lapack_int* liwork, lapack_int *info );
extern "C" int sgetrf_( lapack_int* m, lapack_int* n, float* a, lapack_int* lda, lapack_int* ipiv, lapack_int *info );
extern "C" int dgetrf_( lapack_int* m, lapack_int* n, double* a, lapack_int* lda, lapack_int* ipiv, lapack_int *info );
extern "C" int cgetrf_( lapack_int* m, lapack_int* n, lapack_complex_float* a, lapack_int* lda, lapack_int* ipiv, lapack_int *info );
extern "C" int zgetrf_( lapack_int* m, lapack_int* n, lapack_complex_double* a, lapack_int* lda, lapack_int* ipiv, lapack_int *info );
extern "C" int sgetri_( lapack_int* n, float* a, lapack_int* lda, lapack_int* ipiv, float* work, lapack_int* lwork, lapack_int *info );
extern "C" int dgetri_( lapack_int* n, double* a, lapack_int* lda, lapack_int* ipiv, double* work, lapack_int* lwork, lapack_int *info );
extern "C" int cgetri_( lapack_int* n, lapack_complex_float* a, lapack_int* lda, lapack_int* ipiv, lapack_complex_float* work, lapack_int* lwork, lapack_int *info );
extern "C" int zgetri_( lapack_int* n, lapack_complex_double* a, lapack_int* lda, lapack_int* ipiv, lapack_complex_double* work, lapack_int* lwork, lapack_int *info );
extern "C" int sspsv_( char* uplo, lapack_int* n, lapack_int* nrhs, float* ap, lapack_int* ipiv, float* b, lapack_int* ldb, lapack_int *info );
extern "C" int dspsv_( char* uplo, lapack_int* n, lapack_int* nrhs, double* ap, lapack_int* ipiv, double* b, lapack_int* ldb, lapack_int *info );
extern "C" int cspsv_( char* uplo, lapack_int* n, lapack_int* nrhs, lapack_complex_float* ap, lapack_int* ipiv, lapack_complex_float* b, lapack_int* ldb, lapack_int *info );
extern "C" int zspsv_( char* uplo, lapack_int* n, lapack_int* nrhs, lapack_complex_double* ap, lapack_int* ipiv, lapack_complex_double* b, lapack_int* ldb, lapack_int *info );
extern "C" int cheev_( char* jobz, char* uplo, lapack_int* n, lapack_complex_float* a, lapack_int* lda, float* w, lapack_complex_float* work, lapack_int* lwork, float* rwork, lapack_int *info );
extern "C" int zheev_( char* jobz, char* uplo, lapack_int* n, lapack_complex_double* a, lapack_int* lda, double* w, lapack_complex_double* work, lapack_int* lwork, double* rwork, lapack_int *info );
extern "C" int cheevd_( char* jobz, char* uplo, lapack_int* n, lapack_complex_float* a, lapack_int* lda, float* w, lapack_complex_float* work, lapack_int* lwork, float* rwork, lapack_int* lrwork, lapack_int* iwork, lapack_int* liwork, lapack_int *info );
extern "C" int zheevd_( char* jobz, char* uplo, lapack_int* n, lapack_complex_double* a, lapack_int* lda, double* w, lapack_complex_double* work, lapack_int* lwork, double* rwork, lapack_int* lrwork, lapack_int* iwork, lapack_int* liwork, lapack_int *info );
extern "C" void sgemm_(char * transa, char * transb, lapack_int * m, lapack_int * n, lapack_int * k, float* alpha, float * A, lapack_int * lda, float * B, lapack_int * ldb, float * beta, float *, lapack_int * ldc);
extern "C" void dgemm_(char * transa, char * transb, lapack_int * m, lapack_int * n, lapack_int * k, double * alpha, double * A, lapack_int * lda, double * B, lapack_int * ldb, double * beta, double *, lapack_int * ldc);
extern "C" void cgemm_(char*,char*,lapack_int*,lapack_int*,lapack_int*, lapack_complex_float*, lapack_complex_float*, lapack_int*, lapack_complex_float*,lapack_int*, lapack_complex_float*, lapack_complex_float*,lapack_int*);
extern "C" void zgemm_(char*, char*, lapack_int*, lapack_int*, lapack_int*, lapack_complex_double*, lapack_complex_double*,lapack_int*, lapack_complex_double*,lapack_int*, lapack_complex_double*,lapack_complex_double*,lapack_int*);
extern "C" void cgemm_(char* , char* , lapack_int* ,lapack_int* , lapack_int* , lapack_complex_float*, lapack_complex_float*, lapack_int*, lapack_complex_float*,lapack_int*, lapack_complex_float*, lapack_complex_float*,lapack_int*);
extern "C" void zgemm_(char* , char* , lapack_int* , lapack_int* , lapack_int* , lapack_complex_double*, lapack_complex_double*,lapack_int*, lapack_complex_double*,lapack_int*, lapack_complex_double*,lapack_complex_double*,lapack_int*);
extern "C" int ssymm_(char *side, char *uplo, lapack_int *m, lapack_int *n, float *alpha, float *a, lapack_int *lda, float *b, lapack_int *ldb, float *beta, float *c__, lapack_int *ldc);
extern "C" int dsymm_(char *side, char *uplo, lapack_int *m, lapack_int *n, double *alpha, double *a, lapack_int *lda, double *b, lapack_int *ldb, double *beta, double *c__, lapack_int *ldc);
extern "C" int csymm_(char *side, char *uplo, lapack_int *m, lapack_int *n, lapack_complex_float *alpha, lapack_complex_float *a, lapack_int *lda, lapack_complex_float *b, lapack_int *ldb, lapack_complex_float *beta, lapack_complex_float *c__, lapack_int *ldc);
extern "C" int zsymm_(char *side, char *uplo, lapack_int *m, lapack_int *n, lapack_complex_double *alpha, lapack_complex_double *a, lapack_int *lda, lapack_complex_double *b, lapack_int *ldb, lapack_complex_double *beta, lapack_complex_double *c__, lapack_int *ldc);
extern "C" int chemm_(char *side, char *uplo, lapack_int *m, lapack_int *n, lapack_complex_float *alpha, lapack_complex_float *a, lapack_int *lda, lapack_complex_float *b, lapack_int *ldb, lapack_complex_float *beta, lapack_complex_float *c__, lapack_int *ldc);
extern "C" int zhemm_(char *side, char *uplo, lapack_int *m, lapack_int *n, lapack_complex_double *alpha, lapack_complex_double *a, lapack_int *lda, lapack_complex_double *b, lapack_int *ldb, lapack_complex_double *beta, lapack_complex_double *c__, lapack_int *ldc);
//#endif

154
include/internal/Matrix.h
View File

@ -16,6 +16,7 @@
#include <algorithm>
#include <chrono>
#include <iomanip>
#include <random>
#include "LapackAdapters.h"
#include "StdAdapters.h"
@ -106,7 +107,7 @@ to_string(const T& a_value, int precision)
{
std::ostringstream out;
out << std::setprecision(precision) << a_value.real();
out << "+";
out << (a_value.imag() < 0 ? "" : "+"); //"-" comes from the number itself
out << std::setprecision(precision) << a_value.imag();
out << Poly::ComplexUnit;
return out.str();
@ -133,9 +134,9 @@ const static EIGENS_METHOD METHOD_DIVIDE_CONQUER = 1;
template <typename T, int M>
void MultiplyMatrices(const Matrix<T,M> &matrixLHS, const Matrix<T,M> &matrixRHS, Matrix<T,M> &result);
void MultiplyMatrices(const Matrix<T,M> &matrixLHS, const Matrix<T,M> &matrixRHS, Matrix<T,M> &result, const MATRIX_TYPE& lhs_type = MATRIX_GENERAL, const MATRIX_TYPE& rhs_type = MATRIX_GENERAL);
template <typename T, int M>
void MultiplyMatrices(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs, Matrix<T,M> &to_add_then_result, T alpha, T beta);
void MultiplyMatrices(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs, Matrix<T,M> &to_add_then_result, T alpha, T beta, const MATRIX_TYPE& lhs_type = MATRIX_GENERAL, const MATRIX_TYPE& rhs_type = MATRIX_GENERAL);
template <typename T, int M>
Matrix<T,M> IdentityMatrix(const typename Matrix<T,M>::size_type& size);
template <typename T, typename U, int M>
@ -192,7 +193,7 @@ public:
typename std::enable_if<std::is_same<T, std::complex<U>>::value, void>::type
copyFrom(const Matrix<U, M> & rhs);
void initialize(const std::initializer_list<T>& list, Matrix<T,M>::size_type start_row = 0, Matrix<T,M>::size_type start_column = 0);
void initialize(const std::initializer_list<T>& list, typename Matrix<T,M>::size_type start_row = 0, typename Matrix<T,M>::size_type start_column = 0);
bool operator==(const Matrix<T,M>& rhs);
bool operator!=(const Matrix<T,M>& rhs);
template <typename _T, int _M>
@ -394,9 +395,9 @@ void Matrix<T,M>::resize(size_type rows, size_type columns, const T& initialValu
this->_matEls.resize(newSize);
this->_rows = rows;
this->_columns = columns;
for(size_type i = 0; i < oldRows; i++)
for(size_type i = 0; i < (oldRows > rows ? rows : oldRows); i++)
{
for(size_type j = 0; j < oldColumns; j++)
for(size_type j = 0; j < (oldColumns > columns ? columns : oldColumns); j++)
{
this->at(i,j) = oldMat.at(i,j);
}
@ -482,6 +483,7 @@ bool Matrix<T,M>::operator!=(const Matrix<T,M>& rhs)
return !(*this == rhs);
}
//FIXME: Consider non-square matrices
template <typename T, int M>
void Matrix<T,M>::transpose_inplace()
{
@ -557,7 +559,6 @@ Matrix<std::complex<T>,M> operator*(const T& lhs, const Matrix<std::complex<T>,M
return result;
}
template <typename T, int M>
Matrix<T,M> operator*(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs)
{
@ -689,28 +690,32 @@ Matrix<std::complex<T>,M> operator-(const Matrix<T,M>& lhs, const Matrix<std::co
template <typename T, int M>
Matrix<T,M>& Matrix<T,M>::operator+=(const Matrix<T,M>& rhs)
{
(*this) = (*this) + rhs;
std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::plus<T>());
// (*this) = (*this) + rhs;
return *this;
}
template <typename T, int M>
Matrix<T,M>& Matrix<T,M>::operator-=(const Matrix<T,M>& rhs)
{
(*this) = (*this) - rhs;
std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::minus<T>());
// (*this) = (*this) - rhs;
return *this;
}
template <typename T, int M>
Matrix<T,M>& Matrix<T,M>::operator*=(const Matrix<T,M>& rhs)
{
(*this) = (*this) * rhs;
std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::multiplies<T>());
// (*this) = (*this) * rhs;
return *this;
}
template <typename T, int M>
Matrix<T,M>& Matrix<T,M>::operator*=(const T& rhs)
{
(*this) = (*this) * rhs;
std::transform(this->begin(),this->end(),this->begin(),std::bind2nd(std::multiplies<T>(),rhs));
// (*this) = (*this) * rhs;
return *this;
}
@ -926,7 +931,7 @@ Matrix<T, M>::operator=(const Matrix<U, M> & rhs)
}
template <typename T, int M>
void Matrix<T,M>::initialize(const std::initializer_list<T> &list, Matrix<T,M>::size_type start_row, Matrix<T,M>::size_type start_column)
void Matrix<T,M>::initialize(const std::initializer_list<T> &list, typename Matrix<T,M>::size_type start_row, typename Matrix<T,M>::size_type start_column)
{
size_type start_from = RowColumnToElement(start_row,start_column,this->rows(),this->columns());
if(list.size() + start_from > this->_matEls.size()) //if tried to initalize with an array bigger than
@ -1493,6 +1498,93 @@ std::string _lapack_eigens_exception_converge_error(int info);
std::string _lapack_unsupported_type();
std::string _unsupported_type(const std::string& function_name);
template <typename T, int M>
void _Lapack_multiply_matrix_symmetric(const Matrix<T,M>& lhs_i, const Matrix<T,M>& rhs_i, bool sym_mat_matrix_is_left, Matrix<T,M> &to_add_then_result, T alpha, T beta)
{
char sym_mat_is_left_or_right = (sym_mat_matrix_is_left ? 'L' : 'R');
const Matrix<T,M> *lhs, *rhs;
//whether the first matrix is on the left or right is determined by the first parameter in the LAPACK call
if(sym_mat_matrix_is_left)
{
lhs = &lhs_i;
rhs = &rhs_i;
}
else
{
lhs = &rhs_i;
rhs = &lhs_i;
}
char uplo = 'U';
int m = static_cast<int>(lhs->rows());
int n = static_cast<int>(rhs->columns());
int k = static_cast<int>(lhs->columns());
int lda = std::max({1,k});
int ldb = std::max({1,n});
int ldc = std::max({1,m});
if(std::is_same<T,double>::value)
{
typedef double TP;
dsymm_(&sym_mat_is_left_or_right, &uplo, &m, &n, (TP*)&alpha, (TP*)&lhs->front(), &lda, (TP*)&rhs->front(), &ldb, (TP*)&beta, (TP*)&to_add_then_result.front(), &ldc);
}
else if(std::is_same<T,lapack_complex_double>::value)
{
typedef lapack_complex_double TP;
zsymm_(&sym_mat_is_left_or_right, &uplo, &m, &n, (TP*)&alpha, (TP*)&lhs->front(), &lda, (TP*)&rhs->front(), &ldb, (TP*)&beta, (TP*)&to_add_then_result.front(), &ldc);
}
else if(std::is_same<T,float>::value)
{
typedef float TP;
ssymm_(&sym_mat_is_left_or_right, &uplo, &m, &n, (TP*)&alpha, (TP*)&lhs->front(), &lda, (TP*)&rhs->front(), &ldb, (TP*)&beta, (TP*)&to_add_then_result.front(), &ldc);
}
else if(std::is_same<T,lapack_complex_float>::value)
{
typedef lapack_complex_float TP;
csymm_(&sym_mat_is_left_or_right, &uplo, &m, &n, (TP*)&alpha, (TP*)&lhs->front(), &lda, (TP*)&rhs->front(), &ldb, (TP*)&beta, (TP*)&to_add_then_result.front(), &ldc);
}
else
{
throw std::logic_error(_lapack_unsupported_type());
}
}
template <typename T, int M>
void _Lapack_multiply_matrix_hermitian(const Matrix<T,M>& lhs_i, const Matrix<T,M>& rhs_i, bool herm_mat_matrix_is_left, Matrix<T,M> &to_add_then_result, T alpha, T beta)
{
char herm_mat_is_left_or_right = (herm_mat_matrix_is_left ? 'L' : 'R');
const Matrix<T,M> *lhs, *rhs;
//whether the first matrix is on the left or right is determined by the first parameter in the LAPACK call
if(herm_mat_matrix_is_left)
{
lhs = &lhs_i;
rhs = &rhs_i;
}
else
{
lhs = &rhs_i;
rhs = &lhs_i;
}
char uplo = 'U';
int m = static_cast<int>(lhs->rows());
int n = static_cast<int>(rhs->columns());
int k = static_cast<int>(lhs->columns());
int lda = std::max({1,k});
int ldb = std::max({1,n});
int ldc = std::max({1,m});
if(std::is_same<T,lapack_complex_double>::value)
{
typedef lapack_complex_double TP;
zhemm_(&herm_mat_is_left_or_right, &uplo, &m, &n, (TP*)&alpha, (TP*)&lhs->front(), &lda, (TP*)&rhs->front(), &ldb, (TP*)&beta, (TP*)&to_add_then_result.front(), &ldc);
}
else if(std::is_same<T,lapack_complex_float>::value)
{
typedef lapack_complex_float TP;
chemm_(&herm_mat_is_left_or_right, &uplo, &m, &n, (TP*)&alpha, (TP*)&lhs->front(), &lda, (TP*)&rhs->front(), &ldb, (TP*)&beta, (TP*)&to_add_then_result.front(), &ldc);
}
else
{
throw std::logic_error(_lapack_unsupported_type());
}
}
template <typename T, int M>
void _Lapack_multiply_matrix_general(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs, Matrix<T,M> &to_add_then_result, T alpha, T beta)
@ -1536,7 +1628,13 @@ void _Lapack_multiply_matrix_general(const Matrix<T,M>& lhs, const Matrix<T,M>&
}
template <typename T, int M>
void MultiplyMatrices(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs, Matrix<T,M> &to_add_then_result, T alpha, T beta)
void MultiplyMatrices(const Matrix<T,M>& lhs,
const Matrix<T,M>& rhs,
Matrix<T,M> &to_add_then_result,
T alpha,
T beta,
const MATRIX_TYPE& lhs_type,
const MATRIX_TYPE& rhs_type)
{
// std::cout<<"multiplication: "<<lhs.columns()<<"\t"<<rhs.rows()<<std::endl;
#ifdef POLYMATH_DEBUG
@ -1550,9 +1648,21 @@ void MultiplyMatrices(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs, Matrix<T,M
std::is_same<T,float>::value ||
std::is_same<T,lapack_complex_double>::value ||
std::is_same<T,lapack_complex_float>::value)
{
if(lhs_type == MATRIX_SYMMETRIC || rhs_type == MATRIX_SYMMETRIC)
{
_Lapack_multiply_matrix_symmetric(lhs,rhs,(lhs_type == MATRIX_SYMMETRIC ? true : false),to_add_then_result,alpha,beta);
}
else if((lhs_type == MATRIX_HERMITIAN || rhs_type == MATRIX_HERMITIAN) &&
(std::is_same<T,lapack_complex_float>::value || std::is_same<T,lapack_complex_double>::value))
{
_Lapack_multiply_matrix_hermitian(lhs,rhs,(lhs_type == MATRIX_HERMITIAN ? true : false),to_add_then_result,alpha,beta);
}
else
{
_Lapack_multiply_matrix_general(lhs,rhs,to_add_then_result,alpha,beta);
}
}
else
{
T comp;
@ -1574,9 +1684,11 @@ void MultiplyMatrices(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs, Matrix<T,M
template <typename T, int M>
void MultiplyMatrices(const Matrix<T,M> &matrixLHS,
const Matrix<T,M> &matrixRHS,
Matrix<T,M> &result)
Matrix<T,M> &result,
const MATRIX_TYPE& lhs_type,
const MATRIX_TYPE& rhs_type)
{
Poly::MultiplyMatrices(matrixLHS,matrixRHS,result,T(1),T(0));
Poly::MultiplyMatrices(matrixLHS,matrixRHS,result,T(1),T(0),lhs_type,rhs_type);
}
template <typename T, int M = ColMaj>
@ -1589,24 +1701,20 @@ RandomMatrix(const typename Matrix<T,M>::size_type& rows,
const MATRIX_TYPE& type = Poly::MATRIX_GENERAL)
{
Matrix<T,M> result(rows,columns);
if(std::is_integral<T>::value)
{
if(std::is_integral<T>::value) {
std::uniform_int_distribution<T> distribution(min,max);
std::default_random_engine generator(seed);
std::generate(result.begin(), result.end(), [&]() { return distribution(generator); });
}
if(type == MATRIX_SYMMETRIC || type == MATRIX_HERMITIAN)
{
if(type == MATRIX_SYMMETRIC || type == MATRIX_HERMITIAN) {
result = result/T(2) + Poly::Transpose<T,M>(result/T(2));
}
else if(type == MATRIX_ANTISYMMETRIC || type == MATRIX_ANTIHERMITIAN)
{
else if(type == MATRIX_ANTISYMMETRIC || type == MATRIX_ANTIHERMITIAN) {
result = result/T(2) - Poly::Transpose<T,M>(result/T(2));
}
else if(type == MATRIX_GENERAL) {}
else
{
else {
throw std::logic_error(_unsupported_type(__func__).c_str());
}
return result;