From 07043ecb2d552a26280124675ca4f8463c56d8b4 Mon Sep 17 00:00:00 2001 From: Samer Afach Date: Tue, 8 Nov 2016 16:03:47 +0100 Subject: [PATCH 1/4] Inplace operators use std::transform(); and fixed resize problem when resizing from bigger to smaller matrix. --- include/internal/Matrix.h | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/include/internal/Matrix.h b/include/internal/Matrix.h index d26c19e..d34a91e 100644 --- a/include/internal/Matrix.h +++ b/include/internal/Matrix.h @@ -394,9 +394,9 @@ void Matrix::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 +482,7 @@ bool Matrix::operator!=(const Matrix& rhs) return !(*this == rhs); } +//FIXME: Consider non-square matrices template void Matrix::transpose_inplace() { @@ -689,28 +690,32 @@ Matrix,M> operator-(const Matrix& lhs, const Matrix Matrix& Matrix::operator+=(const Matrix& rhs) { - (*this) = (*this) + rhs; + std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::plus()); +// (*this) = (*this) + rhs; return *this; } template Matrix& Matrix::operator-=(const Matrix& rhs) { - (*this) = (*this) - rhs; + std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::minus()); +// (*this) = (*this) - rhs; return *this; } template Matrix& Matrix::operator*=(const Matrix& rhs) { - (*this) = (*this) * rhs; + std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::multiplies()); +// (*this) = (*this) * rhs; return *this; } template Matrix& Matrix::operator*=(const T& rhs) { - (*this) = (*this) * rhs; + std::transform(this->begin(),this->end(),this->begin(),std::bind2nd(std::multiplies(),rhs)); +// (*this) = (*this) * rhs; return *this; } From 7e6076af68157d9a1162385123e74274065f1927 Mon Sep 17 00:00:00 2001 From: Samer Afach Date: Thu, 10 Nov 2016 20:12:43 +0100 Subject: [PATCH 2/4] Optimizations and fixes for Windows compatibility. --- include/internal/LapackAdapters.h | 32 +++++++++++++++---------------- include/internal/Matrix.h | 4 ++-- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/include/internal/LapackAdapters.h b/include/internal/LapackAdapters.h index 2b646e3..3ebcdf2 100644 --- a/include/internal/LapackAdapters.h +++ b/include/internal/LapackAdapters.h @@ -12,22 +12,22 @@ typedef int lapack_int; typedef std::complex lapack_complex_float; typedef std::complex 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*); diff --git a/include/internal/Matrix.h b/include/internal/Matrix.h index d34a91e..bbcbb8d 100644 --- a/include/internal/Matrix.h +++ b/include/internal/Matrix.h @@ -192,7 +192,7 @@ public: typename std::enable_if>::value, void>::type copyFrom(const Matrix & rhs); - void initialize(const std::initializer_list& list, Matrix::size_type start_row = 0, Matrix::size_type start_column = 0); + void initialize(const std::initializer_list& list, typename Matrix::size_type start_row = 0, typename Matrix::size_type start_column = 0); bool operator==(const Matrix& rhs); bool operator!=(const Matrix& rhs); template @@ -931,7 +931,7 @@ Matrix::operator=(const Matrix & rhs) } template -void Matrix::initialize(const std::initializer_list &list, Matrix::size_type start_row, Matrix::size_type start_column) +void Matrix::initialize(const std::initializer_list &list, typename Matrix::size_type start_row, typename Matrix::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 From d09bf4d2bfd4036b3d52430afb6076f809d880ca Mon Sep 17 00:00:00 2001 From: Samer Afach Date: Fri, 11 Nov 2016 19:18:59 +0100 Subject: [PATCH 3/4] -Added Hermitian and symmetric matrices multiplication -Fixed imaginary numbers printing by removing excess "+" when imaginary part is negative - --- include/internal/LapackAdapters.h | 11 ++- include/internal/Matrix.h | 134 ++++++++++++++++++++++++++---- 2 files changed, 127 insertions(+), 18 deletions(-) diff --git a/include/internal/LapackAdapters.h b/include/internal/LapackAdapters.h index 3ebcdf2..6687219 100644 --- a/include/internal/LapackAdapters.h +++ b/include/internal/LapackAdapters.h @@ -30,8 +30,15 @@ extern "C" int cheevd_( char* jobz, char* uplo, lapack_int* n, lapack_complex_f 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 diff --git a/include/internal/Matrix.h b/include/internal/Matrix.h index bbcbb8d..20282ef 100644 --- a/include/internal/Matrix.h +++ b/include/internal/Matrix.h @@ -106,7 +106,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 +133,9 @@ const static EIGENS_METHOD METHOD_DIVIDE_CONQUER = 1; template -void MultiplyMatrices(const Matrix &matrixLHS, const Matrix &matrixRHS, Matrix &result); +void MultiplyMatrices(const Matrix &matrixLHS, const Matrix &matrixRHS, Matrix &result, const MATRIX_TYPE& lhs_type = MATRIX_GENERAL, const MATRIX_TYPE& rhs_type = MATRIX_GENERAL); template -void MultiplyMatrices(const Matrix& lhs, const Matrix& rhs, Matrix &to_add_then_result, T alpha, T beta); +void MultiplyMatrices(const Matrix& lhs, const Matrix& rhs, Matrix &to_add_then_result, T alpha, T beta, const MATRIX_TYPE& lhs_type = MATRIX_GENERAL, const MATRIX_TYPE& rhs_type = MATRIX_GENERAL); template Matrix IdentityMatrix(const typename Matrix::size_type& size); template @@ -558,7 +558,6 @@ Matrix,M> operator*(const T& lhs, const Matrix,M return result; } - template Matrix operator*(const Matrix& lhs, const Matrix& rhs) { @@ -1498,6 +1497,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 +void _Lapack_multiply_matrix_symmetric(const Matrix& lhs_i, const Matrix& rhs_i, bool sym_mat_matrix_is_left, Matrix &to_add_then_result, T alpha, T beta) +{ + char sym_mat_is_left_or_right = (sym_mat_matrix_is_left ? 'L' : 'R'); + const Matrix *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(lhs->rows()); + int n = static_cast(rhs->columns()); + int k = static_cast(lhs->columns()); + int lda = std::max({1,k}); + int ldb = std::max({1,n}); + int ldc = std::max({1,m}); + if(std::is_same::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::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::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::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 +void _Lapack_multiply_matrix_hermitian(const Matrix& lhs_i, const Matrix& rhs_i, bool herm_mat_matrix_is_left, Matrix &to_add_then_result, T alpha, T beta) +{ + char herm_mat_is_left_or_right = (herm_mat_matrix_is_left ? 'L' : 'R'); + const Matrix *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(lhs->rows()); + int n = static_cast(rhs->columns()); + int k = static_cast(lhs->columns()); + int lda = std::max({1,k}); + int ldb = std::max({1,n}); + int ldc = std::max({1,m}); + if(std::is_same::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::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 void _Lapack_multiply_matrix_general(const Matrix& lhs, const Matrix& rhs, Matrix &to_add_then_result, T alpha, T beta) @@ -1541,7 +1627,13 @@ void _Lapack_multiply_matrix_general(const Matrix& lhs, const Matrix& } template -void MultiplyMatrices(const Matrix& lhs, const Matrix& rhs, Matrix &to_add_then_result, T alpha, T beta) +void MultiplyMatrices(const Matrix& lhs, + const Matrix& rhs, + Matrix &to_add_then_result, + T alpha, + T beta, + const MATRIX_TYPE& lhs_type, + const MATRIX_TYPE& rhs_type) { // std::cout<<"multiplication: "<& lhs, const Matrix& rhs, Matrix::value || std::is_same::value) { - _Lapack_multiply_matrix_general(lhs,rhs,to_add_then_result,alpha,beta); + 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::value || std::is_same::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 { @@ -1579,9 +1683,11 @@ void MultiplyMatrices(const Matrix& lhs, const Matrix& rhs, Matrix void MultiplyMatrices(const Matrix &matrixLHS, const Matrix &matrixRHS, - Matrix &result) + Matrix &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 @@ -1594,24 +1700,20 @@ RandomMatrix(const typename Matrix::size_type& rows, const MATRIX_TYPE& type = Poly::MATRIX_GENERAL) { Matrix result(rows,columns); - if(std::is_integral::value) - { + if(std::is_integral::value) { std::uniform_int_distribution 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(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(result/T(2)); } else if(type == MATRIX_GENERAL) {} - else - { + else { throw std::logic_error(_unsupported_type(__func__).c_str()); } return result; From 65e10e094ab95ecb9f118aa2eaa150a0539326f6 Mon Sep 17 00:00:00 2001 From: Samer Afach Date: Sat, 14 Jan 2017 19:00:36 +0100 Subject: [PATCH 4/4] Update Matrix.h --- include/internal/Matrix.h | 4311 +++++++++++++++++++------------------ 1 file changed, 2156 insertions(+), 2155 deletions(-) diff --git a/include/internal/Matrix.h b/include/internal/Matrix.h index 20282ef..a7c3e51 100644 --- a/include/internal/Matrix.h +++ b/include/internal/Matrix.h @@ -1,2155 +1,2156 @@ -/* - * File: Matrix.h - * Author: Samer Afach - * - * Created on 01. Oktober 2016, 17:12 - */ - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "LapackAdapters.h" -#include "StdAdapters.h" - - -#ifndef MATRIX_H -#define MATRIX_H - -#define ColMaj 0 -#define RowMaj 1 - - -/** - This class is an implementation of a matrix in mathematics. It can be used as a vector implementation as well. - */ - -namespace Poly -{ -extern std::string ComplexUnit; - -template -class Matrix; -} - - -namespace std -{ -template -std::string to_string(std::complex value) -{ - std::stringstream sstr; - - sstr << value.real(); - sstr << "+"; - sstr << value.imag(); - sstr << Poly::ComplexUnit; - return sstr.str(); -} - -template -Poly::Matrix real(const Poly::Matrix,M>& mat) -{ - Poly::Matrix result(mat.rows(),mat.columns()); - std::transform(mat.begin(),mat.end(),result.begin(),[](const std::complex& elem) -> T {return elem.real();}); -} -template -Poly::Matrix imag(const Poly::Matrix,M>& mat) -{ - Poly::Matrix result(mat.rows(),mat.columns()); - std::transform(mat.begin(),mat.end(),result.begin(),[](const std::complex& elem) -> T {return elem.imag();}); -} - -template -void swap(Poly::Matrix&__a, Poly::Matrix&__b) -{ - std::swap(__a._rows,__b._rows); - std::swap(__a._columns,__b._columns); - std::swap(__a._matEls,__b._matEls); - -} -} - -namespace Poly -{ - -template -struct ExtractType; - -template