diff --git a/fourindex.h b/fourindex.h index f05b936..4339a24 100644 --- a/fourindex.h +++ b/fourindex.h @@ -716,12 +716,12 @@ if (IJ<0 || IJ>=(unsigned int)NRSMat::nn || KL<0 || KL>=(unsigned int)NRSMat< template T& fourindex_dense::operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) { -unsigned int I = SMat_index_1(i,j); -unsigned int J = SMat_index_1(k,l); +int I = SMat_index_1(i,j); +int J = SMat_index_1(k,l); //I,J act as indices of a NRSmat #ifdef DEBUG if (*NRSMat::count != 1) laerror("lval (i,j,k,l) with count > 1 in fourindex_dense"); - if (I<0 || I>=(unsigned int)NRSMat::nn || J<0 || J>=(unsigned int)NRSMat::nn) laerror("fourindex_dense index out of range"); + if (I<0 || I>=NRSMat::nn || J<0 || J>=NRSMat::nn) laerror("fourindex_dense index out of range"); if (!NRSMat::v) laerror("access to unallocated fourindex_dense"); #endif return NRSMat::v[SMat_index(I,J)]; @@ -731,11 +731,11 @@ return NRSMat::v[SMat_index(I,J)]; template const T& fourindex_dense::operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const { -unsigned int I = SMat_index_1(i,j); -unsigned int J = SMat_index_1(k,l); +int I = SMat_index_1(i,j); +int J = SMat_index_1(k,l); //I,J act as indices of a NRSmat #ifdef DEBUG - if (I<0 || I>=(unsigned int)NRSMat::nn || J<0 || J>=(unsigned int)NRSMat::nn) laerror("fourindex_dense index out of range"); + if (I<0 || I>=NRSMat::nn || J<0 || J>=NRSMat::nn) laerror("fourindex_dense index out of range"); if (!NRSMat::v) laerror("access to unallocated fourindex_dense"); #endif return NRSMat::v[SMat_index(I,J)]; diff --git a/smat.h b/smat.h index 2c7c15f..dfd7f91 100644 --- a/smat.h +++ b/smat.h @@ -336,6 +336,23 @@ inline T SMat_index_1ilej(T i, T j) return j*(j-1)/2+i-1; } +//indexing for antisymmetric matrix (used by fourindex) + +template +inline T ASMat_index(T i, T j) +{ +if(i==j) return -1; +return (i>j) ? i*(i-1)/2+j : j*(j-1)/2+i; +} + +template +inline T ASMat_index_1(T i, T j) +{ +if(i==j) return -1; +return (i>j)? (i-2)*(i-1)/2+j-1 : (j-2)*(j-1)/2+i-1; +} + + // access the element, 2-dim array case template