cleanup debug prints

bitvector.h

@@ -19,9 +19,11 @@
 #ifndef _BITVECTOR_H_
 #define _BITVECTOR_H_
 
+#include "la_traits.h"
 #include "vec.h"
 #include "numbers.h"
 #include "laerror.h"
+
 #include <stdint.h>
 
 //TODO: if efficiency is requires, make also a monic_bitvector, which will not store the leading 1 explicitly
@@ -77,27 +79,39 @@ public:
 	void set(const unsigned int i) 
 		{
 #ifdef DEBUG
-		if(i>=size()) laerror("bitvector index out of range in");
+		if(i>=size()) 
+			{
+			std::cout<<"bitvector index 0 < "<<i<<" "<<size()<<std::endl;
+			laerror("bitvector index out of range in");
+			}
 #endif
 		v[i/blockbits] |= (1UL<<(i%blockbits));
 		};
 	void reset(const unsigned int i) 
 		{
 #ifdef DEBUG    
-                if(i>=size()) laerror("bitvector index out of range in");
+if(i>=size()) 
+                        {
+                        std::cout<<"bitvector index 0 < "<<i<<" "<<size()<<std::endl;
+                        laerror("bitvector index out of range in");
+                        }
 #endif
 		v[i/blockbits] &= ~(1UL<<(i%blockbits));
 		};
 	void flip(const unsigned int i) 
 		{
 #ifdef DEBUG    
-                if(i>=size()) laerror("bitvector index out of range in");
+		if(i>=size()) 
+                        {
+                        std::cout<<"bitvector index 0 < "<<i<<" "<<size()<<std::endl;
+                        laerror("bitvector index out of range in");
+                        }
 #endif
 		v[i/blockbits] ^= (1UL<<(i%blockbits));
 		};
 	const bool assign(const unsigned int i, const bool r) {if(r) set(i); else reset(i); return r;};
 	void clear() {copyonwrite(true); memset(v,0,nn*sizeof(bitvector_block));};
-	void fill() {memset(v,0xff,nn*sizeof(bitvector_block));};
+	void fill() {copyonwrite(true); memset(v,0xff,nn*sizeof(bitvector_block));};
 	void zero_padding() const; 
 	bool is_zero() const {zero_padding(); for(int i=0; i<nn; ++i) if(v[i]) return false; return true;};
 	bool is_one() const {zero_padding(); if(v[0]!=1) return false; for(int i=1; i<nn; ++i) if(v[i]) return false;return true;};
@@ -251,5 +265,17 @@ public:
 	unsigned int population(const unsigned int before=0) const {return bitvector::population(before?before-1:0);};
 	};
 
+
+//some necessary traits of the non-scalar class to be able to use LA methods
+template<>
+class LA_traits<bitvector> {
+public: 
+        static bool is_plaindata() {return false;};
+        static void copyonwrite(bitvector& x) {x.copyonwrite();};
+	typedef bool elementtype; 
+	typedef bool normtype; 
+};
+
+
 }//namespace
 #endif

conjgrad.h

@@ -67,9 +67,9 @@ for(int iter=0; iter<= itmax; iter++)
 	double err=p.norm();
 	if(verbose) 
 		{
-		std::cout << "conjgrad: iter= "<<iter<<" err= "<<
+		std::cout << "conjgrad: iter= "<<iter<<" err= "<< err<<std::endl;
-		std::setiosflags(std::ios::scientific)<<std::setprecision(8) <<err<<
+		//std::setiosflags(std::ios::scientific)<<std::setprecision(8) <<err<<
-		std::resetiosflags(std::ios::scientific)<<std::setprecision(12)<<"\n";
+		//std::resetiosflags(std::ios::scientific)<<std::setprecision(12)<<"\n";
 		std::cout.flush();
 		}
 	if(err <= tol)

davidson.h

@@ -46,7 +46,7 @@ namespace LA {
 
 
 template <typename T, typename Matrix>
-extern void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile, 
+void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile, 
 		int nroots=1,  const bool verbose=0, const double eps=1e-6,
 	 	const bool incore=1, int maxit=100, const int maxkrylov = 500,
 		void (*initguess)(NRVec<T> &)=NULL, const typename LA_traits<T>::normtype *target=NULL)
@@ -309,6 +309,25 @@ if(incore) {delete[] v0; delete[] v1;}
 else  {delete s0; delete s1;}
 
 if(flag) laerror("no convergence in davidson");
