diff --git a/mat.cc b/mat.cc
index aa1e482..0c0fae2 100644
--- a/mat.cc
+++ b/mat.cc
@@ -755,7 +755,7 @@ void NRMat<double>::diagmultr(const NRVec<double> &rhs)
 	if (mm != rhs.size()) laerror("incompatible matrix dimension in diagmultr");
 #endif
 	copyonwrite();
-	for (int i=0; i<mm; i++) cblas_dscal(nn, rhs[i], (*this)[i], mm);
+	for (int i=0; i<mm; i++) cblas_dscal(nn, rhs[i], &(*this)(0,i), mm);
 }
 
 
@@ -767,7 +767,7 @@ void NRMat< complex<double> >::diagmultr(const NRVec< complex<double> > &rhs)
 	if (mm != rhs.size()) laerror("incompatible matrix dimension in diagmultl");
 #endif
 	copyonwrite();
-	for (int i=0; i<mm; i++) cblas_zscal(nn, &rhs[i], (*this)[i], mm);
+	for (int i=0; i<mm; i++) cblas_zscal(nn, &rhs[i], &(*this)(0,i), mm);
 }
 
 
@@ -1049,62 +1049,11 @@ return divide?NULL:&r[0];
 
 
 
-//////////////////////////////////////////////////////////////////////////////
-//// forced instantization in the corespoding object file
-#define INSTANTIZE(T) \
-template ostream & operator<<(ostream &s, const NRMat< T > &x); \
-template istream & operator>>(istream  &s, NRMat< T > &x); \
-
-INSTANTIZE(double)
-INSTANTIZE(complex<double>)
-INSTANTIZE(int)
-INSTANTIZE(short)
-INSTANTIZE(char)
-INSTANTIZE(unsigned char)
-INSTANTIZE(unsigned long)
-INSTANTIZE(unsigned int)
-
-
-export template <typename T>
-ostream& operator<<(ostream &s, const NRMat<T> &x)
-                {
-                int i,j,n,m;
-                n=x.nrows();
-                m=x.ncols();
-                s << n << ' ' << m << '\n';
-                for(i=0;i<n;i++)
-                        {
-                        for(j=0; j<m;j++) s << (typename LA_traits_io<T>::IOtype) x[i][j] << (j==m-1 ? '\n' : ' '); // endl cannot be used in the conditional expression, since it is an overloaded function
-                        }
-                return s;
-                }
-
-export template <typename T>
-istream& operator>>(istream  &s, NRMat<T> &x)
-                {
-                int i,j,n,m;
-                s >> n >> m;
-                x.resize(n,m);
-		typename LA_traits_io<T>::IOtype tmp;
-                for(i=0;i<n;i++) for(j=0; j<m;j++) { s>>tmp; x[i][j]=tmp;}
-                return s;
-                }
 
 //direct sum and product (oplus, otimes) to be done
 
 
 
-
-
-
-
-
-
-
-
-
-
-
 //////////////////////////////////////////////////////////////////////////////
 //// forced instantization in the corresponding object file
 template class NRMat<double>;
diff --git a/mat.h b/mat.h
index 93209c0..85e4b0c 100644
--- a/mat.h
+++ b/mat.h
@@ -548,8 +548,31 @@ return r;
 
 
 // I/O
-template <typename T> extern ostream& operator<<(ostream &s, const NRMat<T> &x);
-template <typename T> extern istream& operator>>(istream  &s, NRMat<T> &x);
+template <typename T>
+ostream& operator<<(ostream &s, const NRMat<T> &x)
+                {
+                int i,j,n,m;
+                n=x.nrows();
+                m=x.ncols();
+                s << n << ' ' << m << '\n';
+                for(i=0;i<n;i++)
+                        {
+                        for(j=0; j<m;j++) s << (typename LA_traits_io<T>::IOtype) x[i][j] << (j==m-1 ? '\n' : ' '); // endl cannot be used in the conditional expression, since it is an overloaded function
+                        }
+                return s;
+                }
+
+template <typename T>
+istream& operator>>(istream  &s, NRMat<T> &x)
+                {
+                int i,j,n,m;
+                s >> n >> m;
+                x.resize(n,m);
+		typename LA_traits_io<T>::IOtype tmp;
+                for(i=0;i<n;i++) for(j=0; j<m;j++) { s>>tmp; x[i][j]=tmp;}
+                return s;
+                }
+
 
 //optional indexing from 1
 //all possible constructors have to be given explicitly, other stuff is inherited
diff --git a/smat.cc b/smat.cc
index 1d27661..ab4c9dc 100644
--- a/smat.cc
+++ b/smat.cc
@@ -358,50 +358,6 @@ void NRSMat< complex<double> >::axpy(const complex<double> alpha,
 }
 
 
