From 60e8a379f5f8fe7f43b64feffeacc91a743886fd Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Thu, 13 May 2021 16:45:10 +0200 Subject: [PATCH] progressing implementation of permutations --- mat.cc | 53 +++++++++++++++++++++++++++++++++++++++++++ mat.h | 11 +++++++++ permutation.cc | 61 +++++++++++++++++++++++++++++++++++++++++++++++++- permutation.h | 15 +++++++++---- smat.cc | 18 +++++++++++++++ smat.h | 6 +++++ t.cc | 16 ++++++++++++- vec.cc | 16 +++++++++++++ vec.h | 11 +++++++++ 9 files changed, 201 insertions(+), 6 deletions(-) diff --git a/mat.cc b/mat.cc index 350ef1c..0d972f6 100644 --- a/mat.cc +++ b/mat.cc @@ -3075,6 +3075,59 @@ NRMat& NRMat::swap_rows_cols(){ return *this; } +//apply permutations +template +const NRMat NRMat::permute_rows(const NRPerm &p) const +{ +#ifdef DEBUG +if(!p.is_valid()) laerror("invalid permutation of matrix"); +#endif +int n=p.size(); +if(n!=nn) laerror("incompatible permutation and matrix"); +#ifdef CUDALA + if(this->getlocation() != cpu || p.getlocation() != cpu ) laerror("permutations can be done only in CPU memory"); +#endif +NRMat r(nn,mm); +for(int i=1; i<=n; ++i) {int pi=p[i]-1; for(int j=0; j +const NRMat NRMat::permute_cols(const NRPerm &p) const +{ +#ifdef DEBUG +if(!p.is_valid()) laerror("invalid permutation of matrix"); +#endif +int n=p.size(); +if(n!=mm) laerror("incompatible permutation and matrix"); +#ifdef CUDALA + if(this->getlocation() != cpu || p.getlocation() != cpu ) laerror("permutations can be done only in CPU memory"); +#endif +NRMat r(nn,mm); +for(int i=1; i<=n; ++i) {int pi=p[i]-1; for(int j=0; j +const NRMat NRMat::permute_both(const NRPerm &p, const NRPerm &q) const +{ +#ifdef DEBUG +if(!p.is_valid() || !q.is_valid() ) laerror("invalid permutation of matrix"); +#endif +int n=p.size(); +int m=q.size(); +if(n!=nn ||m!=mm) laerror("incompatible permutation and matrix"); +#ifdef CUDALA + if(this->getlocation() != cpu || p.getlocation() != cpu ) laerror("permutations can be done only in CPU memory"); +#endif +NRMat r(nn,mm); +for(int i=1; i<=n; ++i) {int pi=p[i]-1; for(int j=1; j<=m; ++j) r(i-1,j-1) = (*this)(pi,q[j]-1);} +return r; +} + + + + /***************************************************************************//** * forced instantization in the corresponding object file ******************************************************************************/ diff --git a/mat.h b/mat.h index 0fd28fe..34198d2 100644 --- a/mat.h +++ b/mat.h @@ -24,6 +24,10 @@ namespace LA { +//forward declaration +template class NRPerm; + + /***************************************************************************//** * \brief NRMat class template implementing the matrix interface * @see NRVec, NRSMat @@ -111,6 +115,12 @@ public: //! ensure that the data of this matrix are referenced exactly once void copyonwrite(bool detachonly=false); + //! permute matrix elements + const NRMat permute_rows(const NRPerm &p) const; + const NRMat permute_cols(const NRPerm &p) const; + const NRMat permute_both(const NRPerm &p, const NRPerm &q) const; + + /***************************************************************************//** * routines for CUDA related stuff * \li getlocation() gets the protected data member location @@ -376,6 +386,7 @@ public: #include "smat.h" #include "sparsemat.h" #include "sparsesmat.h" +#include "permutation.h" namespace LA { diff --git a/permutation.cc b/permutation.cc index 6106002..cee6b63 100644 --- a/permutation.cc +++ b/permutation.cc @@ -19,6 +19,32 @@ #include "permutation.h" namespace LA { +template +void NRPerm::identity() +{ +T n=this->size(); +#ifdef DEBUG + if(n<0) laerror("invalid permutation size"); +#endif +if(n==0) return; +for(T i=1; i<=n; ++i) (*this)[i]=i; +} + + +template +bool NRPerm::is_identity() const +{ +T n=this->size(); +#ifdef DEBUG + if(n<0) laerror("invalid permutation size"); +#endif +if(n==0) return 1; +for(T i=1; i<=n; ++i) if((*this)[i]!=i) return 0; +return 1; +} + + + template bool NRPerm::is_valid() const { @@ -39,6 +65,10 @@ return 1; template NRPerm NRPerm::inverse() const { +#ifdef DEBUG + if(!this->is_valid()) laerror("inverse of invalid permutation"); +#endif + NRPerm q(this->size()); for(T i=1; i<=this->size(); ++i) q[(*this)[i]]=i; return q; @@ -48,15 +78,44 @@ return q; template NRPerm NRPerm::operator*(const NRPerm q) const { +#ifdef DEBUG + if(!this->is_valid() || !q.is_valid()) laerror("multiplication of invalid permutation"); +#endif + T n=this->size(); if(n!=q.size()) laerror("product of incompatible permutations"); NRPerm r(n); -for(T i=1; i<=n; ++i) r[i] = q[(*this)[i]]; +for(T i=1; i<=n; ++i) r[i] = (*this)[q[i]]; return r; } +template +NRPerm NRPerm::conjugate_by(const NRPerm q) const +{ +#ifdef DEBUG + if(!this->is_valid() || !q.is_valid()) laerror("multiplication of invalid permutation"); +#endif +T n=this->size(); +if(n!=q.size()) laerror("product of incompatible permutations"); +NRPerm qi=q.inverse(); +NRPerm r(n); +for(T i=1; i<=n; ++i) r[i] = qi[(*this)[q[i]]]; +return r; + +} + +template +int NRPerm::parity() const +{ +if(!this->is_valid()) return 0; +T n=this->size(); +if(n==1) return 1; +T count=0; +for(T i=2;i<=n;i++) for(T j=1;j(*this)[i]) count++; +return (count&1)? -1:1; +} diff --git a/permutation.h b/permutation.h index b6e714f..1839fe6 100644 --- a/permutation.h +++ b/permutation.h @@ -28,6 +28,9 @@ namespace LA { +//forward declaration +template class NRVec_from1; + template class NRPerm : public NRVec_from1 { public: @@ -39,17 +42,21 @@ public: NRPerm(const T *a, const int n): NRVec_from1(a, n) {}; //specific operations + void identity(); bool is_valid() const; //is it really a permutation + bool is_identity() const; NRPerm inverse() const; - NRPerm operator*(const NRPerm rhs) const; + NRPerm operator*(const NRPerm q) const; //q is rhs and applied first, this applied second + NRPerm conjugate_by(const NRPerm q) const; //q^-1 p q + int parity() const; + //TODO: - //@@@conjugate by q + //@@@permutation matrix //@@@permgener - //@@@lex rank //@@@next permutation + //@@@lex rank //@@@inversion tables - //@@@parity //@@@conversion to cycle structure and back }; diff --git a/smat.cc b/smat.cc index 178613c..75101bb 100644 --- a/smat.cc +++ b/smat.cc @@ -303,6 +303,24 @@ void NRSMat::fscanf(FILE *f, const char *format) { laerror("NRSMat::fscanf(FILE *, const char *) - unable to read matrix element"); } +//apply permutation +template +const NRSMat NRSMat::permute(const NRPerm &p) const +{ +#ifdef DEBUG +if(!p.is_valid()) laerror("invalid permutation of smatrix"); +#endif +int n=p.size(); +if(n!=(*this).size()) laerror("incompatible permutation and smatrix"); +#ifdef CUDALA + if(this->getlocation() != cpu || p.getlocation() != cpu ) laerror("permutations can be done only in CPU memory"); +#endif +NRSMat r(n); +for(int i=1; i<=n; ++i) {int pi = p[i]-1; r(i-1,i-1) = (*this)(pi,pi);} +return r; +} + + /***************************************************************************//** * multiply this real double-precision symmetric matrix \f$S\f$ stored in packed form diff --git a/smat.h b/smat.h index aba27a1..d6b4a09 100644 --- a/smat.h +++ b/smat.h @@ -27,6 +27,9 @@ namespace LA { #define NN2 ((size_t)nn*(nn+1)/2) +//forward declaration +template class NRPerm; + /***************************************************************************//** * This class implements a general symmetric or hermitian matrix the elements @@ -89,6 +92,8 @@ public: //! assign scalar value to diagonal elements NRSMat & operator=(const T &a); + //! permute matrix elements + const NRSMat permute(const NRPerm &p) const; inline int getcount() const {return count?*count:0;} @@ -176,6 +181,7 @@ public: //due to mutual includes this has to be after full class declaration #include "vec.h" #include "mat.h" +#include "permutation.h" namespace LA { diff --git a/t.cc b/t.cc index 984412c..25d5d59 100644 --- a/t.cc +++ b/t.cc @@ -22,6 +22,7 @@ #include "la.h" #include "vecmat3.h" #include "quaternion.h" +#include "permutation.h" using namespace std; using namespace LA_Vecmat3; @@ -1980,11 +1981,24 @@ cout <<"normquat2euler test "< a,b,c; cin >>a>>b; c=a+b; cout< p; +cin >>p; +int n=p.size(); +NRVec_from1 v(n); +int i; +for(i=1; i<=n; ++i) v[i]=10.*i; +cout < > complexify(const NRVec &rhs) { return r; } +template +const NRVec NRVec::permute(const NRPerm &p) const +{ +#ifdef DEBUG +if(!p.is_valid()) laerror("invalid permutation of vector"); +#endif +int n=p.size(); +if(n!=(*this).size()) laerror("incompatible permutation and vector"); +#ifdef CUDALA + if(this->getlocation() != cpu || p.getlocation() != cpu ) laerror("permutations can be done only in CPU memory"); +#endif +NRVec r(n); +for(int i=1; i<=n; ++i) r[i-1] = v[p[i]-1]; +return r; +} + /***************************************************************************//** * forced instantization in the corespoding object file ******************************************************************************/ diff --git a/vec.h b/vec.h index 47ba7d6..09da04a 100644 --- a/vec.h +++ b/vec.h @@ -30,6 +30,8 @@ namespace LA { template void lawritemat(FILE *file, const T *a, int r, int c, const char *form0, int nodim, int modulo, int issym); +template class NRPerm; + /***************************************************************************//** * static constants used in several cblas-routines ******************************************************************************/ @@ -260,6 +262,9 @@ public: return sum; }; + //! permute vector elements + const NRVec permute(const NRPerm &p) const; + //! compute the sum of the absolute values of the elements of this vector inline const typename LA_traits::normtype asum() const; @@ -382,6 +387,12 @@ public: inline T& operator[] (const int i); }; +}//namespace +//needs NRVec_from1 +#include "permutation.h" +namespace LA { + + /***************************************************************************//** * indexing operator giving the element at given position with range checking in * the DEBUG mode