From 51a26337c6f28f73766f6868c3b3eedaaf669171 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Wed, 26 Jul 2023 21:18:57 +0200
Subject: [PATCH] implemented subvector/submatrix with individual index
 selection

---
 mat.cc  |  45 +++++++++++++++++++++++
 mat.h   |   2 +
 smat.cc | 112 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 smat.h  |   9 +++++
 vec.cc  |  34 +++++++++++++++++
 vec.h   |   2 +
 6 files changed, 204 insertions(+)

diff --git a/mat.cc b/mat.cc
index f7feb29..314d3eb 100644
--- a/mat.cc
+++ b/mat.cc
@@ -763,6 +763,30 @@ const NRMat<T> NRMat<T>::submatrix(const int fromrow, const int torow, const int
 	return r;
 }
 
+
+template <typename T>
+const NRMat<T> NRMat<T>::submatrix(const NRVec<int> &rows, const NRVec<int> &cols) const
+{
+	NOT_GPU(*this);
+        const int n = rows.size();
+        const int m = cols.size();
+        NRMat<T> r(n,m);
+
+        for(int i=0; i<n; ++i)
+		{
+		int ii=rows[i];
+		if(ii<0||ii>=nn) laerror("bad row index in submatrix");
+		for(int j=0; j<m; ++j)
+			{
+			int jj=cols[j];
+			if(jj<0||jj>=mm) laerror("bad col index in submatrix");
+			r(i,j) = (*this)(ii,jj);
+			}
+                }
+return r;
+}
+
+
 /***************************************************************************//**
  * places given matrix as submatrix at given position
  * @param[in] fromrow row-coordinate of top left corner
@@ -799,6 +823,27 @@ void NRMat<T>::storesubmatrix(const int fromrow, const int fromcol, const NRMat
 	}
 }
 
+template <typename T>
+void NRMat<T>::storesubmatrix(const NRVec<int> &rows, const NRVec<int> &cols, const NRMat &rhs)
+{
+        NOT_GPU(*this);
+        const int n = rows.size();
+        const int m = cols.size();
+	if(rhs.nrows()!=n || rhs.ncols()!=m) laerror("incompatible dimensions in storesubmatrix");
+                
+        for(int i=0; i<n; ++i)
+                {
+                int ii=rows[i];
+                if(ii<0||ii>=nn) laerror("bad row index in storesubmatrix");
+                for(int j=0; j<m; ++j)
+                        {
+                        int jj=cols[j];
+                        if(jj<0||jj>=mm) laerror("bad col index in storesubmatrix");
+                        (*this)(ii,jj) = rhs(i,j);
+                        }
+                }       
+}
+
 /***************************************************************************//**
  * compute matrix transposition for a principal leading minor
  * @param[in] _n order of the leading minor
diff --git a/mat.h b/mat.h
index 5a5a968..a97861a 100644
--- a/mat.h
+++ b/mat.h
@@ -326,9 +326,11 @@ public:
 
 	//! extract specified submatrix
         const NRMat submatrix(const int fromrow, const int torow, const int fromcol, const int tocol) const;
+	const NRMat submatrix(const NRVec<int> &rows, const NRVec<int> &cols) const;
 
 	//! store given matrix at given position into the current matrix
 	void storesubmatrix(const int fromrow, const int fromcol, const NRMat &rhs);
+	void storesubmatrix(const NRVec<int> &rows, const NRVec<int> &cols, const NRMat &rhs);
 
 	//! perform the \b gemm operation
 	void gemm(const T &beta, const NRMat &a, const char transa, const NRMat &b, const char transb, const T &alpha);
diff --git a/smat.cc b/smat.cc
index 8d7c3b1..9939a62 100644
--- a/smat.cc
+++ b/smat.cc
@@ -727,6 +727,118 @@ NRSMat<std::complex<double> >::NRSMat(const NRSMat<double> &rhs, bool imagpart):
 #endif
 }
 
+
+
+/***************************************************************************//**
+ * extract block submatrix
+ * @param[in] from starting position
+ * @param[in] to final position
+ * @return extracted block submatrix
+******************************************************************************/
+template <typename T>
+const NRSMat<T> NRSMat<T>::submatrix(const int from, const int to) const
+{
+#ifdef DEBUG
+        if(from<0 || from>=nn|| to<0 || to>=nn || from>to){
+                laerror("invalid submatrix specification");
+        }
+#endif
+
+NOT_GPU(*this);
+
+        const int n = to - from + 1;
+        NRSMat<T> r(n);
+	
+        for(int i=0; i<n; ++i)
+                {
+                int ii= i+from;
+                if(ii<0||ii>=nn) laerror("bad index in submatrix");
+                for(int j=0; j<=i; ++j)
+                        {
+                        int jj= j+from;
+                        r(i,j) = (*this)(ii,jj);
+                        }
+                }
+        return r;
+}
+
+
+template <typename T>
+const NRSMat<T> NRSMat<T>::submatrix(const NRVec<int> &selection) const
+{
+        NOT_GPU(*this);
+        const int n = selection.size();
+        NRSMat<T> r(n);
+
+        for(int i=0; i<n; ++i)
+                {
+                int ii=selection[i];
+                if(ii<0||ii>=nn) laerror("bad index in submatrix");
+		for(int j=0; j<=i; ++j)
+			{
+			int jj=selection[j];
+			r(i,j) = (*this)(ii,jj);
+			}
+                }
+return r;
+
+}
+
+/***************************************************************************//**
+ * places given matrix as submatrix at given position
+ * @param[in] fromrow row-coordinate of top left corner
+ * @param[in] fromcol col-coordinate of top left corner
+ * @param[in] rhs input matrix
+******************************************************************************/
+
+template <typename T>
+void NRSMat<T>::storesubmatrix(const int from, const NRSMat &rhs)
+{
+#ifdef DEBUG
+        if(from<0 || from>=nn){
+                laerror("invalid submatrix specification");
+        }
+#endif
+
+NOT_GPU(*this);
+
+        const int n = rhs.nrows();
+
+        for(int i=0; i<n; ++i)
+                {
+                int ii= i+from;
+                if(ii<0||ii>=nn) laerror("bad index in storesubmatrix");
+                for(int j=0; j<=i; ++j)
+                        {
+                        int jj= j+from;
+                        (*this)(ii,jj) = rhs(i,j);
+                        }
+                }
+}
+
+
+template <typename T>
+void NRSMat<T>::storesubmatrix(const NRVec<int> &selection, const NRSMat &rhs)
+{
+        NOT_GPU(*this);
+        const int n = selection.size();
+	if(rhs.size()!=n) laerror("size mismatch in storesubmatrix");
+
+        for(int i=0; i<n; ++i)
+                {
+                int ii=selection[i];
+                if(ii<0||ii>=nn) laerror("bad index in storesubmatrix");
+                for(int j=0; j<=i; ++j)
+                        {
+                        int jj=selection[j];
+                        (*this)(ii,jj) = rhs(i,j);
+                        }
+                }
+}
+
+
+
+
 /***************************************************************************//**
  * forced instantization in the corresponding object file
  ******************************************************************************/
