tensor: indexmatrix() and zero index detection

la_traits.h

@@ -296,7 +296,7 @@ static void multiput(size_t n, int fd, const std::complex<C> *x, bool dimensions
 	  }
 	while(total < n);
 	}
-static void copy(std::complex<C> *dest, std::complex<C> *src, size_t n) {memcpy(dest,src,n*sizeof(std::complex<C>));}
+static void copy(std::complex<C> *dest, const std::complex<C> *src, size_t n) {memcpy(dest,src,n*sizeof(std::complex<C>));}
 static void clear(std::complex<C> *dest, size_t n) {memset(dest,0,n*sizeof(std::complex<C>));}
 static void copyonwrite(std::complex<C> &x) {};
 static bool is_plaindata() {return true;}
@@ -356,7 +356,7 @@ static void multiput(size_t n, int fd, const C *x, bool dimensions=0)
 	  }
 	while(total < n);
 	}
-static void copy(C *dest, C *src, size_t n) {memcpy(dest,src,n*sizeof(C));}
+static void copy(C *dest, const C *src, size_t n) {memcpy(dest,src,n*sizeof(C));}
 static void clear(C *dest, size_t n) {memset(dest,0,n*sizeof(C));}
 static void copyonwrite(C &x) {};
 static bool is_plaindata() {return true;}
