diff --git a/fourindex.h b/fourindex.h
index 1fc58f8..2a1f65c 100644
--- a/fourindex.h
+++ b/fourindex.h
@@ -103,15 +103,15 @@ public:
 		piterator(matel4<I,T> *pp): symmetry(nosymmetry),p(pp),permindex(0){};
                 ~piterator() {};
 		piterator(const fourindex &x): symmetry(x.symmetry),p(x.list),permindex(0) {setup();};
-		piterator& operator++() {if(++permindex==fourindex_permnumbers[symmetry]) {permindex=0; p=p->next;} setup(); return *this;}
+		piterator& operator++() {if(++permindex>=fourindex_permnumbers[symmetry]) {permindex=0; p=p->next;} setup(); return *this;}
                 const matel4<I,T> & operator*() const {return my;}
                 const matel4<I,T> * operator->() const {return &my;}
 		piterator operator++(int) {laerror("postincrement not possible on permute-iterator");}
-                bool operator==(const piterator &rhs) const {return p==rhs.p && permindex==rhs.permindex && symmetry==rhs.symmetry;}
-                bool operator!=(const piterator &rhs) const {return p!=rhs.p || permindex!=rhs.permindex || symmetry!=rhs.symmetry;}
+                bool operator==(const piterator &rhs) const {return p==rhs.p && (!p || permindex==rhs.permindex);}
+                bool operator!=(const piterator &rhs) const {return p!=rhs.p || p && rhs.p && permindex!=rhs.permindex;}
         };
-        piterator pbegin() const {return *this;}
-        piterator pend() const {return NULL;}
+        piterator pbegin() const {return piterator(*this);}
+        piterator pend() const {return piterator(NULL);}
 
 	//constructors etc.
 	inline fourindex() :nn(0),count(NULL),list(NULL) {};
diff --git a/la_traits.h b/la_traits.h
index 0c221c5..77bb29b 100644
--- a/la_traits.h
+++ b/la_traits.h
@@ -39,14 +39,12 @@ typedef class scalar_true {};
 
 //default is non-scalar
 template<typename C>
-class isscalar {
-	typedef scalar_false scalar_type;
-	};
+class isscalar { public: typedef scalar_false scalar_type;};
 
 //specializations
 #define SCALAR(X) \
 template<>\
-class isscalar<X> {typedef scalar_true scalar_type;};
+class isscalar<X> {public: typedef scalar_true scalar_type;};
 
 //declare what is scalar
 SCALAR(char)
diff --git a/mat.cc b/mat.cc
index fa6d3e8..e6469a4 100644
--- a/mat.cc
+++ b/mat.cc
@@ -3,6 +3,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);
@@ -12,11 +13,11 @@ extern ssize_t write(int, const void *, size_t);
 
 //////////////////////////////////////////////////////////////////////////////
 //// forced instantization in the corresponding object file
-template NRMat<double>;
-template NRMat< complex<double> >;
-template NRMat<int>;
-template NRMat<short>;
-template NRMat<char>;
+template class NRMat<double>;
+template class NRMat< complex<double> >;
+template class NRMat<int>;
+template class NRMat<short>;
+template class NRMat<char>;
 
 
 /*
@@ -288,12 +289,16 @@ void NRMat<T>::fscanf(FILE *f, const char *format)
  */
 
 // Mat *= a
+template<>
 NRMat<double> & NRMat<double>::operator*=(const double &a)
 {
 	copyonwrite();
 	cblas_dscal(nn*mm, a, *this, 1);
 	return *this;
 }
+
+
+template<>
 NRMat< complex<double> > &
 NRMat< complex<double> >::operator*=(const complex<double> &a)
 {
@@ -302,6 +307,7 @@ NRMat< complex<double> >::operator*=(const complex<double> &a)
 	return *this;
 }
 
+
 //and for general type
 template <typename T>
 NRMat<T> & NRMat<T>::operator*=(const T &a)
@@ -318,6 +324,7 @@ NRMat<T> & NRMat<T>::operator*=(const T &a)
 
 
 // Mat += Mat
+template<>
 NRMat<double> & NRMat<double>::operator+=(const NRMat<double>  &rhs)
 {
 #ifdef DEBUG
@@ -328,6 +335,9 @@ NRMat<double> & NRMat<double>::operator+=(const NRMat<double>  &rhs)
 	cblas_daxpy(nn*mm, 1.0, rhs, 1, *this, 1);
 	return *this;
 }