diff --git a/smat.h b/smat.h
index 621e33a..f896fa2 100644
--- a/smat.h
+++ b/smat.h
@@ -167,6 +167,15 @@ public:
 	void resize(const int n);
 	void dealloc(void) {resize(0);}
 
+        //! extract specified submatrix
+        const NRSMat submatrix(const int from, const int to) const;
+        const NRSMat submatrix(const NRVec<int> &selection) const;
+
+        //! store given matrix at given position into the current matrix
+        void storesubmatrix(const int from, const NRSMat &rhs);
+        void storesubmatrix(const NRVec<int> &selection, const NRSMat &rhs);
+
+
 	inline operator T*();
 	inline operator const T*() const;
 	void fprintf(FILE *f, const char *format, const int modulo) const; 
diff --git a/vec.cc b/vec.cc
index 66064c4..08c5109 100644
--- a/vec.cc
+++ b/vec.cc
@@ -837,6 +837,24 @@ const NRVec<T> NRVec<T>::subvector(const int from, const int to) const
 	return r;
 }
 
+
+template <typename T>
+const NRVec<T> NRVec<T>::subvector(const NRVec<int> &selection) const
+{
+        NOT_GPU(*this);
+        const int n = selection.size();
+        NRVec<T> r(n);
+
+        for(int i=0; i<n; ++i)
+                {
+                int ii=selection[i];
+                if(ii<0||ii>=nn) laerror("bad row index in subvector");
+                r[i] = (*this)[ii];
+                }
+return r;
+
+}
+
 /***************************************************************************//**
  * places given vector as subvector at given position
  * @param[in] from coordinate 
@@ -866,6 +884,22 @@ void NRVec<T>::storesubvector(const int from, const NRVec &rhs)
 }
 
 
+template <typename T>
+void NRVec<T>::storesubvector(const NRVec<int> &selection, const NRVec &rhs)
+{
+        NOT_GPU(*this);
+        const int n = selection.size();
+	if(n!=rhs.size()) laerror("size mismatch in storesubvector");
+
+        for(int i=0; i<n; ++i)
+                {
+                int ii=selection[i];
+                if(ii<0||ii>=nn) laerror("bad  index in storesubvector");
+                (*this)[ii] = rhs[i];
+                }
+}
+
+
 /***************************************************************************//**
  * forced instantization in the corespoding object file
  ******************************************************************************/
diff --git a/vec.h b/vec.h
index b3070ed..ae3af15 100644
--- a/vec.h
+++ b/vec.h
@@ -192,9 +192,11 @@ public:
 
         //! extract specified subvector
         const NRVec subvector(const int from, const int to) const;
+	const NRVec subvector(const NRVec<int> &selection) const;
 
         //! store given vector at given position into the current vector
         void storesubvector(const int from, const NRVec &rhs);
+	void storesubvector(const NRVec<int> &selection, const NRVec &rhs);
 
 	//! relational operators 
 	const bool operator!=(const NRVec &rhs) const