-export template <typename T>
-ostream& operator<<(ostream &s, const NRSMat<T> &x)
-                {
-                int i,j,n;
-                n=x.nrows();
-                s << n << ' ' << n << '\n';
-                for(i=0;i<n;i++)
-                        {
-                        for(j=0; j<n;j++) s << (typename LA_traits_io<T>::IOtype)x(i,j) << (j==n-1 ? '\n' : ' ');
-                        }
-                return s;
-                }
-
-
-export template <typename T>
-istream& operator>>(istream  &s, NRSMat<T> &x)
-                {
-                int i,j,n,m;
-                s >> n >> m;
-                if(n!=m) laerror("input symmetric matrix not square");
-                x.resize(n);
-		typename LA_traits_io<T>::IOtype tmp;
-                for(i=0;i<n;i++) for(j=0; j<m;j++) {s>>tmp; x(i,j)=tmp;}
-                return s;
-                }
-
-
-
-
-//////////////////////////////////////////////////////////////////////////////
-//// forced instantization in the corespoding object file
-#define INSTANTIZE(T) \
-template ostream & operator<<(ostream &s, const NRSMat< T > &x); \
-template istream & operator>>(istream  &s, NRSMat< T > &x); \
-
-INSTANTIZE(double)
-INSTANTIZE(complex<double>)
-INSTANTIZE(int)
-INSTANTIZE(short)
-INSTANTIZE(char)
-INSTANTIZE(unsigned int)
-INSTANTIZE(unsigned long)
-
-
 
 //////////////////////////////////////////////////////////////////////////////
 ////// forced instantization in the corresponding object file
diff --git a/smat.h b/smat.h
index f13cf5d..d642936 100644
--- a/smat.h
+++ b/smat.h
@@ -514,8 +514,32 @@ return r;
 }
 
 // I/O
-template <typename T> extern ostream& operator<<(ostream &s, const NRSMat<T> &x);
-template <typename T> extern istream& operator>>(istream  &s, NRSMat<T> &x);
+template <typename T>
+ostream& operator<<(ostream &s, const NRSMat<T> &x)
+                {
+                int i,j,n;
+                n=x.nrows();
+                s << n << ' ' << n << '\n';
+                for(i=0;i<n;i++)
+                        {
+                        for(j=0; j<n;j++) s << (typename LA_traits_io<T>::IOtype)x(i,j) << (j==n-1 ? '\n' : ' ');
+                        }
+                return s;
+                }
+
+
+template <typename T>
+istream& operator>>(istream  &s, NRSMat<T> &x)
+                {
+                int i,j,n,m;
+                s >> n >> m;
+                if(n!=m) laerror("input symmetric matrix not square");
+                x.resize(n);
+		typename LA_traits_io<T>::IOtype tmp;
+                for(i=0;i<n;i++) for(j=0; j<m;j++) {s>>tmp; x(i,j)=tmp;}
+                return s;
+                }
+
 
 // generate operators: SMat + a, a + SMat, SMat * a
 NRVECMAT_OPER(SMat,+)
diff --git a/sparsemat.cc b/sparsemat.cc
index e9ab6b6..dfba048 100644
--- a/sparsemat.cc
+++ b/sparsemat.cc
@@ -19,48 +19,6 @@ static inline void SWAP(T &a, T &b)
 
 
 