+} //davidson
+
+
+//reconstruction of explicit dense matrix from the implicit one (useful for debugging)
+template <typename T, typename Matrix>
+NRMat<T> explicit_matrix(const Matrix &bigmat) 
+{
+NRMat<T> r(bigmat.nrows(), bigmat.ncols());
+for(int i=0; i<bigmat.ncols(); ++i)
+        {
+	NRVec<T> ket(bigmat.ncols());
+        ket.clear();
+        ket[i]=(T)1;
+        NRVec<T> hket(bigmat.nrows());
+        bigmat.gemv((T)0,hket,'n',(T)1,ket);
+	for(int l=0; l<bigmat.nrows(); ++l) r(l,i) = hket[l];
+        }
+
+return r;
 }
 
 }//namespace

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;}\

laerror.cc

@@ -109,8 +109,8 @@ void laerror2(const char *s1, const char *s2)
   std::cerr.flush();
   if(s1)
   {
-    std::cerr <<"LA:ERROR - "<< s2 << ": " << s1 << "\n";
+    std::cerr <<"LA:ERROR - "<< s1 << "\nCalled from " << s2 << "\n";
-    std::cout <<"LA:ERROR - "<< s2 << ": " << s1 << "\n";
+    std::cout <<"LA:ERROR - "<< s1 << "\nCalled from " << s2 << "\n";
   }
   else 
 	{

lanczos.h

@@ -50,7 +50,7 @@ if ( n!= (int)bigmat.ncols()) laerror("non-square matrix in lanczos");
 if(eivals.size()<nroots) laerror("too small eivals dimension in lanczos");
 
 NRVec<T> vec1(n),vec2(n);
-NRVec<T> *v0,*v1;
+NRVec<T> *v0;
 AuxStorage<T> *s0,*s1;
 
 if(incore)
@@ -112,9 +112,27 @@ for(j=1; j<maxkrylov;++j)
 		}
 	else
 		{
-		laerror("zero norm in lanczos");
+		//generate an arbitrary vector and orthonormalize it to all previous v_j
-		//could generate an arbitrary vector and orthonormalize it
+		vec2.randomize(1.);
+		for(int k=0; k<j; ++k)
+			{
+			T f;
+			if(incore) 
+				{
+				f = v0[k].dot(vec2);
+				vec2.axpy(-f,v0[k]);
+				}
+			else
+				{
+				NRVec<T> vec3(n);
+				s0->get(vec3,k);
+                                f = vec3.dot(vec2);
+                                vec2.axpy(-f,vec3);
+				}
+			}
+		vec2.normalize();
 		}
+	//vec2 now stores v_j
 	if(incore) v0[j]=vec2; else s0->put(vec2,j);
 	vec1 *= -beta[j-1];
 	bigmat.gemv(1,vec1,'n',1,vec2);

mat.cc

@@ -97,22 +97,47 @@ const NRMat<T> NRMat<T>::otimes(const NRMat<T> &rhs, bool reversecolumns) const
  * @return extracted elements as a NRVec<T> object
  ******************************************************************************/
 template <typename T>
