diff --git a/la_traits.h b/la_traits.h
index 6191724..0c221c5 100644
--- a/la_traits.h
+++ b/la_traits.h
@@ -45,6 +45,7 @@ class isscalar {
 
 //specializations
 #define SCALAR(X) \
+template<>\
 class isscalar<X> {typedef scalar_true scalar_type;};
 
 //declare what is scalar
diff --git a/mat.h b/mat.h
index eae57dc..548d053 100644
--- a/mat.h
+++ b/mat.h
@@ -24,8 +24,10 @@ public:
 	NRMat(const T *a, const int n, const int m);
 	inline NRMat(const NRMat &rhs);
 	explicit NRMat(const NRSMat<T> &rhs);
-#ifndef MATPTR
-	NRMat(const NRVec<T> &rhs, const int n, const int m);
+#ifdef MATPTR
+	explicit NRMat(const NRVec<T> &rhs, const int n, const int m) :NRMat(&rhs[0][0],n,m) {};
+#else
+	explicit NRMat(const NRVec<T> &rhs, const int n, const int m);
 #endif
 	~NRMat();
 #ifdef MATPTR
@@ -329,6 +331,7 @@ inline NRMat<T>::operator const T* () const
 }
 
 // max element of Mat
+template<>
 inline const double  NRMat<double>::amax() const
 {
 #ifdef MATPTR
@@ -337,6 +340,7 @@ inline const double  NRMat<double>::amax() const
 	return v[cblas_idamax(nn*mm, v, 1)];
 #endif
 }
+template<>
 inline const complex<double>  NRMat< complex<double> >::amax() const
 {
 #ifdef MATPTR
diff --git a/smat.h b/smat.h
index ffd0708..5c139cc 100644
--- a/smat.h
+++ b/smat.h
@@ -14,7 +14,7 @@ public:
 	friend class NRVec<T>;
 	friend class NRMat<T>;
 	
-	inline NRSMat<T>::NRSMat() : nn(0),v(0),count(0) {};
+	inline NRSMat() : nn(0),v(0),count(0) {};
 	inline explicit NRSMat(const int n);			// Zero-based array
 	inline NRSMat(const T &a, const int n);	//Initialize to constant
 	inline NRSMat(const T *a, const int n);	// Initialize to array
@@ -121,12 +121,15 @@ NRSMat<T>::NRSMat(const NRVec<T> &rhs, const int n) // type conversion
 }
 
 // S *= a
+template<>
 inline NRSMat<double> & NRSMat<double>::operator*=(const double & a)
 {
 	copyonwrite();
 	cblas_dscal(NN2, a, v, 1);
 	return *this;
 }
+
+template<>
 inline NRSMat< complex<double> > &
 NRSMat< complex<double> >::operator*=(const complex<double> & a)
 {
@@ -163,6 +166,7 @@ inline NRSMat<T> & NRSMat<T>::operator-=(const T &a)
 }
 
 // S += S
+template<>
 inline NRSMat<double> &
 NRSMat<double>::operator+=(const NRSMat<double> & rhs)
 {
@@ -173,6 +177,8 @@ NRSMat<double>::operator+=(const NRSMat<double> & rhs)
 	cblas_daxpy(NN2, 1.0, rhs.v, 1, v, 1);
 	return *this;
 }