-export template <class T>
-ostream& operator<<(ostream &s, const SparseMat<T> &x)
-                {
-                SPMatindex n,m;
-                n=x.nrows();
-                m=x.ncols();
-                s << (int)n << ' ' << (int)m << '\n';
-		matel<T> *list=x.getlist();
-		while(list)
-                        {
-                        s << (int)list->row << ' ' << (int)list->col  << ' ' << (typename LA_traits_io<T>::IOtype)list->elem << '\n';
-			list=list->next;
-                        }
-		s << "-1 -1\n";
-                return s;
-                }
-
-export template <class T>
-istream& operator>>(istream  &s, SparseMat<T> &x)
-                {
-                int i,j;
-		int n,m;
-		matel<T> *l=NULL;
-                s >> n >> m;
-                x.resize(n,m);
-		s >> i >> j;
-		while(i>=0 && j>=0)
-			{
-			matel<T> *ll = l;
-			l= new matel<T>;
-			l->next= ll;
-			l->row=i;
-			l->col=j;
-			typename LA_traits_io<T>::IOtype tmp;
-			s >> tmp;
-			l->elem=tmp;	
-			s >> i >> j;
-			}
-		x.setlist(l);
-                return s;
-                }
-
 extern "C" {
 extern ssize_t read(int, void *, size_t);
 extern ssize_t write(int, const void *, size_t);
@@ -1309,8 +1267,6 @@ return *this;
 template SparseMat<T> & SparseMat<T>::oplusequal(const SparseMat<T> &rhs);\
 template SparseMat<T> & SparseMat<T>::oplusequal(const NRMat<T> &rhs);\
 template SparseMat<T> & SparseMat<T>::oplusequal(const NRSMat<T> &rhs);\
-template ostream& operator<<(ostream &s, const SparseMat<T> &x); \
-template istream& operator>>(istream &s, SparseMat<T> &x); \
 template void SparseMat<T>::get(int fd, bool dimen, bool transp); \
 template void SparseMat<T>::put(int fd, bool dimen, bool transp) const; \
 template void SparseMat<T>::copyonwrite(); \
diff --git a/sparsemat.h b/sparsemat.h
index 256c258..1cf9ad2 100644
--- a/sparsemat.h
+++ b/sparsemat.h
@@ -137,11 +137,48 @@ template <typename T>
 inline const NRMat<T> SparseMat<T>::operator*(const NRMat<T> &rhs) const
 {NRMat<T> result(nn,rhs.ncols()); result.gemm((T)0,*this,'n',rhs,'n',(T)1); return result;};
 
-template <typename T>
-	extern istream& operator>>(istream  &s, SparseMat<T> &x);
+template <class T>
+ostream& operator<<(ostream &s, const SparseMat<T> &x)
+                {
+                SPMatindex n,m;
+                n=x.nrows();
+                m=x.ncols();
+                s << (int)n << ' ' << (int)m << '\n';
+		matel<T> *list=x.getlist();
+		while(list)
+                        {
+                        s << (int)list->row << ' ' << (int)list->col  << ' ' << (typename LA_traits_io<T>::IOtype)list->elem << '\n';
+			list=list->next;
+                        }
+		s << "-1 -1\n";
+                return s;
+                }
+
+template <class T>
+istream& operator>>(istream  &s, SparseMat<T> &x)
+                {
+                int i,j;
+		int n,m;
+		matel<T> *l=NULL;
+                s >> n >> m;
+                x.resize(n,m);
+		s >> i >> j;
+		while(i>=0 && j>=0)
+			{
+			matel<T> *ll = l;
+			l= new matel<T>;
+			l->next= ll;
+			l->row=i;
+			l->col=j;
+			typename LA_traits_io<T>::IOtype tmp;
+			s >> tmp;
+			l->elem=tmp;	
+			s >> i >> j;
+			}
+		x.setlist(l);
+                return s;
+                }
 
-template <typename T>
-	extern ostream& operator<<(ostream &s, const SparseMat<T> &x);
 
 //destructor
 template <typename T>
diff --git a/vec.cc b/vec.cc
index 9891398..7fe10e9 100644
--- a/vec.cc
+++ b/vec.cc
@@ -15,8 +15,6 @@ extern ssize_t write(int, const void *, size_t);
 //////////////////////////////////////////////////////////////////////////////
 //// forced instantization in the corespoding object file
 #define INSTANTIZE(T) \
-template ostream & operator<<(ostream &s, const NRVec< T > &x); \
-template istream & operator>>(istream  &s, NRVec< T > &x); \
 template void NRVec<T>::put(int fd, bool dim, bool transp) const; \
 template void NRVec<T>::get(int fd, bool dim, bool transp); \
 
@@ -54,30 +52,6 @@ NRVec<T>::NRVec(const NRMat<T> &rhs)
 
 
 
-// formatted I/O
-template <typename T>
-ostream & operator<<(ostream &s, const NRVec<T> &x)
-{
-  int i, n;
-
-  n = x.size();
-  s << n << endl;
-  for(i=0; i<n; i++) s << (typename LA_traits_io<T>::IOtype)x[i] << (i == n-1 ? '\n' : ' ');
-  return s;
-}
-
-template <typename T>
-istream & operator>>(istream &s, NRVec<T> &x)
-{
-  int i,n;
-
-  s >> n;
-  x.resize(n);
-  typename LA_traits_io<T>::IOtype tmp;
-  for(i=0; i<n; i++) {s >> tmp; x[i]=tmp;}
-  return s;
-}
-
 
 //raw I/O
 template <typename T>
diff --git a/vec.h b/vec.h
index e01fb7b..a8c40c5 100644
--- a/vec.h
+++ b/vec.h
@@ -119,8 +119,30 @@ public:
 #include "smat.h"
 #include "sparsemat.h"
 
-template <typename T> ostream & operator<<(ostream &s, const NRVec<T> &x);
-template <typename T> istream & operator>>(istream  &s, NRVec<T> &x);
+// formatted I/O
+template <typename T>
+ostream & operator<<(ostream &s, const NRVec<T> &x)
+{
+  int i, n;
+
+  n = x.size();
+  s << n << endl;
+  for(i=0; i<n; i++) s << (typename LA_traits_io<T>::IOtype)x[i] << (i == n-1 ? '\n' : ' ');
+  return s;
+}
+
+template <typename T>
+istream & operator>>(istream &s, NRVec<T> &x)
+{
+  int i,n;
+
+  s >> n;
+  x.resize(n);
+  typename LA_traits_io<T>::IOtype tmp;
+  for(i=0; i<n; i++) {s >> tmp; x[i]=tmp;}
+  return s;
+}
+
 
 // INLINES