diff --git a/mat.cc b/mat.cc
index 0c0fae2..45febb9 100644
--- a/mat.cc
+++ b/mat.cc
@@ -243,6 +243,24 @@ for(int i=fromrow; i<=torow; ++i)
 return r;
 }
 
+template <typename T>
+void NRMat<T>::storesubmatrix(const int fromrow, const int fromcol, const NRMat &rhs)
+{
+int tocol=fromcol+rhs.ncols()-1;
+int torow=fromrow+rhs.nrows()-1;
+#ifdef DEBUG
+if(fromrow <0 ||fromrow >=nn||torow >=nn ||fromcol<0||fromcol>=mm||tocol>=mm) laerror("bad indices in storesubmatrix");
+#endif
+int n=torow-fromrow+1;
+int m=tocol-fromcol+1;
+for(int i=fromrow; i<=torow; ++i)
+#ifdef MATPTR
+        memcpy(v[i]+fromcol,rhs.v[i-fromrow],m*sizeof(T));
+#else
+        memcpy(v+i*mm+fromcol,rhs.v+(i-fromrow)*m,m*sizeof(T));
+#endif
+}
+
 
 // transpose Mat
 template <typename T>
diff --git a/mat.h b/mat.h
index 85e4b0c..97b132d 100644
--- a/mat.h
+++ b/mat.h
@@ -92,6 +92,7 @@ public:
 	const NRMat transpose(bool conj=false) const;
 	const NRMat conjugate() const;
         const NRMat submatrix(const int fromrow, const int torow, const int fromcol, const int tocol) const; //there is also independent less efficient routine for generally indexed submatrix
+	void storesubmatrix(const int fromrow, const int fromcol, const NRMat &rhs); //overwrite a block with external matrix
 	void gemm(const T &beta, const NRMat &a, const char transa, const NRMat &b,
 			const char transb, const T &alpha);//this = alpha*op( A )*op( B ) + beta*this
 	void fprintf(FILE *f, const char *format, const int modulo) const;