-const NRVec<T> NRMat<T>::row(const int i, int l) const {
+const NRVec<T> NRMat<T>::row(const int i, int l, int offset) const {
+	if(l < 0) l = mm;
 #ifdef DEBUG
 	if(i < 0 || i >= nn) laerror("illegal index");
+	if(offset<0||l+offset>mm) laerror("illegal len/offset");
 #endif
-	if(l < 0) l = mm;
 	NRVec<T> r(l);
 	LA_traits<T>::copy(&r[0],
 #ifdef MATPTR
-	        v[i]
+	        v[i]+offset
 #else
-		v + i*(size_t)l
+		v + i*(size_t)mm + offset
 #endif
 		, l);
 	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, int offset) {
+	if(l < 0) l = mm;
+#ifdef DEBUG
+	if(i < 0 || i >= nn) laerror("illegal index");
+	if(offset<0||l+offset>mm) laerror("illegal len/offset");
+#endif
+	LA_traits<T>::copy(
+#ifdef MATPTR
+	        v[i]+offset
+#else
+		v + i*(size_t)mm  + offset
+#endif
+		, &r[0], l);
+}
+
+
+
 /***************************************************************************//**
  * routine for raw output
  * @param[in] fd file descriptor for output

mat.h

@@ -35,6 +35,9 @@ template<typename T, typename R> class WeightPermutation;
 template<typename T, typename R> class PermutationAlgebra;
 template<typename T> class CyclePerm;
 template<typename T> class NRMat_from1;
+template<typename T> class SparseMat;
+template<typename T> class SparseSMat;
+template<typename T> class SemiSparseMat;
 
 
 /***************************************************************************//**
@@ -264,7 +267,10 @@ public:
 	void orthonormalize(const bool rowcol, const NRSMat<T> *metric = NULL);
 
 	//! get the i<sup>th</sup> row
-	const NRVec<T> row(const int i, int l = -1) const;
+	const NRVec<T> row(const int i, int len = -1, int offset=0) const;
+
+	//! set the i<sup>th</sup> row
+	void rowset(const NRVec<T> &r, const int i, int len = -1, int offset=0);
 
 	//! get the j<sup>th</sup> column
 	const NRVec<T> column(const int j, int l = -1) const {
@@ -275,6 +281,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
@@ -394,6 +407,8 @@ public:
 	explicit NRMat(const SparseSMat<T> &rhs);
 	//! explicit constructor converting sparse CSR matrix into \c NRMat<T> object
         explicit NRMat(const CSRMat<T> &rhs);
+	//! explicit constructor converting sparse matrix into \c NRMat<T> object
+	explicit NRMat(const SemiSparseMat<T> &rhs);
 
 	//! add up given sparse matrix
 	NRMat & operator+=(const SparseMat<T> &rhs);
@@ -743,6 +758,7 @@ inline const NRMat<T> NRMat<T>::operator-(const NRSMat<T> &rhs) const {
 template <typename T>
 inline T* NRMat<T>::operator[](const int i) {
 #ifdef DEBUG
+	if(nn==0||mm==0) laerror("operator[] on unallocated matrix");
 	if (_LA_count_check && *count != 1) laerror("matrix with *count>1 used as l-value");
 	if(i<0||i>=nn) std::cout<<"INDEX PROBLEM "<<0<<" "<<i<<" "<<nn-1<<std::endl;
 	if (i < 0) laerror("Mat [] out of range - low");
@@ -786,6 +802,7 @@ inline const T* NRMat<T>::operator[](const int i) const {
 template <typename T>
 inline T& NRMat<T>::operator()(const int i, const int j){
 #ifdef DEBUG
+	if(nn==0||mm==0) laerror("operator() of unallocated matrix");
 	if (_LA_count_check && *count != 1) laerror("NRMat::operator(,) used as l-value for a matrix with count > 1");
 	if(i<0||i>=nn||j<0||j>mm) std::cout<<"INDEX PROBLEM "<<0<<" "<<i<<" "<<nn-1<<" "<<0<<" "<<j<<" "<<mm-1<<std::endl;
 	if (i < 0) laerror("first index out of range - low");

noncblas.cc

@@ -648,6 +648,7 @@ int clapack_dgesv(const CBLAS_ORDER Order, const int N, const int NRHS,
 		double *A, const int lda, int *ipiv,
 		double *B, const int ldb)
 {
+//std::cout <<"In MY clapack_dgesv, N and NRHS = "<< N<<" "<< NRHS<<"\n";
 	FINT INFO=0;
 	if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted");
 	//B should be in the same physical order, just transpose A in place and the LU result on output
@@ -657,8 +658,9 @@ int clapack_dgesv(const CBLAS_ORDER Order, const int N, const int NRHS,
 	const FINT nrhstmp=NRHS;
 	const FINT ldatmp=lda;
 	const FINT ldbtmp=ldb;
-	FINT ipivtmp=*ipiv;
+	FINT ipivtmp[N];
-	FORNAME(dgesv) (&ntmp,&nrhstmp,A,&ldatmp,&ipivtmp,B,&ldbtmp,&INFO);
+	FORNAME(dgesv) (&ntmp,&nrhstmp,A,&ldatmp,ipivtmp,B,&ldbtmp,&INFO);
+	for(int i=0; i<N; ++i) ipiv[i]=ipivtmp[i];
 #else
 	FORNAME(dgesv) (&N,&NRHS,A,&lda,ipiv,B,&ldb,&INFO);
 #endif
@@ -672,6 +674,7 @@ int clapack_sgesv(const CBLAS_ORDER Order, const int N, const int NRHS,
 		float *A, const int lda, int *ipiv,
 		float *B, const int ldb)
 {
+std::cout <<"In my clapack_sgesv\n";
 	FINT INFO=0;
 	if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted");
 	//B should be in the same physical order, just transpose A in place and the LU result on output
@@ -681,8 +684,9 @@ int clapack_sgesv(const CBLAS_ORDER Order, const int N, const int NRHS,
 	const FINT nrhstmp=NRHS;
 	const FINT ldatmp=lda;
 	const FINT ldbtmp=ldb;
-	FINT ipivtmp=*ipiv;
+	FINT ipivtmp[N];
-	FORNAME(sgesv) (&ntmp,&nrhstmp,A,&ldatmp,&ipivtmp,B,&ldbtmp,&INFO);
+	FORNAME(sgesv) (&ntmp,&nrhstmp,A,&ldatmp,ipivtmp,B,&ldbtmp,&INFO);
+	for(int i=0; i<N; ++i) ipiv[i]=ipivtmp[i];
 #else
 	FORNAME(sgesv) (&N,&NRHS,A,&lda,ipiv,B,&ldb,&INFO);
 #endif

nonclass.cc

@@ -26,7 +26,7 @@
 #include "qsort.h"
 #include "fortran.h"
 
-#undef IPIV_DEBUG
+//#define IPIV_DEBUG
 
 
 namespace LA {
@@ -140,11 +140,23 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
 	int r, *ipiv;
 	int iswap=0;
 
-	
+	if(nrhs==0) std::cout<<"Warning: some dgesv implementations might skip LU decomposition when nrhs==0\n";	
 	if (n==A.nrows() && A.nrows() != A.ncols()) laerror("linear_solve() call for non-square matrix");
 	A.copyonwrite();
 	ipiv = new int[A.nrows()];
-	r = clapack_dgesv(CblasRowMajor, n, nrhs, A[0], A.ncols(), ipiv, B , ldb);
+#ifdef IPIV_DEBUG
+	for(int i=0; i<A.nrows(); ++i) ipiv[i]=123456789;
+	std::cout <<"canary ipiv initialized\n";
+	std::cout <<"A before clapack_dgesv = "<<A<<std::endl;
+	std::cout <<"Active dimension n= "<<n<<std::endl;
+#endif
+	r = clapack_dgesv(CblasRowMajor, n, nrhs, &A(0,0), A.ncols(), ipiv, B , ldb);
+#ifdef IPIV_DEBUG
+	std::cout <<"A after clapack_dgesv = "<<A<<std::endl;
+	std::cout <<"ipiv = ";
+	for (int i=0; i<n; ++i) std::cout <<ipiv[i]<<" ";
+	std::cout <<std::endl;
+#endif
 	if (r < 0) {
 		delete[] ipiv;
 		laerror("illegal argument in lapack_gesv");
@@ -158,7 +170,11 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
 	        for (int i=0; i<n; ++i)
 			{
 			if(ipiv[i]==0) shift=0;
-			if(ipiv[i]<0 || ipiv[i]>n) laerror("problem with ipiv in clapack_dgesv");
+			if(ipiv[i]<0 || ipiv[i]>n) 
+				{
+				std::cout <<"IPIV["<<i<<"] = "<<ipiv[i]<<std::endl;
+				laerror("problem with ipiv in clapack_dgesv");
+				}
 			}
 #ifdef IPIV_DEBUG
        std::cout <<"shift = "<<shift<<std::endl;
@@ -169,10 +185,6 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
 	if(det && r>0) *det = 0;
 #ifdef IPIV_DEBUG
        std::cout <<"iswap = "<<iswap<<std::endl;
-
-	std::cout <<"ipiv = ";
-	for (int i=0; i<n; ++i) std::cout <<ipiv[i]<<" ";
-	std::cout <<std::endl;
 #endif
 
 	delete [] ipiv;
@@ -204,6 +216,7 @@ extern "C" void FORNAME(dspsv)(const char *UPLO, const FINT *N, const FINT *NRHS
 
 static void linear_solve_do(NRSMat<double> &a, double *b, const int nrhs, const int ldb, double *det, int n)
 {
+ if(nrhs==0) std::cout<<"Warning: some dspsv implementations might skip LU decomposition when nrhs==0\n";
 	FINT r, *ipiv;
 	a.copyonwrite();
 	ipiv = new FINT[n];
@@ -216,13 +229,14 @@ static void linear_solve_do(NRSMat<double> &a, double *b, const int nrhs, const
 #else	
 	FORNAME(dspsv)(&U, &n, &nrhs, a, ipiv, b, &ldb,&r);
 #endif	
+//	std::cout <<"A after dspsv = "<<a<<std::endl;
 	if (r < 0) {
 		delete[] ipiv;
 		laerror("illegal argument in spsv() call of linear_solve()");
 	}
 	if (det && r == 0) {
 		*det = 1.;
-		for (int i=1; i<n; i++) {double t=a(i,i); if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *= t;}
+		for (int i=0; i<n; i++) {double t=a(i,i); if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *= t;}
 		//do not use ipiv, since the permutation matrix occurs twice in the decomposition and signs thus cancel (man dspsv)
 	}
         if (det && r>0)  *det = 0;
@@ -282,7 +296,11 @@ void linear_solve(NRMat< std::complex<double> > &A, NRMat< std::complex<double>
                 for (int i=0; i<n; ++i)
                         {
                         if(ipiv[i]==0) shift=0;
-                        if(ipiv[i]<0 || ipiv[i]>n) laerror("problem with ipiv in zgesv");
+                        if(ipiv[i]<0 || ipiv[i]>n) 
+                                {
+                                std::cout <<"IPIV["<<i<<"] = "<<ipiv[i]<<std::endl;
+                                laerror("problem with ipiv in zgesv");
+                                }
                         }
 #ifdef IPIV_DEBUG
        std::cout <<"shift = "<<shift<<std::endl;
@@ -918,9 +936,6 @@ void singular_decomposition(NRMat<std::complex<double> > &a, NRMat<std::complex<
 }
 
 
-
-
-
 //QR decomposition
 //extern "C" void FORNAME(dgeqrf)(const int *M, const int *N, double *A, const int *LDA, double *TAU, double *WORK, int *LWORK, int *INFO);
 

nonclass.h

@@ -158,13 +158,60 @@ extern void linear_solve(NRMat<T> &a, NRVec<T> &b, double *det=0, int n=0); \
 extern void linear_solve(NRSMat<T> &a, NRVec<T> &b, double *det=0, int n=0); \
 extern void diagonalize(NRMat<T> &a, NRVec<LA_traits<T>::normtype> &w, const bool eivec=1, const bool corder=1, int n=0, NRMat<T> *b=NULL, const int itype=1); \
 extern void diagonalize(NRSMat<T> &a, NRVec<LA_traits<T>::normtype> &w, NRMat<T> *v, const bool corder=1, int n=0, NRSMat<T> *b=NULL, const int itype=1);\
-extern void singular_decomposition(NRMat<T> &a, NRMat<T> *u, NRVec<LA_traits<T>::normtype> &s, NRMat<T> *v, const bool vnotdagger=0, int m=0, int n=0);
+extern void singular_decomposition(NRMat<T> &a, NRMat<T> *u, NRVec<LA_traits<T>::normtype> &s, NRMat<T> *v, const bool vnotdagger=0, int m=0, int n=0); \
 
 /*NOTE!!! all versions of diagonalize DESTROY A and generalized diagonalize also B matrix */
 
 declare_la(double)
 declare_la(std::complex<double>)
 
+
+
+//svd_solve without contamination from nullspace, b is n x nrhs, result returns in b, A is destroyed
+template<typename T>
+NRVec<T> svd_solve(NRMat<T> &a, const NRVec<T> &b, double thres=1e-12, int *nullity = NULL)
+{
+if(b.size()!=a.nrows()) laerror("size mismatch in svd_solve");
+if(nullity) *nullity=0;
+NRVec<double> w(a.ncols());
+NRMat<T> u(a.nrows(),a.ncols());
+NRMat<T> vt(a.ncols(),a.ncols());
+singular_decomposition(a,&u,w,&vt,false);
+NRVec<T> utb = b*u;
+for(int i=0; i<w.size(); ++i) 
+	{
+	if(w[i]>thres) utb[i] /= w[i];
+	else {utb[i]=0; if(nullity) ++*nullity;}
+	}
+return utb*vt;
+}
+
+
+template<typename T>
+NRMat<T> svd_solve(NRMat<T> &a, const NRMat<T> &b, double thres=1e-12, int *nullity = NULL)
+{
+if(b.nrows()!=a.nrows()) laerror("size mismatch in svd_solve");
+if(nullity) *nullity=0;
+NRVec<double> w(a.ncols());
+NRMat<T> u(a.nrows(),a.ncols());
+NRMat<T> vt(a.ncols(),a.ncols());
+singular_decomposition(a,&u,w,&vt,false);
+NRVec<T> ww(w.size());
+for(int i=0; i<w.size(); ++i)
+	{
+	if(w[i]>thres)   ww[i] =  1./w[i];
+	else {ww[i]=0; if(nullity) ++*nullity;}
+	}
+NRMat<T> utb(a.ncols(),b.ncols());
+utb.gemm((T)0,u,'c',b,'n',(T)1);
+utb.diagmultl(ww);
+NRMat<T> res(a.ncols(),b.ncols());
+res.gemm((T)0,vt,'c',utb,'n',(T)1);
+return res;
+}
+
+
+
 // Separate declarations
 //general nonsymmetric matrix and generalized diagonalization
 //corder =0 ... C rows are eigenvectors, =1 ... C columns are eigenvectors
@@ -225,16 +272,9 @@ typename LA_traits<MAT>::normtype MatrixNorm(const MAT  &A, const char norm);
 template<class MAT>
 typename LA_traits<MAT>::normtype CondNumber(const MAT  &A, const char norm);
 
-
+#ifdef HAS_MKL
-//general determinant
+#define NO_OPENBLAS_WORKAROUND
-template<class MAT>
+#endif
-const typename LA_traits<MAT>::elementtype determinant(MAT a)//passed by value
-{
-typename LA_traits<MAT>::elementtype det;
-if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix");
-linear_solve(a,NULL,&det);
-return det;
-}
 
 //general determinant destructive on input 
 template<class MAT>
@@ -242,10 +282,27 @@ const typename LA_traits<MAT>::elementtype determinant_destroy(MAT &a) //passed
 {
 typename LA_traits<MAT>::elementtype det;
 if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix");
+
+//for openblas 0.3.31 we have to fake some RHS otherwise LU decomp. is not performed
+#ifdef NO_OPENBLAS_WORKAROUND
 linear_solve(a,NULL,&det);
+#else
+//fake rhs
+NRVec<typename LA_traits<MAT>::elementtype> b(a.ncols());
+for(int i=0; i<b.size(); ++i) b[i] = (typename LA_traits<MAT>::elementtype)i;
+linear_solve(a,b,&det);
+#endif
 return det;
 }
 
+//general determinant
+template<class MAT>
+const typename LA_traits<MAT>::elementtype determinant(MAT a)//passed by value
+{
+a.copyonwrite();
+return determinant_destroy(a);
+}
+
 
 //------------------------------------------------------------------------------
 // solves set of linear equations using gesvx

smat.h

@@ -770,6 +770,7 @@ return (i>j)? (i-2)*(i-1)/2+j-1 : (j-2)*(j-1)/2+i-1;
 template <typename T>
 inline T & NRSMat<T>::operator()(const int i, const int j) {
 #ifdef DEBUG
+	if(nn==0) laerror("operator() of unallocated smatrix");
 	if(_LA_count_check && *count != 1) laerror("T & NRSMat<T>::operator()(const int, const int) used for matrix with count > 1");
 	if(i<0||i>=nn||j<0||j>=nn) std::cout<<"INDEX PROBLEM "<<0<<" "<<i<<" "<<nn-1<<" "<<0<<" "<<j<<" "<<nn-1<<std::endl;
 	if(i<0) laerror("T & NRSMat<T>::operator()(const int, const int) first index out of range - low");

sparsemat.cc

@@ -557,6 +557,16 @@ for(i=0;i<nn;++i)
 		}
 }
 