+
+
+template<>
 NRMat< complex<double> > &
 NRMat< complex<double> >::operator+=(const NRMat< complex<double> >  &rhs)
 {
@@ -340,6 +350,8 @@ NRMat< complex<double> >::operator+=(const NRMat< complex<double> >  &rhs)
 	 return *this;
 }
 
+
+
 //and for general type
 template <typename T>
 NRMat<T> & NRMat<T>::operator+=(const NRMat<T>  &rhs)
@@ -359,6 +371,7 @@ NRMat<T> & NRMat<T>::operator+=(const NRMat<T>  &rhs)
 
 
 // Mat -= Mat
+template<>
 NRMat<double> & NRMat<double>::operator-=(const NRMat<double>  &rhs)
 {
 #ifdef DEBUG
@@ -369,6 +382,10 @@ NRMat<double> & NRMat<double>::operator-=(const NRMat<double>  &rhs)
 	cblas_daxpy(nn*mm, -1.0, rhs, 1, *this, 1);
 	return *this;
 }
+
+
+
+template<>
 NRMat< complex<double> > &
 NRMat< complex<double> >::operator-=(const NRMat< complex<double> >  &rhs)
 {
@@ -381,6 +398,8 @@ NRMat< complex<double> >::operator-=(const NRMat< complex<double> >  &rhs)
 	 return *this;
 }
 
+
+
 //and for general type
 template <typename T>
 NRMat<T> & NRMat<T>::operator-=(const NRMat<T>  &rhs)
@@ -400,6 +419,7 @@ NRMat<T> & NRMat<T>::operator-=(const NRMat<T>  &rhs)
 
 
 // Mat += SMat
+template<>
 NRMat<double> & NRMat<double>::operator+=(const NRSMat<double> &rhs)
 {
 #ifdef DEBUG
@@ -418,6 +438,10 @@ NRMat<double> & NRMat<double>::operator+=(const NRSMat<double> &rhs)
 	}
 	return *this;
 }
+
+
+
+template<>
 NRMat< complex<double> > & 
 NRMat< complex<double> >::operator+=(const NRSMat< complex<double> > &rhs)
 {
@@ -438,6 +462,9 @@ NRMat< complex<double> >::operator+=(const NRSMat< complex<double> > &rhs)
 	return *this;
 }
 
+
+
+
 //and for general type
 template <typename T>
 NRMat<T> & NRMat<T>::operator+=(const NRSMat<T> &rhs)
@@ -461,6 +488,7 @@ NRMat<T> & NRMat<T>::operator+=(const NRSMat<T> &rhs)
 
 
 // Mat -= SMat
+template<>
 NRMat<double> & NRMat<double>::operator-=(const NRSMat<double> &rhs)
 {
 #ifdef DEBUG
@@ -479,6 +507,10 @@ NRMat<double> & NRMat<double>::operator-=(const NRSMat<double> &rhs)
 	}
 	return *this;
 }
+
+
+
+template<>
 NRMat< complex<double> > & 
 NRMat< complex<double> >::operator-=(const NRSMat< complex<double> > &rhs)
 {
@@ -521,7 +553,11 @@ NRMat<T> & NRMat<T>::operator-=(const NRSMat<T> &rhs)
         return *this;
 }
 
+
+
+
 // Mat.Mat - scalar product
+template<>
 const double NRMat<double>::dot(const NRMat<double> &rhs) const
 {
 #ifdef DEBUG
@@ -529,6 +565,10 @@ const double NRMat<double>::dot(const NRMat<double> &rhs) const
 #endif
 	return cblas_ddot(nn*mm, (*this)[0], 1, rhs[0], 1);
 }
+
+
+
+template<>
 const complex<double>
 NRMat< complex<double> >::dot(const NRMat< complex<double> > &rhs) const
 {
@@ -542,6 +582,7 @@ NRMat< complex<double> >::dot(const NRMat< complex<double> > &rhs) const
 }
 
 // Mat * Mat
+template<>
 const NRMat<double> NRMat<double>::operator*(const NRMat<double> &rhs) const
 {
 #ifdef DEBUG
@@ -552,6 +593,10 @@ const NRMat<double> NRMat<double>::operator*(const NRMat<double> &rhs) const
 			*this, mm, rhs, rhs.mm, 0.0, result, rhs.mm);
 	return result;
 }
