From f7ef134fc54bc2e5b01407d8c4280e636a4ec5ec Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Mon, 26 Oct 2026 15:11:45 +0100
Subject: [PATCH] bugfixes in vec mat smat - mainly missing copyonwrite()s

---
 mat.cc  | 25 ++++++++++++++++++-------
 mat.h   |  1 +
 smat.cc |  9 +++++++--
 vec.cc  |  2 ++
 vec.h   |  4 +++-
 5 files changed, 31 insertions(+), 10 deletions(-)

diff --git a/mat.cc b/mat.cc
index 40c1762..ace1db3 100644
--- a/mat.cc
+++ b/mat.cc
@@ -127,6 +127,7 @@ void NRMat<T>::rowset(const NRVec<T> &r, const int i, int l, int offset) {
 	if(i < 0 || i >= nn) laerror("illegal index");
 	if(offset<0||l+offset>mm) laerror("illegal len/offset");
 #endif
+	copyonwrite();
 	LA_traits<T>::copy(
 #ifdef MATPTR
 	        v[i]+offset
@@ -893,6 +894,7 @@ void NRMat<T>::storesubmatrix(const int fromrow, const int fromcol, const NRMat
 #endif
 	SAME_LOC(*this, rhs);
 	if(!LA_traits<T>::is_plaindata()) laerror("only implemented for plain data");
+	copyonwrite();
 
 	const int m = tocol - fromcol + 1;
 	for(register int i = fromrow; i <= torow; ++i){
@@ -921,7 +923,8 @@ void NRMat<T>::storesubmatrix(const NRVec<int> &rows, const NRVec<int> &cols, co
         const int n = rows.size();
         const int m = cols.size();
 	if(rhs.nrows()!=n || rhs.ncols()!=m) laerror("incompatible dimensions in storesubmatrix");
-                
+	copyonwrite();
+
         for(int i=0; i<n; ++i)
                 {
                 int ii=rows[i];
@@ -942,7 +945,7 @@ void NRMat<T>::storesubmatrix(const NRVec<int> &rows, const NRVec<int> &cols, co
  ******************************************************************************/
 template <typename T>
 NRMat<T>& NRMat<T>::transposeme(const int _n) {
-	const int n = (n <= 0)?nn:_n;//!< transpose the entire matrix
+	const int n = (_n <= 0)?nn:_n;//!< transpose the entire matrix
 #ifdef DEBUG
 	if (n==nn && nn != mm || n>mm || n>nn ) laerror("NRMat<T>::transposeme() - invalid parameter n. Non-square matrix?");
 #endif
@@ -1184,23 +1187,27 @@ void NRMat<T>::fscanf(FILE *f, const char *format) {
 	int n(0), m(0);
 	if (::fscanf(f, "%d %d", &n, &m) != 2) laerror("cannot read matrix dimensions");
 
+	resize(n, m);
 #ifdef CUDALA
 	if(location == cpu){
-		p = *this;
+#endif
+		p = *this; //resize detached any shared data
+#ifdef CUDALA
 	}else{
-		tmp = new NRMat<T>(n, m, this->location);
+		tmp = new NRMat<T>(n, m, cpu); //read on cpu, then copied to the device
 		p = *tmp;
 	}
 #endif
-	resize(n, m);
 	for(int i=0; i<n; i++){
-		for(int j=0; j<n; j++){
+		for(int j=0; j<m; j++){
 			if(::fscanf(f,format, p++) != 1) laerror("cannot read matrix element");
 		}
 	}
 
 #ifdef CUDALA
 	if(location != cpu){
+		tmp->moveto(this->location);
+		*this |= *tmp;
 		delete tmp;
 	}
 #endif
@@ -1402,6 +1409,7 @@ void NRMat<double>::randomize(const double &x) {
 #ifdef CUDALA
 	if(location == cpu){
 #endif
+		copyonwrite();
 		for(register int i=0; i<nn; ++i){
 			for(register int j=0; j<mm; ++j){
 				(*this)(i,j) = x*RANDDOUBLESIGNED();
@@ -1431,6 +1439,7 @@ void NRMat<std::complex<double> >::randomize(const double &x) {
 #ifdef CUDALA
 	if(location == cpu){
 #endif
+		copyonwrite();
 		for(register int i=0; i<nn; ++i){
 			for(register int j=0; j<mm; ++j){
 				const double re = x*RANDDOUBLESIGNED();
@@ -2641,7 +2650,8 @@ int nnmin= nn<=mm?nn:mm;
 template<>
 void NRMat<double>::orthonormalize(const bool rowcol, const NRSMat<double> *metric) {
 
-	SAME_LOC(*this, *metric);
+	if(metric) SAME_LOC(*this, *metric);
+	copyonwrite();
 	if(metric){
 		if(rowcol){ //vectors are stored in rows
 			if((*metric).nrows() != mm) laerror("incompatible metric in NRMat<double>::orthonormalize(rowcol = true)");
@@ -3271,6 +3281,7 @@ template<typename T>
 void NRMat<T>::axpy(const T alpha, const NRPerm<int> &p, const bool direction)
 {
 int n=p.size();
+copyonwrite();
 for(int i=0; i<n; ++i)
 	{
 	if(direction) (*this)(i,p[i+1]-1) +=alpha;
diff --git a/mat.h b/mat.h
index d10e31d..d728bd5 100644
--- a/mat.h
+++ b/mat.h
@@ -291,6 +291,7 @@ public:
 	void columnset(const NRVec<T> &r, const int j, int l = -1) {
 		NOT_GPU(*this);
 		if(l < 0) l = nn;
+		copyonwrite();
 		for(register int i=0; i<l; ++i) (*this)(i,j) = r[i];
 	};
 
diff --git a/smat.cc b/smat.cc
index 0e51207..bbb3b9e 100644
--- a/smat.cc
+++ b/smat.cc
@@ -258,6 +258,7 @@ const T NRSMat<T>::trace() const {
 template<>
 void NRSMat<double>::randomize(const double &x) {
 	NOT_GPU(*this);
+	copyonwrite();
 
 	for(size_t i=0; i<NN2; ++i){
 		v[i] = x*RANDDOUBLESIGNED();
@@ -270,6 +271,7 @@ void NRSMat<double>::randomize(const double &x) {
  ******************************************************************************/
 template<>
 void NRSMat<std::complex<double> >::randomize(const double &x) {
+	copyonwrite();
 	for(register size_t i=0; i<NN2; ++i) v[i].real(x*RANDDOUBLESIGNED());
 	for(register size_t i=0; i<NN2; ++i) v[i].imag(x*RANDDOUBLESIGNED());
 	for(register int i=0; i<nn; ++i){
@@ -308,7 +310,8 @@ void NRSMat<T>::fscanf(FILE *f, const char *format) {
 		laerror("cannot read matrix dimensions in NRSMat<T>::fscanf(FILE *, const char *)");
 	if (n != m) laerror("different dimensions in NRSMat<T>::fscanf(FILE *, const char *)");
 	resize(n);
-	for (int i=0; i<n; i++) 
+	//the whole square matrix is read on purpose, the data are assumed to be symmetric (elements (i,j) and (j,i) share the storage)
+	for (int i=0; i<n; i++)
 		for (int j=0; j<n; j++)
 			if (::fscanf(f,format,&((*this)(i,j))) != 1)
 				laerror("NRSMat<T>::fscanf(FILE *, const char *) - unable to read matrix element");
@@ -873,6 +876,7 @@ void NRSMat<T>::storesubmatrix(const int from, const NRSMat &rhs)
 #endif
 
 NOT_GPU(*this);
+copyonwrite();
 
         const int n = rhs.nrows();
 
@@ -894,7 +898,8 @@ 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");
+	if(rhs.nrows()!=n) laerror("size mismatch in storesubmatrix");
+	copyonwrite();
 
         for(int i=0; i<n; ++i)
                 {
diff --git a/vec.cc b/vec.cc
index b861853..6b34ca8 100644
--- a/vec.cc
+++ b/vec.cc
@@ -192,6 +192,7 @@ const bool NRVec<T>::operator<(const NRVec &rhs) const {
 template<>
 void NRVec<double>::randomize(const double &x){
 	NOT_GPU(*this);
+	copyonwrite();
 
 	for(register int i=0; i<nn; ++i){
 		v[i] = x*RANDDOUBLESIGNED();
@@ -209,6 +210,7 @@ void NRVec<double>::randomize(const double &x){
 template<>
 void NRVec<std::complex<double> >::randomize(const double &x) {
 	NOT_GPU(*this);
+	copyonwrite();
 
 	for(register int i=0; i<nn; ++i){
 		v[i] = std::complex<double>(x*RANDDOUBLESIGNED(), x*RANDDOUBLESIGNED());
diff --git a/vec.h b/vec.h
index 245a783..3992a7f 100644
--- a/vec.h
+++ b/vec.h
@@ -202,7 +202,7 @@ public:
         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
+        //! store given vector at given position into the current vector (scatter)
         void storesubvector(const int from, const NRVec &rhs);
 	void storesubvector(const NRVec<int> &selection, const NRVec &rhs);
 
@@ -2294,6 +2294,7 @@ void NRVec<T>::storesubvector(const int from, const NRVec &rhs)
 #endif
 	SAME_LOC(*this, rhs);
 	if(!LA_traits<T>::is_plaindata()) laerror("only implemented for plain data");
+	copyonwrite();
 
 	#ifdef CUDALA
 		if(location == cpu){
@@ -2315,6 +2316,7 @@ 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");
+	copyonwrite();
 
         for(int i=0; i<n; ++i)
                 {