fixed permutation matrices

This commit is contained in:
Jiri Pittner 2024-01-18 17:56:59 +01:00
parent 680fa93425
commit 2cb5258cd0
3 changed files with 44 additions and 10 deletions

31
mat.cc
View File

@ -1834,6 +1834,25 @@ NRMat<std::complex<double> >::dot(const NRMat<std::complex<double> > &rhs) const
return ret; return ret;
} }
template<>
const NRMat<int> NRMat<int>::operator*(const NRMat<int> &rhs) const {
#ifdef DEBUG
if(mm != rhs.nn) laerror("incompatible matrices in NRMat<int>::operator*(const NRMat<int>&)");
if(mm<=0 ||rhs.mm <= 0) laerror("illegal matrix dimension in gemm");
#endif
NOT_GPU(rhs);
NOT_GPU(*this);
NRMat<int> result(nn, rhs.mm);
result.clear();
for(int i=0;i<nn;++i)
for(int j=0; j<mm; ++j)
for(int k=0; k<rhs.mm; ++k)
result(i,k) += (*this)(i,j)*rhs(j,k);
return result;
}
/***************************************************************************//** /***************************************************************************//**
* compute product of this matrix \f$A\f$ with given real matrix \f$B\f$ * compute product of this matrix \f$A\f$ with given real matrix \f$B\f$
* @param[in] rhs matrix \f$B\f$ * @param[in] rhs matrix \f$B\f$
@ -3167,7 +3186,7 @@ template<typename T>
NRMat<T>::NRMat(const NRPerm<int> &p, const bool direction, const bool parity) NRMat<T>::NRMat(const NRPerm<int> &p, const bool direction, const bool parity)
{ {
int n=p.size(); int n=p.size();
resize(n,n); nn=mm=0; count=0; v=0; resize(n,n);
clear(); clear();
T alpha= parity? p.parity():1; T alpha= parity? p.parity():1;
axpy(alpha,p,direction); axpy(alpha,p,direction);
@ -3177,18 +3196,18 @@ template<typename T>
NRMat<T>::NRMat(const WeightPermutation<int,T> &wp, const bool direction) NRMat<T>::NRMat(const WeightPermutation<int,T> &wp, const bool direction)
{ {
int n=wp.size(); int n=wp.size();
resize(n,n); nn=mm=0; count=0; v=0; resize(n,n);
clear(); clear();
axpy(wp.weight,wp.perm,direction); axpy(wp.weight,wp.perm,direction);
} }
template<typename T> template<typename T>
NRMat<T>::NRMat(const PermutationAlgebra<int,T> &ap, const bool direction) NRMat<T>::NRMat(const PermutationAlgebra<int,T> &ap, const bool direction , int nforce)
{ {
int na= ap.size(); int na= ap.size();
if(na<=0) laerror("cannot deduce matrix size from empty PermutationAlgebra"); if(na<=0 && nforce<=0) laerror("cannot deduce matrix size from empty PermutationAlgebra");
int n=ap[0].size(); int n= nforce>0?nforce:ap[0].size();
resize(n,n); nn=mm=0; count=0; v=0; resize(n,n);
clear(); clear();
for(int i=0; i<na; ++i) axpy(ap[i].weight,ap[i].perm,direction); for(int i=0; i<na; ++i) axpy(ap[i].weight,ap[i].perm,direction);
} }

3
mat.h
View File

@ -139,7 +139,7 @@ public:
void axpy(const T alpha, const NRPerm<int> &p, const bool direction); void axpy(const T alpha, const NRPerm<int> &p, const bool direction);
explicit NRMat(const NRPerm<int> &p, const bool direction, const bool parity=false); //permutation matrix explicit NRMat(const NRPerm<int> &p, const bool direction, const bool parity=false); //permutation matrix
explicit NRMat(const WeightPermutation<int,T> &p, const bool direction); explicit NRMat(const WeightPermutation<int,T> &p, const bool direction);
explicit NRMat(const PermutationAlgebra<int,T> &p, const bool direction); explicit NRMat(const PermutationAlgebra<int,T> &p, const bool direction, const int nforce=0); //note that one cannot represent e.g. young projectors in this way, since the representation of S(n) by permutation matrices is reducible just to two irreps [n] and [n-1,1]
/***************************************************************************//** /***************************************************************************//**
@ -1167,6 +1167,7 @@ void NRMat<T>::copyonwrite(bool detachonly) {
* @see count, NRMat<T>::copyonwrite(), NRMat<T>::operator|=() * @see count, NRMat<T>::copyonwrite(), NRMat<T>::operator|=()
* @return reference to the newly copied matrix * @return reference to the newly copied matrix
******************************************************************************/ ******************************************************************************/
//NOTE it must not be used in constructors when the data were not initialized yet
template <typename T> template <typename T>
void NRMat<T>::resize(int n, int m) { void NRMat<T>::resize(int n, int m) {
#ifdef DEBUG #ifdef DEBUG

20
t.cc
View File

@ -91,6 +91,7 @@ cout<<p;
static int unitary_n; static int unitary_n;
static PERM_RANK_TYPE space_dim; static PERM_RANK_TYPE space_dim;
static NRVec<PermutationAlgebra<int,int> > allyoung; static NRVec<PermutationAlgebra<int,int> > allyoung;
static NRVec<NRMat<int> >allyoungmat;
static NRVec<int> allyoung_irrep; static NRVec<int> allyoung_irrep;
int current_irrep; int current_irrep;
int allyoung_index; int allyoung_index;
@ -101,6 +102,8 @@ cout <<y;
if(!y.is_standard()) laerror("internal error in young"); if(!y.is_standard()) laerror("internal error in young");
allyoung[allyoung_index] = y.young_operator(); allyoung[allyoung_index] = y.young_operator();
cout <<"Young "<<allyoung_index<<" (irrep "<<current_irrep<<") = "<<allyoung[allyoung_index]<<endl; cout <<"Young "<<allyoung_index<<" (irrep "<<current_irrep<<") = "<<allyoung[allyoung_index]<<endl;
allyoungmat[allyoung_index] = NRMat<int>(allyoung[allyoung_index],false);
cout <<"Matrix representation = "<<allyoungmat[allyoung_index];
allyoung_irrep[allyoung_index]=current_irrep; allyoung_irrep[allyoung_index]=current_irrep;
allyoung_index++; allyoung_index++;
} }
@ -2245,6 +2248,7 @@ cout <<Sn;
if(!Sn.is_valid()) laerror("internal error in Sn character calculation"); if(!Sn.is_valid()) laerror("internal error in Sn character calculation");
cout <<"allyoung.resize "<<Sn.sumirrepdims()<<endl; cout <<"allyoung.resize "<<Sn.sumirrepdims()<<endl;
allyoung.resize(Sn.sumirrepdims()); allyoung.resize(Sn.sumirrepdims());
allyoungmat.resize(Sn.sumirrepdims());
allyoung_irrep.resize(Sn.sumirrepdims()); allyoung_irrep.resize(Sn.sumirrepdims());
allyoung_index=0; allyoung_index=0;
@ -2257,15 +2261,20 @@ if(tot!=partitions(n)) laerror("internal error in partition generation or enumer
if(space_dim!=longpow(unitary_n,n)) {cout<<space_dim<<" "<<ipow(unitary_n,n)<<endl;laerror("integer overflow or internal error in space dimensions");} if(space_dim!=longpow(unitary_n,n)) {cout<<space_dim<<" "<<ipow(unitary_n,n)<<endl;laerror("integer overflow or internal error in space dimensions");}
for(int i=0; i<allyoung.size(); ++i) for(int i=0; i<allyoung.size(); ++i)
{
for(int j=0; j<allyoung.size(); ++j) for(int j=0; j<allyoung.size(); ++j)
{ {
PermutationAlgebra<int,int> r=allyoung[i]*allyoung[j];
//cout <<"Young "<<i<<" "<<allyoung[i]<<endl; //cout <<"Young "<<i<<" "<<allyoung[i]<<endl;
//cout <<"Young "<<j<<" "<<allyoung[j]<<endl; //cout <<"Young "<<j<<" "<<allyoung[j]<<endl;
cout <<"Product of Young "<<i<<" and "<<j<<" = "<<r<<"\n"; PermutationAlgebra<int,int> r=allyoung[i]*allyoung[j];
NRMat<int> rm(r,false,n);
NRMat<int> rm2 = allyoungmat[i]*allyoungmat[j];
cout <<"Product of Young "<<i<<" and "<<j<<" = "<<r<<"\n"<<"matrix "<<rm<<endl;
if(rm!=rm2) laerror("internal error in matrix representation of permutationalgebra");
if(i!=j && !r.is_zero()) cout <<"NONORTHOGONAL Young operators found "<<i<< " "<<j<<" (irreps "<<allyoung_irrep[i]<<" "<<allyoung_irrep[j]<<")\n"; if(i!=j && !r.is_zero()) cout <<"NONORTHOGONAL Young operators found "<<i<< " "<<j<<" (irreps "<<allyoung_irrep[i]<<" "<<allyoung_irrep[j]<<")\n";
if(allyoung_irrep[i]!=allyoung_irrep[j] && !r.is_zero()) laerror("internal error in PermutationAlgebra"); if(allyoung_irrep[i]!=allyoung_irrep[j] && !r.is_zero()) laerror("internal error in PermutationAlgebra");
} }
}
} }
@ -3093,7 +3102,7 @@ for(int i=0; i<a.size(); ++i)
if(a[i].weight!=a[i].perm.parity()) laerror("internal error in parity"); if(a[i].weight!=a[i].perm.parity()) laerror("internal error in parity");
} }
if(1) if(0)
{ {
/*this is kucharskiP antisymmetrizer just without parsing the input /*this is kucharskiP antisymmetrizer just without parsing the input
generate antisymmetrization operator for Brandow diagrams in the Kucharski convention generate antisymmetrization operator for Brandow diagrams in the Kucharski convention
@ -3122,4 +3131,9 @@ for(int i=0; i<a.size(); ++i)
} }
} }
if(1)
{
}
} }