+
+
+
+template<>
 const NRMat< complex<double> > 
 NRMat< complex<double> >::operator*(const NRMat< complex<double> > &rhs) const
 {
@@ -565,7 +610,9 @@ NRMat< complex<double> >::operator*(const NRMat< complex<double> > &rhs) const
 	return result;
 }
 
+
 // Multiply by diagonal from L
+template<>
 void NRMat<double>::diagmultl(const NRVec<double> &rhs)
 {
 #ifdef DEBUG
@@ -574,6 +621,10 @@ void NRMat<double>::diagmultl(const NRVec<double> &rhs)
 	copyonwrite();
 	for(int i=0; i<nn; i++) cblas_dscal(mm, rhs[i], (*this)[i], 1);
 }
+
+
+
+template<>
 void NRMat< complex<double> >::diagmultl(const NRVec< complex<double> > &rhs)
 {
 #ifdef DEBUG
@@ -583,7 +634,11 @@ void NRMat< complex<double> >::diagmultl(const NRVec< complex<double> > &rhs)
 	for (int i=0; i<nn; i++) cblas_zscal(mm, &rhs[i], (*this)[i], 1);
 }
 
+
+
+
 // Multiply by diagonal from R
+template<>
 void NRMat<double>::diagmultr(const NRVec<double> &rhs)
 {
 #ifdef DEBUG
@@ -592,6 +647,10 @@ void NRMat<double>::diagmultr(const NRVec<double> &rhs)
 	copyonwrite();
 	for (int i=0; i<mm; i++) cblas_dscal(nn, rhs[i], (*this)[i], mm);
 }
+
+
+
+template<>
 void NRMat< complex<double> >::diagmultr(const NRVec< complex<double> > &rhs)
 {
 #ifdef DEBUG
@@ -601,7 +660,11 @@ void NRMat< complex<double> >::diagmultr(const NRVec< complex<double> > &rhs)
 	for (int i=0; i<mm; i++) cblas_zscal(nn, &rhs[i], (*this)[i], mm);
 }
 
+
+
+
 // Mat * Smat, decomposed to nn x Vec * Smat
+template<>
 const NRMat<double> 
 NRMat<double>::operator*(const NRSMat<double> &rhs) const
 {
@@ -614,6 +677,10 @@ NRMat<double>::operator*(const NRSMat<double> &rhs) const
 				(*this)[i], 1, 0.0, result[i], 1);
 	return result;
 }
+
+
+
+template<>
 const NRMat< complex<double> >
 NRMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
 {
@@ -629,6 +696,7 @@ NRMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
 
 
 // sum of rows
+template<>
 const NRVec<double> NRMat<double>::rsum() const
 {
 	NRVec<double> result(mm);
@@ -637,6 +705,7 @@ const NRVec<double> NRMat<double>::rsum() const
 }
 
 // sum of columns
+template<>
 const NRVec<double> NRMat<double>::csum() const
 {
 	NRVec<double> result(nn);
@@ -645,8 +714,10 @@ const NRVec<double> NRMat<double>::csum() const
 }
 
 // complex conjugate of Mat
+template<>
 NRMat<double> &NRMat<double>::conjugateme() {return *this;}
 
+template<>
 NRMat< complex<double> > & NRMat< complex<double> >::conjugateme()
 {
 	copyonwrite();
@@ -655,12 +726,14 @@ NRMat< complex<double> > & NRMat< complex<double> >::conjugateme()
 }
 
 // transpose and optionally conjugate
+template<>
 const NRMat<double> NRMat<double>::transpose(bool conj) const
 {
 	NRMat<double> result(mm,nn);
 	for(int i=0; i<nn; i++) cblas_dcopy(mm, (*this)[i], 1, result[0]+i, nn);
 	return result;
 }
+template<>
 const NRMat< complex<double> >
 NRMat< complex<double> >::transpose(bool conj) const
 {
@@ -672,6 +745,7 @@ NRMat< complex<double> >::transpose(bool conj) const
 }
 
 // gemm : this = alpha*op( A )*op( B ) + beta*this
+template<>
 void NRMat<double>::gemm(const double &beta, const NRMat<double> &a,
 		const char transa, const NRMat<double> &b, const char transb, 
 		const double &alpha)
@@ -691,6 +765,10 @@ void NRMat<double>::gemm(const double &beta, const NRMat<double> &a,
 			(transb=='n' ? CblasNoTrans : CblasTrans), nn, mm, k, alpha, a,
 			a.mm, b , b.mm, beta, *this , mm);
 }
