From 2cb5258cd03fb1200cb086affe4ed70678e2a04a Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Thu, 18 Jan 2024 17:56:59 +0100
Subject: [PATCH] fixed permutation matrices

---
 mat.cc | 31 +++++++++++++++++++++++++------
 mat.h  |  3 ++-
 t.cc   | 20 +++++++++++++++++---
 3 files changed, 44 insertions(+), 10 deletions(-)

diff --git a/mat.cc b/mat.cc
index 469246e..78d0750 100644
--- a/mat.cc
+++ b/mat.cc
@@ -1834,6 +1834,25 @@ NRMat<std::complex<double> >::dot(const NRMat<std::complex<double> > &rhs) const
 	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$
  * @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)
 {
 int n=p.size();
-resize(n,n);
+nn=mm=0; count=0; v=0; resize(n,n);
 clear();
 T alpha= parity? p.parity():1;
 axpy(alpha,p,direction);
@@ -3177,18 +3196,18 @@ template<typename T>
 NRMat<T>::NRMat(const WeightPermutation<int,T> &wp, const bool direction)
 {
 int n=wp.size();
-resize(n,n);
+nn=mm=0; count=0; v=0; resize(n,n);
 clear(); 
 axpy(wp.weight,wp.perm,direction);
 }
 
 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();
-if(na<=0) laerror("cannot deduce matrix size from empty PermutationAlgebra");
-int n=ap[0].size();
-resize(n,n);
+if(na<=0 && nforce<=0) laerror("cannot deduce matrix size from empty PermutationAlgebra");
+int n= nforce>0?nforce:ap[0].size();
+nn=mm=0; count=0; v=0; resize(n,n);
 clear();
 for(int i=0; i<na; ++i) axpy(ap[i].weight,ap[i].perm,direction);
 }
diff --git a/mat.h b/mat.h
index f7c86ec..3c2086a 100644
--- a/mat.h
+++ b/mat.h
@@ -139,7 +139,7 @@ public:
 	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 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|=() 
  * @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>
 void NRMat<T>::resize(int n, int m) {
 #ifdef DEBUG
diff --git a/t.cc b/t.cc
index f8a79ec..f101b09 100644
--- a/t.cc
+++ b/t.cc
@@ -91,6 +91,7 @@ cout<<p;
 static int unitary_n;
 static PERM_RANK_TYPE space_dim;
 static NRVec<PermutationAlgebra<int,int> > allyoung;
+static NRVec<NRMat<int> >allyoungmat;
 static NRVec<int> allyoung_irrep;
 int current_irrep;
 int allyoung_index;
@@ -101,6 +102,8 @@ cout <<y;
 if(!y.is_standard()) laerror("internal error in young");
 allyoung[allyoung_index] = y.young_operator();
 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_index++;
 }
@@ -2245,6 +2248,7 @@ cout <<Sn;
 if(!Sn.is_valid()) laerror("internal error in Sn character calculation");
 cout <<"allyoung.resize "<<Sn.sumirrepdims()<<endl;
 allyoung.resize(Sn.sumirrepdims());
+allyoungmat.resize(Sn.sumirrepdims());
 allyoung_irrep.resize(Sn.sumirrepdims());
 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");}
 
 for(int i=0; i<allyoung.size(); ++i)
+	{
 	for(int j=0; j<allyoung.size(); ++j)
 		{
-		PermutationAlgebra<int,int> r=allyoung[i]*allyoung[j];
 		//cout <<"Young "<<i<<" "<<allyoung[i]<<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(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(1)
+if(0)
 {
 /*this is kucharskiP antisymmetrizer just without parsing the input
 generate antisymmetrization operator for Brandow diagrams in the Kucharski convention
@@ -3122,4 +3131,9 @@ for(int i=0; i<a.size(); ++i)
 	}
 }
 
+if(1)
+{
+}
+
+
 }