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

Cleaned size checks.

Browse Source 
Added functions: zeros, ones, trace, elementWiseProduct and
elementWiseProduct_inplace
renamed functions to be with _inplace when they replace "this".
This commit is contained in:
Samer Afach 2016-10-17 23:24:04 +02:00
parent f65ea9ccbd
commit 03ff1f639a
4 changed files with 120 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

216
include/internal/Matrix.h
View File

@ -116,7 +116,9 @@ void MultiplyMatrices(const Matrix<T,M> &matrixLHS, const Matrix<T,M> &matrixRHS
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);
template <typename T, int M>
Matrix<T,M> IdentityMatrix(const typename Poly::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>
void _check_equal_matrices_sizes(const Matrix<T,M>&, const Matrix<U,M>&);
template <typename T, int M>
class Matrix
@ -138,6 +140,8 @@ private:
inline size_type SizeToReserve(const size_type &rows, const size_type &columns) const;
void invert_manual();
void _check_square_matrix();
public:
Matrix<T, M>();
@ -203,16 +207,20 @@ public:
inline typename std::vector<T>::iterator end() noexcept;
inline typename std::vector<T>::const_iterator end() const noexcept;
Matrix<T,M>& operator+=(const Matrix<T,M>& rhs);
Matrix<T,M>& operator-=(const Matrix<T,M>& rhs);
Matrix<T,M>& operator*=(const Matrix<T,M>& rhs);
Matrix<T,M>& operator/=(const Matrix<T,M>& rhs);
void applyInverse();
Matrix<T,M> getInverse();
void elementWiseProduct_inplace(const Matrix<T,M>& rhs);
Matrix<T,M> elementWiseProduct(const Matrix<T,M>& rhs);
T trace();
void zeros();
void ones();
void inverse_inplace();
Matrix<T,M> inverse();
Matrix<T,M> getExp();
void applyConjugate();
void applyConjugateTranspose();
void conjugate_inplace();
void conjugateTranspose_inplace();
std::vector<T> vectorize(const VECTORIZATION_MODE& mode = VECMODE_FULL) const;
const T& at(size_type row, size_type column) const;
T& at(size_type row, size_type column);
@ -224,7 +232,7 @@ public:
const size_type& rows() const;
const size_type& columns() const;
void clear();
void applyTranspose();
void transpose_inplace();
T getDeterminant();
Matrix<T,M> getSVD();
std::string asString(int precision = 7, char open = '[', char close = ']', char sep = ',');
@ -420,13 +428,13 @@ bool Matrix<T,M>::operator!=(const Matrix<T,M>& rhs)
}
template <typename T, int M>
void Matrix<T,M>::applyTranspose()
void Matrix<T,M>::transpose_inplace()
{
*this = Transpose(*this);
}
template <typename T, int M>
void Matrix<T,M>::applyConjugateTranspose()
void Matrix<T,M>::conjugateTranspose_inplace()
{
*this = ConjugateTranspose(*this);
}
@ -509,12 +517,7 @@ Matrix<std::complex<T>,M> operator*(const Matrix<T,M>& lhs, const std::complex<T
template <typename T, int M>
Matrix<T,M> operator+(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs)
{
#ifdef POLYMATH_DEBUG
if((lhs.columns() != rhs.columns()) || (lhs.rows() != rhs.rows()))
{
throw std::length_error("Invalid matrix sizes for addition");
}
#endif
_check_equal_matrices_sizes(lhs,rhs);
Matrix<T,M> result(rhs.rows(),rhs.columns());
std::transform( lhs.begin(), lhs.end(),
rhs.begin(), result.begin(),
@ -526,12 +529,7 @@ template <typename T, int M>
Matrix<T,M> operator-(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs)
{
Matrix<T,M> result;
#ifdef POLYMATH_DEBUG
if((lhs.columns() != rhs.columns()) || (lhs.rows() != rhs.rows()))
{
throw std::length_error("Invalid matrix sizes for subtraction");
}
#endif
_check_equal_matrices_sizes(lhs,rhs);
result.resize(rhs.rows(),rhs.columns());
std::transform( lhs.begin(), lhs.end(),
rhs.begin(), result.begin(),
@ -542,12 +540,7 @@ Matrix<T,M> operator-(const Matrix<T,M>& lhs, const Matrix<T,M>& rhs)
template <typename T, int M>
Matrix<std::complex<T>,M> operator+(const Matrix<std::complex<T>,M>& lhs, const Matrix<T,M>& rhs)
{
#ifdef POLYMATH_DEBUG
if((lhs.columns() != rhs.columns()) || (lhs.rows() != rhs.rows()))
{
throw std::length_error("Invalid matrix sizes for addition");
}
#endif
_check_equal_matrices_sizes(lhs,rhs);
Matrix<std::complex<T>,M> result(rhs.rows(),rhs.columns());
std::transform( lhs.begin(), lhs.end(),
rhs.begin(), result.begin(),
@ -558,12 +551,7 @@ Matrix<std::complex<T>,M> operator+(const Matrix<std::complex<T>,M>& lhs, const
template <typename T, int M>
Matrix<std::complex<T>,M> operator+(const Matrix<T,M>& lhs, const Matrix<std::complex<T>,M>& rhs)
{
#ifdef POLYMATH_DEBUG
if((lhs.columns() != rhs.columns()) || (lhs.rows() != rhs.rows()))
{
throw std::length_error("Invalid matrix sizes for addition");
}
#endif
_check_equal_matrices_sizes(lhs,rhs);
Matrix<std::complex<T>,M> result(rhs.rows(),rhs.columns());
std::transform( lhs.begin(), lhs.end(),
rhs.begin(), result.begin(),
@ -573,12 +561,7 @@ Matrix<std::complex<T>,M> operator+(const Matrix<T,M>& lhs, const Matrix<std::co
template <typename T, int M>
Matrix<std::complex<T>,M> operator-(const Matrix<std::complex<T>,M>& lhs, const Matrix<T,M>& rhs)
{
#ifdef POLYMATH_DEBUG
if((lhs.columns() != rhs.columns()) || (lhs.rows() != rhs.rows()))
{
throw std::length_error("Invalid matrix sizes for addition");
}
#endif
_check_equal_matrices_sizes(lhs,rhs);
Matrix<std::complex<T>,M> result(rhs.rows(),rhs.columns());
std::transform( lhs.begin(), lhs.end(),
rhs.begin(), result.begin(),
@ -589,12 +572,7 @@ Matrix<std::complex<T>,M> operator-(const Matrix<std::complex<T>,M>& lhs, const
template <typename T, int M>
Matrix<std::complex<T>,M> operator-(const Matrix<T,M>& lhs, const Matrix<std::complex<T>,M>& rhs)
{
#ifdef POLYMATH_DEBUG
if((lhs.columns() != rhs.columns()) || (lhs.rows() != rhs.rows()))
{
throw std::length_error("Invalid matrix sizes for addition");
}
#endif
_check_equal_matrices_sizes(lhs,rhs);
Matrix<std::complex<T>,M> result(rhs.rows(),rhs.columns());
std::transform( lhs.begin(), lhs.end(),
rhs.begin(), result.begin(),
@ -624,8 +602,52 @@ Matrix<T,M>& Matrix<T,M>::operator*=(const Matrix<T,M>& rhs)
}
template <typename T, int M>
void Matrix<T,M>::applyInverse()
void Matrix<T,M>::elementWiseProduct_inplace(const Matrix<T,M> &rhs)
{
_check_equal_matrices_sizes(*this,rhs);
std::transform( this->begin(), this->end(),
rhs.begin(), this->begin(),
std::multiplies<T>());
}
template <typename T, int M>
Matrix<T,M> Matrix<T,M>::elementWiseProduct(const Matrix<T,M> &rhs)
{
Matrix<T,M> result(this->rows(),this->columns());
std::transform( this->begin(), this->end(),
rhs.begin(), result.begin(),
std::multiplies<T>());
return result;
}
template <typename T, int M>
T Matrix<T,M>::trace()
{
_check_square_matrix();
T result = T(0);
for(size_type i = 0; i < this->rows(); i++)
{
result += this->at(i,i);
}
return result;
}
template <typename T, int M>
void Matrix<T,M>::zeros()
{
std::fill(this->begin(), this->end(), T(0));
}
template <typename T, int M>
void Matrix<T,M>::ones()
{
std::fill(this->begin(), this->end(), T(1));
}
template <typename T, int M>
void Matrix<T,M>::inverse_inplace()
{
//FIXME: Make size error functions all in one place
if(rows() != columns())
{
throw std::length_error("Matrix inverse is only for square matrices.");
@ -714,6 +736,38 @@ void Matrix<T,M>::invert_manual()
(*this) = sideMat;
}
template <typename T, int M>
void Matrix<T,M>::_check_square_matrix()
{
#ifdef POLYMATH_DEBUG
if(this->columns() != this->rows())
{
throw std::length_error("This operation requires square matrices.");
}
#endif
}
template <typename T, typename U, int M>
void _check_equal_matrices_sizes(const Matrix<T,M>&
#ifdef POLYMATH_DEBUG
mat1
#endif
,
const Matrix<U,M>&
#ifdef POLYMATH_DEBUG
mat2
#endif
)
{
#ifdef POLYMATH_DEBUG
if((mat1.columns() != mat2.columns()) || (mat1.rows() != mat2.rows()))
{
throw std::length_error("Matrices sizes should be equal");
}
#endif
}
template <typename T, int M>
void Matrix<T,M>::copyFrom(const Matrix<T, M> &rhs)
{
@ -755,10 +809,10 @@ void Matrix<T,M>::initialize(const std::initializer_list<T> &list, Matrix<T,M>::
}
template <typename T, int M>
Matrix<T,M> Matrix<T,M>::getInverse()
Matrix<T,M> Matrix<T,M>::inverse()
{
Matrix<T,M> tmpMat = *this;
tmpMat.applyInverse();
tmpMat.inverse_inplace();
return tmpMat;
}
@ -768,7 +822,7 @@ Matrix<T,M> Matrix<T,M>::getExp()
}
template <typename T, int M>
void Matrix<T,M>::applyConjugate()
void Matrix<T,M>::conjugate_inplace()
{
transform(_matEls.begin(),_matEls.end(), _matEls.begin(), [](std::complex<typename T::value_type>& c) -> std::complex<typename T::value_type> { return std::conj(c); });
}
@ -786,12 +840,7 @@ std::vector<T> Matrix<T, M>::vectorize(const VECTORIZATION_MODE &mode) const
}
else if(mode == Poly::VECMODE_UPPER_TRIANGULAR_NO_DIAGONAL)
{
#ifdef POLYMATH_DEBUG
if(this->columns() != this->rows())
{
throw std::length_error("You cannot vectorize a non-square matrix to upper/lower triangulars.");
}
#endif
_check_square_matrix();
output.resize((this->columns()*(this->columns()-1))/2);
typename std::vector<T>::const_iterator it = _matEls.begin() + this->columns(); //iterator at rows to skip in every step, starts at second column
long size_to_copy = 0;
@ -804,12 +853,7 @@ std::vector<T> Matrix<T, M>::vectorize(const VECTORIZATION_MODE &mode) const
}
else if(mode == Poly::VECMODE_UPPER_TRIANGULAR_WITH_DIAGONAL)
{
#ifdef POLYMATH_DEBUG
if(this->columns() != this->rows())
{
throw std::length_error("You cannot vectorize a non-square matrix to upper/lower triangulars.");
}
#endif
_check_square_matrix();
output.resize((this->columns()*(this->columns()+1))/2);
typename std::vector<T>::const_iterator it = _matEls.begin();
long size_to_copy = 0;
@ -822,12 +866,7 @@ std::vector<T> Matrix<T, M>::vectorize(const VECTORIZATION_MODE &mode) const
}
else if(mode == Poly::VECMODE_LOWER_TRIANGULAR_NO_DIAGONAL)
{
#ifdef POLYMATH_DEBUG
if(this->columns() != this->rows())
{
throw std::length_error("You cannot vectorize a non-square matrix to upper/lower triangulars.");
}
#endif
_check_square_matrix();
output.resize((this->columns()*(this->columns()-1))/2);
typename std::vector<T>::const_iterator it = _matEls.begin(); //iterator at rows to skip in every step, starts at second column
long size_to_copy = 0;
@ -841,12 +880,7 @@ std::vector<T> Matrix<T, M>::vectorize(const VECTORIZATION_MODE &mode) const
}
else if(mode == Poly::VECMODE_LOWER_TRIANGULAR_WITH_DIAGONAL)
{
#ifdef POLYMATH_DEBUG
if(this->columns() != this->rows())
{
throw std::length_error("You cannot vectorize a non-square matrix to upper/lower triangulars.");
}
#endif
_check_square_matrix();
output.resize((this->columns()*(this->columns()+1))/2);
typename std::vector<T>::const_iterator it = _matEls.begin(); //iterator at rows to skip in every step, starts at second column
long size_to_copy = 0;
@ -870,12 +904,7 @@ std::vector<T> Matrix<T, M>::vectorize(const VECTORIZATION_MODE &mode) const
}
else if(mode == Poly::VECMODE_LOWER_TRIANGULAR_NO_DIAGONAL)
{
#ifdef POLYMATH_DEBUG
if(this->columns() != this->rows())
{
throw std::length_error("You cannot vectorize a non-square matrix to upper/lower triangulars.");
}
#endif
_check_square_matrix();
output.resize((this->columns()*(this->columns()-1))/2);
typename std::vector<T>::const_iterator it = _matEls.begin() + this->columns(); //iterator at rows to skip in every step, starts at second column
long size_to_copy = 0;
@ -888,12 +917,7 @@ std::vector<T> Matrix<T, M>::vectorize(const VECTORIZATION_MODE &mode) const
}
else if(mode == Poly::VECMODE_LOWER_TRIANGULAR_WITH_DIAGONAL)
{
#ifdef POLYMATH_DEBUG
if(this->columns() != this->rows())
{
throw std::length_error("You cannot vectorize a non-square matrix to upper/lower triangulars.");
}
#endif
_check_square_matrix();
output.resize((this->columns()*(this->columns()+1))/2);
typename std::vector<T>::const_iterator it = _matEls.begin();
long size_to_copy = 0;
@ -906,12 +930,7 @@ std::vector<T> Matrix<T, M>::vectorize(const VECTORIZATION_MODE &mode) const
}
else if(mode == Poly::VECMODE_UPPER_TRIANGULAR_NO_DIAGONAL)
{
#ifdef POLYMATH_DEBUG
if(this->columns() != this->rows())
{
throw std::length_error("You cannot vectorize a non-square matrix to upper/lower triangulars.");
}
#endif
_check_square_matrix();
output.resize((this->columns()*(this->columns()-1))/2);
typename std::vector<T>::const_iterator it = _matEls.begin(); //iterator at rows to skip in every step, starts at second column
long size_to_copy = 0;
@ -925,12 +944,7 @@ std::vector<T> Matrix<T, M>::vectorize(const VECTORIZATION_MODE &mode) const
}
else if(mode == Poly::VECMODE_UPPER_TRIANGULAR_WITH_DIAGONAL)
{
#ifdef POLYMATH_DEBUG
if(this->columns() != this->rows())
{
throw std::length_error("You cannot vectorize a non-square matrix to upper/lower triangulars.");
}
#endif
_check_square_matrix();
output.resize((this->columns()*(this->columns()+1))/2);
typename std::vector<T>::const_iterator it = _matEls.begin(); //iterator at rows to skip in every step, starts at second column
long size_to_copy = 0;
@ -1556,7 +1570,7 @@ int _Lapack_EigenVals_hermitian_divconq_double(Poly::Matrix<T>& eigenValues, Pol
iworkspace.get(), &liwork, &info);
if(M != ColMaj)
{
input.applyConjugateTranspose();
input.conjugateTranspose_inplace();
}
if(info > 0)
{
@ -1609,7 +1623,7 @@ int _Lapack_EigenVals_hermitian_divconq_float(Poly::Matrix<T>& eigenValues, Poly
iworkspace.get(), &liwork, &info);
if(M != ColMaj)
{
input.applyConjugateTranspose();
input.conjugateTranspose_inplace();
}
if(info > 0)
{
@ -1646,7 +1660,7 @@ int _Lapack_EigenVals_hermitian_double(Poly::Matrix<T>& eigenValues, Poly::Matri
zheev_(&jobz, &uplo, &n, (CX_TP*)&input.front(), &lda, (TP*)&eigenValues.front(), (CX_TP*)workspace.get(), &lwork, rwork.get(), &info);
if(M != ColMaj)
{
input.applyConjugateTranspose();
input.conjugateTranspose_inplace();
}
if(info > 0)
{
@ -1683,7 +1697,7 @@ int _Lapack_EigenVals_hermitian_float(Poly::Matrix<T>& eigenValues, Poly::Matrix
cheev_(&jobz, &uplo, &n, (CX_TP*)&input.front(), &lda, (TP*)&eigenValues.front(), (CX_TP*)workspace.get(), &lwork, rwork.get(), &info);
if(M != ColMaj)
{
input.applyConjugateTranspose();
input.conjugateTranspose_inplace();
}
if(info > 0)
{

1
src/Matrix.cpp
View File

@ -26,4 +26,5 @@ std::string _unsupported_type(const std::string& function_name)
{
return std::string("The type you is not supported for the function " + function_name + ".");
}
}

3
tests/tests.cpp
View File

@ -44,6 +44,9 @@ int RunTests()
int main()
{
int len = 3;
Poly::Matrix<double> mat_d = Poly::RandomMatrix<double>(len,len,0,10,std::random_device{}());
Poly::Matrix<double> mat_e = Poly::RandomMatrix<double>(len,len,0,10,std::random_device{}());
RunTests();
std::cout<<"Tests program exited with no errors."<<std::endl;
return 0;

2
tests/tests.h
View File

@ -48,7 +48,7 @@ bool TestInverse(int len, bool print)
if(print) std::cout<<mat_d<<std::endl;
if(print) std::cout<<mat_d.asString(32,'[',']',',')<<std::endl;
int prec = precision_per_type<T>();
auto mat_d_i = mat_d.getInverse();
auto mat_d_i = mat_d.inverse();
PyObject *main = PyImport_AddModule("__main__");
PyRun_SimpleString(std::string("data={}").c_str());
PyRun_SimpleString(std::string("data['a']=np.matrix(" + mat_d.asString(32,'[',']',',') + ",dtype="+python_type_per_type<T>()+")").c_str());