@@ -396,7 +396,7 @@ static void put(int fd, const X<C> &x, bool dimensions=1, bool transp=0, bool or
 static void get(int fd, X<C> &x, bool dimensions=1, bool transp=0, bool orcaformat=false) {x.get(fd,dimensions,transp,orcaformat);} \
 static void multiput(size_t n,int fd, const X<C> *x, bool dimensions=1, bool orcaformat=false) {for(size_t i=0; i<n; ++i) x[i].put(fd,dimensions,false,orcaformat);} \
 static void multiget(size_t n,int fd, X<C> *x, bool dimensions=1, bool orcaformat=false) {for(size_t i=0; i<n; ++i) x[i].get(fd,dimensions,false,orcaformat);} \
-static void copy(X<C> *dest, X<C> *src, size_t n) {for(size_t i=0; i<n; ++i) dest[i]=src[i];} \
+static void copy(X<C> *dest, const X<C> *src, size_t n) {for(size_t i=0; i<n; ++i) dest[i]=src[i];} \
 static void clear(X<C> *dest, size_t n) {for(size_t i=0; i<n; ++i) dest[i].clear();}\
 static void copyonwrite(X<C> &x) {x.copyonwrite();}\
 static bool is_plaindata() {return false;}\
@@ -439,7 +439,7 @@ static void put(int fd, const X<C> &x, bool dimensions=1, bool transp=0, bool or
 static void get(int fd, X<C> &x, bool dimensions=1, bool transp=0, bool orcaformat=false) {x.get(fd,dimensions,false,orcaformat);}  \
 static void multiput(size_t n,int fd, const X<C> *x, bool dimensions=1, bool orcaformat=false) {for(size_t i=0; i<n; ++i) x[i].put(fd,dimensions,false,orcaformat);}  \
 static void multiget(size_t n,int fd, X<C> *x, bool dimensions=1, bool orcaformat=false) {for(size_t i=0; i<n; ++i) x[i].get(fd,dimensions,false,orcaformat);}  \
-static void copy(C *dest, C *src, size_t n) {for(size_t i=0; i<n; ++i) dest[i]=src[i];}  \
+static void copy(C *dest, const C *src, size_t n) {for(size_t i=0; i<n; ++i) dest[i]=src[i];}  \
 static void clear(C *dest, size_t n) {for(size_t i=0; i<n; ++i) dest[i].clear();} \
 static void copyonwrite(X<C> &x) {x.copyonwrite();} \
 static bool is_plaindata() {return false;}\

mat.cc

@@ -113,6 +113,29 @@ const NRVec<T> NRMat<T>::row(const int i, int l) const {
 	return r;
 }
 
+
+/***************************************************************************//**
+ * store given row of this matrix of general type <code>T</code>
+ * @param[in] i row index starting from zero
+ * @param[in] l consider this value as the count of columns
+ ******************************************************************************/
+template <typename T>
+void NRMat<T>::rowset(const NRVec<T> &r, const int i, int l) {
+#ifdef DEBUG
+	if(i < 0 || i >= nn) laerror("illegal index");
+#endif
+	if(l < 0) l = mm;
+	LA_traits<T>::copy(
+#ifdef MATPTR
+	        v[i]
+#else
+		v + i*(size_t)l
+#endif
+		, &r[0], l);
+}
+
+
+
 /***************************************************************************//**
  * routine for raw output
  * @param[in] fd file descriptor for output

mat.h

@@ -266,6 +266,9 @@ public:
 	//! get the i<sup>th</sup> row
 	const NRVec<T> row(const int i, int l = -1) const;
 
+	//! set the i<sup>th</sup> row
+	void rowset(const NRVec<T> &r, const int i, int l = -1);
+
 	//! get the j<sup>th</sup> column
 	const NRVec<T> column(const int j, int l = -1) const {
 		NOT_GPU(*this);
@@ -275,6 +278,13 @@ public:
 		return r;
 	};
 
+	//! set the j<sup>th</sup> column
+	void columnset(const NRVec<T> &r, const int j, int l = -1) {
+		NOT_GPU(*this);
+		if(l < 0) l = nn;
+		for(register int i=0; i<l; ++i) (*this)(i,j) = r[i];
+	};
+
 	//! extract the digonal elements of this matrix and store them into a vector
 	const T* diagonalof(NRVec<T> &, const bool divide = 0, bool cache = false) const;
 	//! set diagonal elements

t.cc

@@ -4597,6 +4597,8 @@ Tensor<double> x(shape); x.randomize(1.);
 x.defaultnames();
 cout <<"x= "<<x.shape << " "<<x.names<<endl;
 
+cout <<"indexmatrix of x = "<<x.indexmatrix();
+
 NRVec<INDEXGROUP> yshape(2);
         yshape[0].number=1; 
         yshape[0].symmetry=0;

tensor.cc

@@ -417,7 +417,7 @@ calcsize();
 
 
 template<typename T>
-void loopingroups(Tensor<T> &t, int ngroup, int igroup, T **p, SUPERINDEX &I, void (*callback)(const SUPERINDEX &, T *))
+void loopingroups(Tensor<T> &t, int ngroup, int igroup, T **p, SUPERINDEX &I, void (*callback)(const SUPERINDEX &, T *), bool skipzeros)
 {
 LA_index istart,iend;
 const INDEXGROUP *sh = &(* const_cast<const NRVec<INDEXGROUP> *>(&t.shape))[ngroup];
@@ -443,6 +443,8 @@ switch(sh->symmetry)
 
 for(LA_index i = istart; i<=iend; ++i)
 	{
+	if(skipzeros && i==0) continue;
+
 	I[ngroup][igroup]=i;
 	if(ngroup==0 && igroup==0) 
 		{
@@ -460,14 +462,14 @@ for(LA_index i = istart; i<=iend; ++i)
 			const INDEXGROUP *sh2 = &(* const_cast<const NRVec<INDEXGROUP> *>(&t.shape))[newngroup];
 			newigroup=sh2->number-1;
 			}
-		loopingroups(t,newngroup,newigroup,p,I,callback);
+		loopingroups(t,newngroup,newigroup,p,I,callback,skipzeros);
 		}
 	}
 }
 
 
 template<typename T>
-void Tensor<T>::loopover(void (*callback)(const SUPERINDEX &, T *))
+void Tensor<T>::loopover(void (*callback)(const SUPERINDEX &, T *), bool skipzeros)
 {
 SUPERINDEX I(shape.size());
 for(int i=0; i<I.size(); ++i)
@@ -479,12 +481,12 @@ for(int i=0; i<I.size(); ++i)
 T *pp=&data[0];
 int ss=shape.size()-1;
 const INDEXGROUP *sh = &(* const_cast<const NRVec<INDEXGROUP> *>(&shape))[ss];
-loopingroups(*this,ss,sh->number-1,&pp,I,callback);
+loopingroups(*this,ss,sh->number-1,&pp,I,callback,skipzeros);
 }
 
 
 template<typename T>
-void constloopingroups(const Tensor<T> &t, int ngroup, int igroup, const T **p, SUPERINDEX &I, void (*callback)(const SUPERINDEX &, const T *))
+void constloopingroups(const Tensor<T> &t, int ngroup, int igroup, const T **p, SUPERINDEX &I, void (*callback)(const SUPERINDEX &, const T *), bool skipzeros)
 {
 LA_index istart,iend;
 const INDEXGROUP *sh = &t.shape[ngroup];
@@ -510,6 +512,8 @@ switch(sh->symmetry)
 
 for(LA_index i = istart; i<=iend; ++i)
 	{
+        if(skipzeros && i==0) continue;
+
 	I[ngroup][igroup]=i;
 	if(ngroup==0 && igroup==0) 
 		{
@@ -527,14 +531,14 @@ for(LA_index i = istart; i<=iend; ++i)
 			const INDEXGROUP *sh2 = &(* const_cast<const NRVec<INDEXGROUP> *>(&t.shape))[newngroup];
 			newigroup=sh2->number-1;
 			}
-		constloopingroups(t,newngroup,newigroup,p,I,callback);
+		constloopingroups(t,newngroup,newigroup,p,I,callback,skipzeros);
 		}
 	}
 }
 
 
 template<typename T>
-void Tensor<T>::constloopover(void (*callback)(const SUPERINDEX &, const T *)) const
+void Tensor<T>::constloopover(void (*callback)(const SUPERINDEX &, const T *), bool skipzeros) const
 {
 SUPERINDEX I(shape.size());
 for(int i=0; i<I.size(); ++i)
@@ -546,7 +550,28 @@ for(int i=0; i<I.size(); ++i)
 const T *pp=&data[0];
 int ss=shape.size()-1;
 const INDEXGROUP *sh = &shape[ss];
-constloopingroups(*this,ss,sh->number-1,&pp,I,callback);
+constloopingroups(*this,ss,sh->number-1,&pp,I,callback,skipzeros);
+}
+
+static INDEXMATRIX *indexmat_p;
+static LA_largeindex indexmat_row;
+
+template<typename T>
+static void indexmatrix_callback(const SUPERINDEX &I, const T *p)
+{
+FLATINDEX f=superindex2flat(I);
+indexmat_p->rowset(f,indexmat_row++);
+}
+
+template<typename T>
+INDEXMATRIX  Tensor<T>::indexmatrix() const
+{
+INDEXMATRIX r;
+r.resize(size(),rank());
+indexmat_p = &r;
+indexmat_row = 0;
+constloopover(indexmatrix_callback<T>,false);
+return r;
 }
 
 
@@ -2439,6 +2464,27 @@ return true;
 }
 
 
+bool zero_in_index(const FLATINDEX &I)
+{
+for(int i=0; i<I.size(); ++i) if(I[i]==0) return true;
+return false;
+}
+
+bool zero_in_index(const INDEXMATRIX &m, const LA_largeindex row) 
+{
+const LA_index *p = &m(row,0);
+for(int i=0; i<m.ncols(); ++i) if(p[i]==0) return true;
+return false;
+}
+
+
+
+bool zero_in_index(const SUPERINDEX &I)
+{
+for(int i=0; i<I.size(); ++i) if(zero_in_index(I[i])) return true;
+return false;
+}
+
 
 template class Tensor<double>;
 template class Tensor<std::complex<double> >;

tensor.h

@@ -190,8 +190,9 @@ class LA_traits<INDEXGROUP> {
 
 
 typedef NRVec<LA_index> FLATINDEX; //all indices but in a single vector
-typedef NRVec<NRVec<LA_index> > SUPERINDEX; //all indices in the INDEXGROUP structure
+typedef NRVec<FLATINDEX> SUPERINDEX; //all indices in the INDEXGROUP structure
 typedef NRVec<LA_largeindex> GROUPINDEX; //set of indices in the symmetry groups
+typedef NRMat<LA_index> INDEXMATRIX; //list of FLATINDEXes (rows of the matrix) of all tensor elements - convenient to be able to run over the whole tensor in a for loop rather than via recursive loopovers with a callback
 struct  INDEX
 {
 int group;
@@ -230,6 +231,11 @@ int flatposition(int group, int index, const NRVec<INDEXGROUP> &shape);
 int flatposition(const INDEX &i, const NRVec<INDEXGROUP> &shape); //position of that index in FLATINDEX
 INDEX indexposition(int flatindex, const NRVec<INDEXGROUP> &shape); //inverse to flatposition
 
+//useful for negative offsets and 0 index excluded
+bool zero_in_index(const FLATINDEX &);
+bool zero_in_index(const SUPERINDEX &);
+bool zero_in_index(const INDEXMATRIX &, const LA_largeindex row);
+
 FLATINDEX superindex2flat(const SUPERINDEX &I);
 
 template<typename T>
@@ -361,11 +367,14 @@ public:
 	bool fulfills_hermiticity() const; //check it is so
 	inline void randomize(const typename LA_traits<T>::normtype &x) {data.randomize(x); enforce_hermiticity();};
 
-	void loopover(void (*callback)(const SUPERINDEX &, T *)); //loop over all elements
+
-	void constloopover(void (*callback)(const SUPERINDEX &, const T *)) const; //loop over all elements
+	void loopover(void (*callback)(const SUPERINDEX &, T *), bool skipzeros=false); //loop over all elements, optionally skip zero indices (i.e. run over ...-2,-1,1,2...) which is useful for special applications
+	void constloopover(void (*callback)(const SUPERINDEX &, const T *), bool skipzeros=false) const; //loop over all elements
 	void grouploopover(void (*callback)(const GROUPINDEX &, T *)); //loop over all elements disregarding the internal structure of index groups
 	void constgrouploopover(void (*callback)(const GROUPINDEX &, const T *)) const; //loop over all elements disregarding the internal structure of index groups
 
+	INDEXMATRIX indexmatrix() const; //get indexmatrix - rows store FLATINDEXes matching data[]
+
 	Tensor permute_index_groups(const NRPerm<int> &p) const; //rearrange the tensor storage permuting index groups as a whole: result_i = source_p_i
 	Tensor permute_index_groups(const NRVec<INDEXNAME> &names) const; //permute to requested order of group's first indices (or permute individual indices of a flat tensor)