+
+
+
+template<>
 void NRMat< complex<double> >::gemm(const complex<double> & beta,
 		const NRMat< complex<double> > & a, const char transa, 
 		const NRMat< complex<double> > & b, const char transb, 
@@ -714,6 +792,7 @@ void NRMat< complex<double> >::gemm(const complex<double> & beta,
 }
 
 // norm of Mat
+template<>
 const double  NRMat<double>::norm(const double scalar) const
 {
 	if (!scalar) return cblas_dnrm2(nn*mm, (*this)[0], 1);
@@ -731,6 +810,10 @@ const double  NRMat<double>::norm(const double scalar) const
 		}
 	return sqrt(sum);
 }
+
+
+
+template<>
 const double  NRMat< complex<double> >::norm(const complex<double> scalar) const
 {
 	if (scalar == CZERO) return cblas_dznrm2(nn*mm, (*this)[0], 1);
@@ -749,7 +832,11 @@ const double  NRMat< complex<double> >::norm(const complex<double> scalar) const
 	return sqrt(sum);
 }
 
+
+
+
 // axpy: this = a * Mat
+template<>
 void NRMat<double>::axpy(const double alpha, const NRMat<double> &mat)
 {
 #ifdef DEBUG
@@ -758,6 +845,10 @@ void NRMat<double>::axpy(const double alpha, const NRMat<double> &mat)
 	copyonwrite();
 	cblas_daxpy(nn*mm, alpha, mat, 1, *this, 1);
 }
