SemiSparseMat introduced

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;
 
 
 /***************************************************************************//**
@@ -404,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);

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,46 @@ 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>;
+protected:
+	NRVec<T> diagonal;
+	SparseMat<T> offdiagonal;
+public:
+        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];
+			}
+		}
+	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);}
+};
+
+
 }//namespace
 #endif

t.cc

@@ -4572,7 +4572,7 @@ for(int i=0; i<m; ++i)
 
 }
 
-if(1)
+if(0)
 {
 int r,n,sym;
 cin>>r>>n>>sym;
@@ -4624,4 +4624,30 @@ cout <<z.shape;
 
 }
 
+#undef sparsity
+#define sparsity (n*2)
+if(1)
+{
+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)
+	{
+	NRMat a(aa);
+	NRVec<double> b(n);
+	cout <<a;
+	diagonalize(a,b);
+	cout <<a;
+	cout <<b;
+
+	}
+}
+
+
 }//main

vec.h

@@ -1436,7 +1436,7 @@ return;
 
 
 /***************************************************************************//**
- * perfrom deep copy
+ * perform a deep copy
  * @param[in] rhs vector being copied
  * @see NRVec<T>::copyonwrite()
  ******************************************************************************/