From b1f1be84579ecd10eaabdc1c29a1af2c729e42d8 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Tue, 3 Feb 2026 15:46:39 +0100
Subject: [PATCH] SemiSparseMat introduced

---
 mat.h        |  5 +++++
 sparsemat.cc | 14 ++++++++++++++
 sparsemat.h  | 41 +++++++++++++++++++++++++++++++++++++++++
 t.cc         | 28 +++++++++++++++++++++++++++-
 vec.h        |  2 +-
 5 files changed, 88 insertions(+), 2 deletions(-)

diff --git a/mat.h b/mat.h
index abf2f13..26c868f 100644
--- a/mat.h
+++ b/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);
diff --git a/sparsemat.cc b/sparsemat.cc
index 586014d..d300882 100644
--- a/sparsemat.cc
+++ b/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); \
diff --git a/sparsemat.h b/sparsemat.h
index dcfd509..3542a0b 100644
--- a/sparsemat.h
+++ b/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
diff --git a/t.cc b/t.cc
index dd1612c..d72edd6 100644
--- a/t.cc
+++ b/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
diff --git a/vec.h b/vec.h
index 110927e..08c2d2c 100644
--- a/vec.h
+++ b/vec.h
@@ -1436,7 +1436,7 @@ return;
 
 
 /***************************************************************************//**
- * perfrom deep copy
+ * perform a deep copy
  * @param[in] rhs vector being copied
  * @see NRVec<T>::copyonwrite()
  ******************************************************************************/