+template <class T>
+NRMat<T>::NRMat(const SemiSparseMat<T> &rhs)
+		: NRMat(rhs.offdiagonal)
+		{
+		if(rhs.diagonal.size()>0)
+			{
+			for(int i=0; i<nn; ++i) (*this)(i,i) += rhs.diagonal[i];
+			}
+		}
+
 
 //get diagonal, do not construct a new object, but store in existing one, important for huge CI matrices
 // with the divide option is used as a preconditioner, another choice of preconditioner is possible
@@ -1397,8 +1407,12 @@ INSTANTIZE(std::complex<double>) //some functions are not OK for hermitean matri
 //// forced instantization in the corresponding object file
 template class SparseMat<double>;
 template class SparseMat<std::complex<double> >;
+template class SemiSparseMat<double>;
+template class SemiSparseMat<std::complex<double> >;
+
 
 #define INSTANTIZE(T) \
+template NRMat<T>::NRMat(const SemiSparseMat<T> &rhs); \
 template NRMat<T>::NRMat(const SparseMat<T> &rhs); \
 template NRSMat<T>::NRSMat(const SparseMat<T> &rhs); \
 template NRVec<T>::NRVec(const SparseMat<T> &rhs); \

sparsemat.h

@@ -321,5 +321,54 @@ while(l)
 return *this;
 }
 