+
+template<>
 NRSMat< complex<double> > &
 NRSMat< complex<double> >::operator+=(const NRSMat< complex<double> > & rhs)
 {
@@ -197,6 +203,7 @@ inline NRSMat<T> & NRSMat<T>::operator+=(const NRSMat<T> & rhs)
 
 
 // S -= S
+template<>
 inline NRSMat<double> &
 NRSMat<double>::operator-=(const NRSMat<double> & rhs)
 {
@@ -207,6 +214,8 @@ NRSMat<double>::operator-=(const NRSMat<double> & rhs)
 	cblas_daxpy(NN2, -1.0, rhs.v, 1, v, 1);
 	return *this;
 }
+
+template<>
 inline NRSMat< complex<double> > &
 NRSMat< complex<double> >::operator-=(const NRSMat< complex<double> > & rhs)
 {
@@ -306,10 +315,12 @@ inline int NRSMat<T>::size() const
 
 
 // max value
+template<>
 inline const double NRSMat<double>::amax() const
 {
 	return v[cblas_idamax(NN2, v, 1)];
 }
+template<>
 inline const complex<double> NRSMat< complex<double> >::amax() const
 {
 	return v[cblas_izamax(NN2, (void *)v, 1)];
diff --git a/vec.cc b/vec.cc
index 377938c..c115c04 100644
--- a/vec.cc
+++ b/vec.cc
@@ -4,6 +4,7 @@
 #include <sys/types.h>
 #include <sys/stat.h>
 #include <fcntl.h>
+#include <errno.h>
 extern "C" {
 extern ssize_t read(int, void *, size_t);
 extern ssize_t write(int, const void *, size_t);
@@ -171,6 +172,7 @@ return nn<rhs.nn;
 
 
 // axpy call for T = double (not strided)
+template<>
 void NRVec<double>::axpy(const double alpha, const NRVec<double> &x)
 {
 #ifdef DEBUG
@@ -181,6 +183,7 @@ void NRVec<double>::axpy(const double alpha, const NRVec<double> &x)
 }
 
 // axpy call for T = complex<double> (not strided)
+template<>
 void NRVec< complex<double> >::axpy(const complex<double> alpha, 
 			const NRVec< complex<double> > &x)
 {
@@ -192,6 +195,7 @@ void NRVec< complex<double> >::axpy(const complex<double> alpha,
 }
 
 // axpy call for T = double (strided)
+template<>
 void NRVec<double>::axpy(const double alpha, const double *x, const int stride)
 {
 	copyonwrite();
@@ -199,6 +203,7 @@ void NRVec<double>::axpy(const double alpha, const double *x, const int stride)
 }
 
 // axpy call for T = complex<double> (strided)
+template<>
 void NRVec< complex<double> >::axpy(const complex<double> alpha, 
 			const complex<double> *x, const int stride)
 {
@@ -207,6 +212,7 @@ void NRVec< complex<double> >::axpy(const complex<double> alpha,
 }
 
 // unary minus
+template<>
 const NRVec<double> NRVec<double>::operator-() const
 {
 	NRVec<double> result(*this);
@@ -214,6 +220,8 @@ const NRVec<double> NRVec<double>::operator-() const
 	cblas_dscal(nn, -1.0, result.v, 1);
 	return result;
 }
+
+template<>
 const NRVec< complex<double> > 
 NRVec< complex<double> >::operator-() const
 {
@@ -236,6 +244,7 @@ NRVec<T> & NRVec<T>::operator=(const T &a)
 }
 
 // Normalization of NRVec<double>
+template<>
 NRVec<double> & NRVec<double>::normalize()
 {
 	double tmp;
@@ -251,6 +260,7 @@ NRVec<double> & NRVec<double>::normalize()
 }
 
 // Normalization of NRVec< complex<double> >
+template<>
 NRVec< complex<double> > & NRVec< complex<double> >::normalize()
 {
 	complex<double> tmp;
@@ -265,9 +275,13 @@ NRVec< complex<double> > & NRVec< complex<double> >::normalize()
 }
 
 //stubs for linkage
+template<>
 NRVec<int> &  NRVec<int>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
+template<>
 NRVec<short> &  NRVec<short>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
+template<>
 NRVec<char> &  NRVec<char>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
+template<>
 void NRVec<int>::gemv(const int beta,
                 const NRSMat<int> &A, const char trans,
                 const int alpha, const NRVec &x)
@@ -275,6 +289,7 @@ void NRVec<int>::gemv(const int beta,
 laerror("not yet implemented");
 }
 
+template<>
 void NRVec<short>::gemv(const short beta,
                 const NRSMat<short> &A, const char trans,
                 const short alpha, const NRVec &x)
@@ -283,6 +298,7 @@ laerror("not yet implemented");
 }
 
 
+template<>
 void NRVec<char>::gemv(const char beta,
                 const NRSMat<char> &A, const char trans,
                 const char alpha, const NRVec &x)
@@ -290,6 +306,7 @@ void NRVec<char>::gemv(const char beta,
 laerror("not yet implemented");
 }
 
+template<>
 void NRVec<int>::gemv(const int beta,
                 const NRMat<int> &A, const char trans,
                 const int alpha, const NRVec &x)
@@ -297,6 +314,7 @@ void NRVec<int>::gemv(const int beta,
 laerror("not yet implemented");
 }
 
+template<>
 void NRVec<short>::gemv(const short beta,
                 const NRMat<short> &A, const char trans,
                 const short alpha, const NRVec &x)
@@ -305,6 +323,7 @@ laerror("not yet implemented");
 }
 
 
+template<>
 void NRVec<char>::gemv(const char beta,
                 const NRMat<char> &A, const char trans,
                 const char alpha, const NRVec &x)
@@ -312,6 +331,7 @@ void NRVec<char>::gemv(const char beta,
 laerror("not yet implemented");
 }
 
+template<>
 void NRVec<int>::gemv(const int beta,
                 const SparseMat<int> &A, const char trans,
                 const int alpha, const NRVec &x)
@@ -319,6 +339,7 @@ void NRVec<int>::gemv(const int beta,
 laerror("not yet implemented");
 }
 
+template<>
 void NRVec<short>::gemv(const short beta,
                 const SparseMat<short> &A, const char trans,
                 const short alpha, const NRVec &x)
@@ -327,6 +348,7 @@ laerror("not yet implemented");
 }
 
 
+template<>
 void NRVec<char>::gemv(const char beta,
                 const SparseMat<char> &A, const char trans,
                 const char alpha, const NRVec &x)
@@ -338,6 +360,7 @@ laerror("not yet implemented");
 
 
 // gemv calls 
+template<>
 void NRVec<double>::gemv(const double beta, const NRMat<double> &A, 
 		const char trans, const double alpha, const NRVec &x)
 {
@@ -349,6 +372,7 @@ void NRVec<double>::gemv(const double beta, const NRMat<double> &A,
 			A.nrows(), A.ncols(), alpha, A, A.ncols(), x.v, 1, beta, v, 1);
 }
 
+template<>
 void NRVec< complex<double> >::gemv(const complex<double> beta,
 		const NRMat< complex<double> > &A, const char trans, 
 		const complex<double> alpha, const NRVec &x)
@@ -363,6 +387,7 @@ void NRVec< complex<double> >::gemv(const complex<double> beta,
 }
 
 
+template<>
 void NRVec<double>::gemv(const double beta, const NRSMat<double> &A,
                 const char trans, const double alpha, const NRVec &x)
 {
@@ -374,6 +399,7 @@ void NRVec<double>::gemv(const double beta, const NRSMat<double> &A,
 }
 
 
+template<>
 void NRVec< complex<double> >::gemv(const complex<double> beta,
                 const NRSMat< complex<double> > &A, const char trans,
                 const complex<double> alpha, const NRVec &x)
@@ -391,12 +417,14 @@ void NRVec< complex<double> >::gemv(const complex<double> beta,
 
 
 // Direc product Mat = Vec | Vec
+template<>
 const NRMat<double> NRVec<double>::operator|(const NRVec<double> &b) const
 {
 	NRMat<double> result(0.,nn,b.nn);
 	cblas_dger(CblasRowMajor, nn, b.nn, 1., v, 1, b.v, 1, result, b.nn);
 	return result;
 }
+template<>
 const NRMat< complex<double> > 
 NRVec< complex<double> >::operator|(const NRVec< complex<double> > &b) const
 {
diff --git a/vec.h b/vec.h
index ef306aa..8e59b97 100644
--- a/vec.h
+++ b/vec.h
@@ -45,7 +45,9 @@ public:
 	inline NRVec(const T *a, const int n);
 	inline NRVec(const NRVec &rhs);
 	inline explicit NRVec(const NRSMat<T> & S);
-#ifndef MATPTR
+#ifdef MATPTR
+	explicit NRVec(const NRMat<T> &rhs) : NRVec(&rhs[0][0],rhs.nrows()*rhs.ncols()) {};
+#else
 	explicit NRVec(const NRMat<T> &rhs);
 #endif
 	NRVec & operator=(const NRVec &rhs);
@@ -152,6 +154,7 @@ inline NRVec<T>::NRVec(const NRSMat<T> &rhs)
 }
 
 // x += a
+template<>
 inline NRVec<double> & NRVec<double>::operator+=(const double &a)
 {
 	copyonwrite();
@@ -159,6 +162,7 @@ inline NRVec<double> & NRVec<double>::operator+=(const double &a)
 	return *this;
 }
 
+template<>
 inline NRVec< complex<double> > &
 NRVec< complex<double> >::operator+=(const complex<double> &a)
 {
@@ -179,12 +183,15 @@ inline NRVec<T> & NRVec<T>::operator+=(const T &a)
 
 
 // x -= a
+template<>
 inline NRVec<double> & NRVec<double>::operator-=(const double &a)
 {
 	copyonwrite();
 	cblas_daxpy(nn, 1.0, &a, 0, v, 1);
 	return *this;
 }
+
+template<>
 inline NRVec< complex<double> > &
 NRVec< complex<double> >::operator-=(const complex<double> &a)
 {
@@ -205,6 +212,7 @@ inline NRVec<T> & NRVec<T>::operator-=(const T &a)
 
 
 // x += x
+template<>
 inline NRVec<double> & NRVec<double>::operator+=(const NRVec<double> &rhs)
 {
 #ifdef DEBUG
@@ -214,6 +222,8 @@ inline NRVec<double> & NRVec<double>::operator+=(const NRVec<double> &rhs)
 	cblas_daxpy(nn, 1.0, rhs.v, 1, v, 1);
 	return *this;
 }
+
+template<>
 inline NRVec< complex<double> > &
 NRVec< complex<double> >::operator+=(const NRVec< complex<double> > &rhs)
 {
@@ -240,6 +250,7 @@ inline NRVec<T> & NRVec<T>::operator+=(const NRVec<T> &rhs)
 
 
 // x -= x
+template<>
 inline NRVec<double> & NRVec<double>::operator-=(const NRVec<double> &rhs)
 {
 #ifdef DEBUG
@@ -249,6 +260,8 @@ inline NRVec<double> & NRVec<double>::operator-=(const NRVec<double> &rhs)
 	cblas_daxpy(nn, -1.0, rhs.v, 1, v, 1);
 	return *this;
 }
+
+template<>
 inline NRVec< complex<double> > &
 NRVec< complex<double> >::operator-=(const NRVec< complex<double> > &rhs)
 {
@@ -275,12 +288,15 @@ inline NRVec<T> & NRVec<T>::operator-=(const NRVec<T> &rhs)
 
 
 // x *= a
+template<>
 inline NRVec<double> & NRVec<double>::operator*=(const double &a)
 {
 	copyonwrite();
 	cblas_dscal(nn, a, v, 1);
 	return *this;
 }
+
+template<>
 inline NRVec< complex<double> > &
 NRVec< complex<double> >::operator*=(const complex<double> &a)
 {
@@ -301,6 +317,7 @@ inline NRVec<T> & NRVec<T>::operator*=(const T &a)
 
 
 // scalar product x.y
+template<>
 inline const double NRVec<double>::operator*(const NRVec<double> &rhs) const
 {
 #ifdef DEBUG
@@ -310,6 +327,7 @@ inline const double NRVec<double>::operator*(const NRVec<double> &rhs) const
 }
 
 
+template<>
 inline const complex<double>
 NRVec< complex<double> >::operator*(const NRVec< complex<double> > &rhs) const
 {
@@ -335,10 +353,12 @@ inline const T NRVec<T>::operator*(const NRVec<T> &rhs) const
 
 
 // Sum of elements
+template<>
 inline const double NRVec<double>::sum() const
 {
 	return cblas_dasum(nn, v, 1);
 }
+template<>
 inline const complex<double>
 NRVec< complex<double> >::sum() const
 {
@@ -348,10 +368,12 @@ NRVec< complex<double> >::sum() const
 }
 
 // Dot product: x * y
+template<>
 inline const double NRVec<double>::dot(const double *y, const int stride) const
 {
 	return cblas_ddot(nn, y, stride, v, 1);
 }
+template<>
 inline const complex<double>
 NRVec< complex<double> >::dot(const complex<double> *y, const int stride) const
 {
@@ -407,20 +429,24 @@ inline NRVec<T>::operator const T*() const
 }
 
 // return norm of the Vec
+template<>
 inline const double  NRVec<double>::norm() const
 {
 	return cblas_dnrm2(nn, v, 1);
 }
+template<>
 inline const double NRVec< complex<double> >::norm() const
 {
 	return cblas_dznrm2(nn, (void *)v, 1);
 }
 
 // Max element of the array
+template<>
 inline const double  NRVec<double>::amax() const
 {
 	return v[cblas_idamax(nn, v, 1)];
 }
+template<>
 inline const complex<double> NRVec< complex<double> >::amax() const
 {
 	return v[cblas_izamax(nn, (void *)v, 1)];