+
+
+
+template<>
 void NRMat< complex<double> >::axpy(const complex<double> alpha, 
 		const NRMat< complex<double> > & mat)
 {
@@ -768,14 +859,24 @@ void NRMat< complex<double> >::axpy(const complex<double> alpha,
 	cblas_zaxpy(nn*mm, (void *)&alpha, mat, 1, (void *)(*this)[0], 1);
 }
 
+
+
+
 // trace of Mat
+template<>
 const double NRMat<double>::trace() const
 {
 #ifdef DEBUG
 	if (nn != mm) laerror("no-square matrix in Mat::trace()");
 #endif
 	return cblas_dasum(nn, (*this)[0], nn+1);
+
+
 }
+
+
+
+template<>
 const complex<double> NRMat< complex<double> >::trace() const
 {
 #ifdef DEBUG
@@ -795,6 +896,7 @@ const complex<double> NRMat< complex<double> >::trace() const
 
 //get diagonal; for compatibility with large matrices do not return newly created object
 //for non-square get diagonal of A^TA, will be used as preconditioner 
+template<>
 void NRMat<double>::diagonalof(NRVec<double> &r, const bool divide) const
 {
 #ifdef DEBUG
diff --git a/smat.cc b/smat.cc
index fe13588..ecfc7f7 100644
--- a/smat.cc
+++ b/smat.cc
@@ -3,6 +3,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);
@@ -13,11 +14,11 @@ extern ssize_t write(int, const void *, size_t);
 
 //////////////////////////////////////////////////////////////////////////////
 ////// forced instantization in the corresponding object file
-template NRSMat<double>;
-template NRSMat< complex<double> >;
-template NRSMat<int>;
-template NRSMat<short>;
-template NRSMat<char>;
+template class NRSMat<double>;
+template class NRSMat< complex<double> >;
+template class NRSMat<int>;
+template class NRSMat<short>;
+template class NRSMat<char>;
 
 
 
@@ -146,6 +147,7 @@ void NRSMat<T>::fscanf(FILE *f, const char *format)
  */
 
 // SMat * Mat
+template<>
 const NRMat<double> NRSMat<double>::operator*(const NRMat<double> &rhs) const
 {
 #ifdef DEBUG
@@ -157,6 +159,9 @@ const NRMat<double> NRSMat<double>::operator*(const NRMat<double> &rhs) const
 				0.0, result[0]+k, rhs.ncols());
 	return result;
 }
+
+
+template<>
 const NRMat< complex<double> >
 NRSMat< complex<double> >::operator*(const NRMat< complex<double> > &rhs) const
 {
@@ -170,7 +175,9 @@ NRSMat< complex<double> >::operator*(const NRMat< complex<double> > &rhs) const
 	return result;
 }
 
+
 // SMat * SMat
+template<>
 const NRMat<double> NRSMat<double>::operator*(const NRSMat<double> &rhs) const
 {
 #ifdef DEBUG
@@ -218,6 +225,10 @@ const NRMat<double> NRSMat<double>::operator*(const NRSMat<double> &rhs) const
 
 	return result;
 }
+
+
+
+template<>
 const NRMat< complex<double> > 
 NRSMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
 {
@@ -230,7 +241,12 @@ NRSMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
 	return result;
 //	laerror("complex SMat*Smat not implemented");
 }
+
+
+
+
 // S dot S
+template<>
 const double NRSMat<double>::dot(const NRSMat<double> &rhs) const
 {
 #ifdef DEBUG
@@ -238,6 +254,10 @@ const double NRSMat<double>::dot(const NRSMat<double> &rhs) const
 #endif
 	return cblas_ddot(NN2, v, 1, rhs.v, 1);
 }
+
+
+
+template<>
 const complex<double> 
 NRSMat< complex<double> >::dot(const NRSMat< complex<double> > &rhs) const
 {
@@ -251,6 +271,7 @@ NRSMat< complex<double> >::dot(const NRSMat< complex<double> > &rhs) const
 
 
 // norm of the matrix
+template<>
 const double  NRSMat<double>::norm(const double scalar) const
 {
 	if (!scalar) return cblas_dnrm2(NN2, v, 1);
@@ -265,8 +286,11 @@ const double  NRSMat<double>::norm(const double scalar) const
 		}
 	return sqrt(sum);
 }
-const double
-NRSMat< complex<double> >::norm(const complex<double> scalar) const
+
+
+
+template<>
+const double NRSMat< complex<double> >::norm(const complex<double> scalar) const
 {
 	if (!(scalar.real()) && !(scalar.imag()))
 		return cblas_dznrm2(NN2, (void *)v, 1);
@@ -282,7 +306,12 @@ NRSMat< complex<double> >::norm(const complex<double> scalar) const
 	return sqrt(sum);
 }
 
+
+
+
+
 // axpy: S = S * a
+template<>
 void NRSMat<double>::axpy(const double alpha, const NRSMat<double> & x)
 {
 #ifdef DEBUG
@@ -291,6 +320,10 @@ void NRSMat<double>::axpy(const double alpha, const NRSMat<double> & x)
 	copyonwrite();
 	cblas_daxpy(NN2, alpha, x.v, 1, v, 1);
 }
+
+
+
+template<>
 void NRSMat< complex<double> >::axpy(const complex<double> alpha,
 			const NRSMat< complex<double> > & x)
 {
diff --git a/sparsemat.cc b/sparsemat.cc
index 6ed5d6f..8f7c18f 100644
--- a/sparsemat.cc
+++ b/sparsemat.cc
@@ -4,12 +4,13 @@
 #include <sys/types.h>
 #include <sys/stat.h>
 #include <fcntl.h>
+#include <errno.h>
 #include "sparsemat.h"
 
 //////////////////////////////////////////////////////////////////////////////
 //// forced instantization in the corresponding object file
-template SparseMat<double>;
-template SparseMat<complex<double> >;
+template class SparseMat<double>;
+template class SparseMat<complex<double> >;
 
 
 
@@ -790,6 +791,7 @@ while(l)
 return *this;
 }
 
+template<>
 const double SparseMat<double>::dot(const NRMat<double> &rhs) const
 {
 double r=0;
@@ -803,12 +805,14 @@ while(l)
 return r;
 }
 
+template<>
 void NRMat<complex<double> >::gemm(const complex<double>  &beta, const SparseMat<complex<double> > &a, const char transa, const NRMat<complex<double> > &b, const char transb, const complex<double>  &alpha)
 {
 laerror("not implemented yet");
 }
 
 
+template<>
 void NRMat<double>::gemm(const double &beta, const SparseMat<double> &a, const char transa, const NRMat<double> &b, const char transb, const double &alpha)
 {
 bool transpa = tolower(transa)!='n'; //not OK for complex
diff --git a/vec.cc b/vec.cc
index c115c04..80c824c 100644
--- a/vec.cc
+++ b/vec.cc
@@ -29,11 +29,11 @@ INSTANTIZE(short)
 INSTANTIZE(unsigned short)
 INSTANTIZE(char)
 INSTANTIZE(unsigned char)
-template NRVec<double>;
-template NRVec<complex<double> >;
-template NRVec<char>;
-template NRVec<short>;
-template NRVec<int>;
+template class NRVec<double>;
+template class NRVec<complex<double> >;
+template class NRVec<char>;
+template class NRVec<short>;
+template class NRVec<int>;