+
+//a rudimentary class for sparse matrix with a dense diagonal - implemented just enough to call davidson etc. on it
+template <typename T>
+class SemiSparseMat {
+	friend class NRMat<T>;
+public:
+	NRVec<T> diagonal;
+	SparseMat<T> offdiagonal;
+        SemiSparseMat() {};
+        SemiSparseMat(const SPMatindex n, const SPMatindex m) : offdiagonal(n,m) {if(n==m) diagonal.resize(n);};
+        SPMatindex nrows() const {return offdiagonal.nrows();}
+        SPMatindex ncols() const {return offdiagonal.ncols();}
+        SPMatindex issymmetric() const {return offdiagonal.issymmetric();}
+	void setsymmetric() {offdiagonal.setsymmetric();}
+	void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x, bool treat_as_symmetric=false) const
+		{
+		offdiagonal.gemv(beta,r,trans,alpha,x,treat_as_symmetric);
+		if(diagonal.size()>0 && alpha!=(T)0)
+			{
+			if(diagonal.size()!=r.size()) laerror("mismatch in diagonal size");
+			for(int i=0; i<diagonal.size(); ++i) r[i] += alpha * diagonal[i]*x[i];
+			}
+		}
+	NRVec<T> operator*(const NRVec<T> &rhs) const {NRVec<T> result(nrows()); gemv((T)0,result,'n',(T)1,rhs); return result;};
+	const T* diagonalof(NRVec<T> &x, const bool divide=0, bool cache=false) const
+		{
+		if(diagonal.size()>0) //we ASSUME the diagonal is ONLY in the vector
+			{
+			if(diagonal.size()!=x.size()) laerror("mismatch in diagonal size");
+			if(divide) {x.copyonwrite(); for (int i=1; i<x.size(); ++i) x[i]/=diagonal[i]; return NULL;}
+			else	{x|=diagonal; return &x[0];}
+			}
+		return offdiagonal.diagonalof(x,divide,cache);
+		}
+	void add(const SPMatindex n, const SPMatindex m, const T elem) {if(diagonal.size()>0 && n==m) diagonal[n]+=elem; else offdiagonal.add(n,m,elem);}
+	void clear() {diagonal.clear(); offdiagonal.clear();}
+	void copyonwrite() {diagonal.copyonwrite(); offdiagonal.copyonwrite();}
+	void resize(const SPMatindex n, const SPMatindex m) {offdiagonal.resize(n,m); diagonal.resize(n==m?n:0);}
+};
+
+template <class T>
+std::ostream& operator<<(std::ostream &s, const SemiSparseMat<T> &x)
+                {
+		s << "diagonal= "<<x.diagonal;
+		s << "offdiagonal= "<<x.offdiagonal;
+                return s;
+                }
+
+
 }//namespace
 #endif

