-Added Hermitian and symmetric matrices multiplication
-Fixed imaginary numbers printing by removing excess "+" when imaginary part is negative -
This commit is contained in:
parent
7e6076af68
commit
d09bf4d2bf
@ -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" 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 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 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 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 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
|
//#endif
|
||||||
|
|
||||||
|
|
||||||
|
@ -106,7 +106,7 @@ to_string(const T& a_value, int precision)
|
|||||||
{
|
{
|
||||||
std::ostringstream out;
|
std::ostringstream out;
|
||||||
out << std::setprecision(precision) << a_value.real();
|
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 << std::setprecision(precision) << a_value.imag();
|
||||||
out << Poly::ComplexUnit;
|
out << Poly::ComplexUnit;
|
||||||
return out.str();
|
return out.str();
|
||||||
@ -133,9 +133,9 @@ const static EIGENS_METHOD METHOD_DIVIDE_CONQUER = 1;
|
|||||||
|
|
||||||
|
|
||||||
template <typename T, int M>
|
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>
|
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>
|
template <typename T, int M>
|
||||||
Matrix<T,M> IdentityMatrix(const typename Matrix<T,M>::size_type& size);
|
Matrix<T,M> IdentityMatrix(const typename Matrix<T,M>::size_type& size);
|
||||||
template <typename T, typename U, int M>
|
template <typename T, typename U, int M>
|
||||||
@ -558,7 +558,6 @@ Matrix<std::complex<T>,M> operator*(const T& lhs, const Matrix<std::complex<T>,M
|
|||||||
return result;
|
return result;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
template <typename T, int M>
|
template <typename T, int M>
|
||||||
Matrix<T,M> operator*(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs)
|
Matrix<T,M> operator*(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs)
|
||||||
{
|
{
|
||||||
@ -1498,6 +1497,93 @@ std::string _lapack_eigens_exception_converge_error(int info);
|
|||||||
std::string _lapack_unsupported_type();
|
std::string _lapack_unsupported_type();
|
||||||
std::string _unsupported_type(const std::string& function_name);
|
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>
|
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)
|
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)
|
||||||
@ -1541,7 +1627,13 @@ void _Lapack_multiply_matrix_general(const Matrix<T,M>& lhs, const Matrix<T,M>&
|
|||||||
}
|
}
|
||||||
|
|
||||||
template <typename T, int 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;
|
// std::cout<<"multiplication: "<<lhs.columns()<<"\t"<<rhs.rows()<<std::endl;
|
||||||
#ifdef POLYMATH_DEBUG
|
#ifdef POLYMATH_DEBUG
|
||||||
@ -1555,9 +1647,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,float>::value ||
|
||||||
std::is_same<T,lapack_complex_double>::value ||
|
std::is_same<T,lapack_complex_double>::value ||
|
||||||
std::is_same<T,lapack_complex_float>::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);
|
_Lapack_multiply_matrix_general(lhs,rhs,to_add_then_result,alpha,beta);
|
||||||
}
|
}
|
||||||
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
T comp;
|
T comp;
|
||||||
@ -1579,9 +1683,11 @@ void MultiplyMatrices(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs, Matrix<T,M
|
|||||||
template <typename T, int M>
|
template <typename T, int M>
|
||||||
void MultiplyMatrices(const Matrix<T,M> &matrixLHS,
|
void MultiplyMatrices(const Matrix<T,M> &matrixLHS,
|
||||||
const Matrix<T,M> &matrixRHS,
|
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>
|
template <typename T, int M = ColMaj>
|
||||||
@ -1594,24 +1700,20 @@ RandomMatrix(const typename Matrix<T,M>::size_type& rows,
|
|||||||
const MATRIX_TYPE& type = Poly::MATRIX_GENERAL)
|
const MATRIX_TYPE& type = Poly::MATRIX_GENERAL)
|
||||||
{
|
{
|
||||||
Matrix<T,M> result(rows,columns);
|
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::uniform_int_distribution<T> distribution(min,max);
|
||||||
std::default_random_engine generator(seed);
|
std::default_random_engine generator(seed);
|
||||||
std::generate(result.begin(), result.end(), [&]() { return distribution(generator); });
|
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));
|
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));
|
result = result/T(2) - Poly::Transpose<T,M>(result/T(2));
|
||||||
}
|
}
|
||||||
else if(type == MATRIX_GENERAL) {}
|
else if(type == MATRIX_GENERAL) {}
|
||||||
else
|
else {
|
||||||
{
|
|
||||||
throw std::logic_error(_unsupported_type(__func__).c_str());
|
throw std::logic_error(_unsupported_type(__func__).c_str());
|
||||||
}
|
}
|
||||||
return result;
|
return result;
|
||||||
|
Loading…
Reference in New Issue
Block a user