t.cc

@@ -61,6 +61,11 @@ inline int randind(const int n)
 
 complex<double> mycident (const complex<double>&x) {return x;}
 
+void randomguess(NRVec<double> &v)
+{
+v.randomize(1);
+}
+
 void printme(const NRPerm<int> &p)
 {
 PERM_RANK_TYPE rank=p.rank();
@@ -1079,13 +1084,13 @@ NRMat<complex<double> > b=exp(a);
 cout <<b;
 }
 
-if(0)
+if(1)
 {
 int n;
 double d;
 cin >>n;
-//NRMat<double> a(n,n);
+NRMat<double> a(n,n);
-NRSMat<double> a(n);
+//NRSMat<double> a(n);
 for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
         {
         a(j,i)=a(i,j)=RANDDOUBLE()*(i==j?10.:1.);
@@ -4572,7 +4577,7 @@ for(int i=0; i<m; ++i)
 
 }
 
-if(1)
+if(0)
 {
 int r,n,sym;
 cin>>r>>n>>sym;
@@ -4597,6 +4602,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;
@@ -4622,4 +4629,122 @@ cout <<z.shape;
 
 }
 
+#undef sparsity
+#define sparsity (n*2)
+if(0)
+{
+int n,m;
+cin >>n>>m;
+	SemiSparseMat<double> aa(n,n);
+	aa.setsymmetric();
+	for(int i=0; i<sparsity;i++) aa.add(randind(n),randind(n),RANDDOUBLE());
+	for(int i=0; i<n; ++i) aa.add(i,i,500*RANDDOUBLE());
+NRVec<double> r(m);
+davidson(aa,r,(NRVec<double> *)NULL,"eivecs",m,1,1e-5,0,300,300);
+cout <<r;
+if(n<=20)
+	{
+	cout<<aa;
+	NRMat a(aa);
+	NRVec<double> b(n);
+	cout <<a;
+	diagonalize(a,b);
+	cout <<a;
+	cout <<b;
+
+	}
+}
+
+if(0)
+{
+int n,m;
+cin>> n>>m;
+NRMat<double> a(n,m);
+a.randomize(1.);
+NRMat<double> b = explicit_matrix<double,NRMat<double> >(a);
+cout <<"Error = "<<(a-b).norm()<<endl;
+}
+
+if(0)
+{
+int m;
+int which;
+cin>>m >>which;
+cout <<"here\n";
+NRSMat<double> a;
+cin>>a;
+int n=a.nrows();
+NRVec<double> rr(n);
+
+NRSMat<double> aa(a);
+NRMat<double> vv(n,n);
+diagonalize(aa,rr,&vv);
+NRVec<double> r(m);
+NRVec<double> *eivecs = new NRVec<double>[m];
+cout <<"Exact energies " <<rr<<endl;
+
+if(which) lanczos(a,r,eivecs,NULL,m,1,1e-6,1,500,500,randomguess);
+else davidson(a,r,eivecs,NULL,m,1,1e-6,1,500,500,randomguess);
+cout <<"Iterative energies " <<r;
+cout <<"Eigenvectors compare:\n";
+for(int i=0; i<m; ++i)
+        {
+        cout <<eivecs[i];
+        for(int j=0; j<n;++j) cout <<vv[j][i]<<" ";
+        cout <<endl;
+        }
+
+}
+
+if(0)
+{
+NRMat<double> a;
+NRVec<double> b;
+cin >>a>>b;
+NRMat<double> aa(a);
+NRVec<double> x = svd_solve(aa,b);
+cout <<x;
+cout <<"Error = "<< (a*x-b).norm()<<endl;
+}
+
+if(0)
+{
+NRMat<double> a;
+NRMat<double> b;
+cin >>a>>b;
+NRMat<double> aa(a);
+NRMat<double> x = svd_solve(aa,b);
+cout <<x;
+cout <<"Error = "<< (a*x-b).norm()<<endl;
+}
+
+if(0)
+{
+int n;
+cin>>n;
+NRMat<double> a(n,n);
+a.randomize(.5);
+NRMat<double>  ae = exp(a);
+NRMat<double> aa = log(ae);
+cout <<a;
+cout <<aa;
+cout <<"Error = "<<(a-aa).norm()<<endl;
+}
+
+
+if(0)
+{
+int n;
+cin>>n;
+NRMat<double> a(n,n);
+a.randomize(.5);
+NRMat<double>  ax = a.transpose()*a;
+NRMat<double> al = log(ax);
+NRMat<double> aa = exp(al);
+cout <<ax;
+cout <<aa;
+cout <<"Error = "<<(ax-aa).norm()<<endl;
+}
+
+
 }//main

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,31 @@ 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());
+if(rank()>0)
+	{
+	indexmat_p = &r;
+	indexmat_row = 0;
+	constloopover(indexmatrix_callback<T>,false);
+	}
+return r;
 }
 
 
@@ -606,7 +634,7 @@ std::ostream & operator<<(std::ostream &s, const Tensor<T> &x)
 {
 s<<x.shape;
 s<<x.names;
-if(x.rank()==0) {s<<x.data[0]; return s;}
+if(x.rank()==0) {s<<x.data[0]<<std::endl; return s;}
 sout= &s;
 x.constloopover(&outputcallback<T>);
 return s;
@@ -2439,6 +2467,34 @@ 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;
+}
+
+bool  shares_index(const FLATINDEX &I, const FLATINDEX &J)
+{
+for(int i=0; i<I.size(); ++i) for(int j=0; j<J.size(); ++j)
+	if(I[i]==J[j]) 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,17 @@ 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
 
+LA_largeindex subindex(int *sign, const INDEXGROUP &g, const NRVec<LA_index> &I); //index of one subgroup
+NRVec<LA_index>  inverse_subindex(const INDEXGROUP &g, LA_largeindex s);
+
+
+//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);
+
+bool shares_index(const FLATINDEX &I, const FLATINDEX &J);
+
 FLATINDEX superindex2flat(const SUPERINDEX &I);
 
 template<typename T>
@@ -361,11 +373,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)
 

vec.h

@@ -1101,6 +1101,7 @@ inline const T NRVec<T>::dot(const T *a, const int stride , bool conjugate) cons
 template <typename T>
 inline T& NRVec<T>::operator[](const int i) {
 #ifdef DEBUG
+	if(nn==0) laerror("operator[] of unallocated vector");
 	if(_LA_count_check && *count != 1) laerror("possible use of NRVec[] with count>1 as l-value");
 	if(i<0||i >= nn) std::cout<<"INDEX PROBLEM "<<0<<" "<<i<<" "<<nn-1<<std::endl;
 	if(i < 0) laerror("out of range - low");
@@ -1436,7 +1437,7 @@ return;
 
 
 /***************************************************************************//**
- * perfrom deep copy
+ * perform a deep copy
  * @param[in] rhs vector being copied
  * @see NRVec<T>::copyonwrite()
  ******************************************************************************/