diff --git a/auxstorage.h b/auxstorage.h
index 04fc3ac..0ea219b 100644
--- a/auxstorage.h
+++ b/auxstorage.h
@@ -17,6 +17,7 @@
 */
 #ifndef _AUXSTORAGE_H_
 #define _AUXSTORAGE_H_
+
 #include "vec.h"
 #include "mat.h"
 #include "smat.h"
@@ -27,6 +28,8 @@
 #include <fcntl.h>
 #include <unistd.h>
 
+namespace LA {
+
 //CAUTION:
 //it will not work if T is itself a class with dynamically allocated components
 //it cannot be implemented for SparseMat, which lacks fixed record length
@@ -136,6 +139,6 @@ if(0>write(fd,&x(0,0),recl)) {perror(""); laerror("write failed in AuxStorage");
 }
 
 
-
+}//namespace
 
 #endif
diff --git a/bisection.h b/bisection.h
index b3ae9e7..349e31d 100644
--- a/bisection.h
+++ b/bisection.h
@@ -17,6 +17,9 @@
 */
 #ifndef _BISECTION_H
 #define _BISECTION_H
+
+namespace LA {
+
 //general bisection search between dm and hm
 //returns dm-1 on failure, otherwise number between dm and hm
 //cmp returns 0 on equal, >0 if first > second argument
@@ -75,5 +78,5 @@ int interpolation_find(INDEX dm0, INDEX high, const SUBJECT *x, const SUBJECT *b
 }
 
 
-
+}//namespace
 #endif
diff --git a/bitvector.cc b/bitvector.cc
index 9c19a5f..d61c0b9 100644
--- a/bitvector.cc
+++ b/bitvector.cc
@@ -18,14 +18,16 @@
 
 #include "bitvector.h"
 
+namespace LA {
+
 //inefficient I/O operators
-ostream & operator<<(ostream &s, const bitvector &x)
+std::ostream & operator<<(std::ostream &s, const bitvector &x)
 {
 for(unsigned int i=0; i<x.size(); ++i) s<< x[i];
 return s;
 }
 
-istream & operator>>(istream  &s, bitvector &x)
+std::istream & operator>>(std::istream  &s, bitvector &x)
 {
 bool a;
 x.copyonwrite();
@@ -168,3 +170,4 @@ if(modulo)
 return s+word_popul(a);
 }
 
+}//namespace
diff --git a/bitvector.h b/bitvector.h
index d930d3b..bdcfa9c 100644
--- a/bitvector.h
+++ b/bitvector.h
@@ -21,6 +21,8 @@
 
 #include "vec.h"
 
+namespace LA {
+
 //compressed storage of large bit vectors
 //any reasonable compiler changes the dividions and modulos to shift/and instructions
 
@@ -75,7 +77,8 @@ public:
 	};
 
 
-extern ostream & operator<<(ostream &s, const bitvector &x);
-extern istream & operator>>(istream  &s, bitvector &x);
+extern std::ostream & operator<<(std::ostream &s, const bitvector &x);
+extern std::istream & operator>>(std::istream  &s, bitvector &x);
 
+}//namespace
 #endif
diff --git a/conjgrad.h b/conjgrad.h
index 4264384..56b27a6 100644
--- a/conjgrad.h
+++ b/conjgrad.h
@@ -15,6 +15,8 @@
     You should have received a copy of the GNU General Public License
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
+#ifndef _CONJGRAD_H_
+#define _CONJGRAD_H_
 #include "vec.h"
 #include "smat.h"
 #include "mat.h"
@@ -22,6 +24,8 @@
 #include "nonclass.h"
 #include <iomanip>
 
+namespace LA {
+
 //conjugate gradient solution of a linear system
 
 //matrix can be any class which has nrows(), ncols(), diagonalof() and gemv() available
@@ -63,10 +67,10 @@ for(int iter=0; iter<= itmax; iter++)
 	double err=p.norm();
 	if(verbose) 
 		{
-		cout << "conjgrad: iter= "<<iter<<" err= "<<
-		setiosflags(ios::scientific)<<setprecision(8) <<err<<
-		resetiosflags(ios::scientific)<<setprecision(12)<<"\n";
-		cout.flush();
+		std::cout << "conjgrad: iter= "<<iter<<" err= "<<
+		std::setiosflags(std::ios::scientific)<<std::setprecision(8) <<err<<
+		std::resetiosflags(std::ios::scientific)<<std::setprecision(12)<<"\n";
+		std::cout.flush();
 		}
 	if(err <= tol)
 		{
@@ -94,3 +98,5 @@ if(!issquare) delete r;
 return false;
 }
 
+}//namespace
+#endif
diff --git a/davidson.h b/davidson.h
index 20c0a36..0bcb06e 100644
--- a/davidson.h
+++ b/davidson.h
@@ -24,6 +24,8 @@
 #include "nonclass.h"
 #include "auxstorage.h"
 
+namespace LA {
+
 //Davidson diagonalization of real symmetric matrix (modified Lanczos)
 
 //matrix can be any class which has nrows(), ncols(), diagonalof(), issymmetric(), and gemv() available
@@ -146,21 +148,21 @@ else
 
 T eival_n=r[nroot];
 oldnroot=nroot;
-typename LA_traits<T>::normtype test=abs(smallV(krylovsize,nroot));
+typename LA_traits<T>::normtype test=std::abs(smallV(krylovsize,nroot));
 if(test<eps) nroot++;
 if(it==0) nroot = 0;
-for(int iroot=0; iroot<=min(krylovsize,nroots-1); ++iroot)
+for(int iroot=0; iroot<=std::min(krylovsize,nroots-1); ++iroot)
 	{
-        test = abs(smallV(krylovsize,iroot));
-        if(test>eps) nroot=min(nroot,iroot);
-        if(verbose && iroot<=max(oldnroot,nroot)) 
+        test = std::abs(smallV(krylovsize,iroot));
+        if(test>eps) nroot=std::min(nroot,iroot);
+        if(verbose && iroot<=std::max(oldnroot,nroot)) 
 		{
-		cout <<"Davidson: iter="<<it <<" dim="<<krylovsize<<" root="<<iroot<<" energy="<<r[iroot]<<"\n";
-		cout.flush();
+		std::cout <<"Davidson: iter="<<it <<" dim="<<krylovsize<<" root="<<iroot<<" energy="<<r[iroot]<<"\n";
+		std::cout.flush();
 		}
 	}
 
-if(verbose && oldnroot!=nroot) {cout <<"root no. "<<oldnroot<<" converged\n"; cout.flush();}
+if(verbose && oldnroot!=nroot) {std::cout <<"root no. "<<oldnroot<<" converged\n"; std::cout.flush();}
 if (nroot>=nroots) goto converged;
 if (it==maxit-1) break; //not converged
 
@@ -195,7 +197,7 @@ eival_n = r[nroot];
 for(i=0; i<n; ++i)
 	{
 	T denom = diagonal[i] - eival_n;
-        denom = denom<0?-max(0.1,abs(denom)):max(0.1,abs(denom));
+        denom = denom<0?-std::max(0.1,std::abs(denom)):std::max(0.1,std::abs(denom));
         vec1[i] /= denom;
 	}
 
@@ -224,7 +226,7 @@ goto finished;
 converged:
 AuxStorage<typename LA_traits<T>::elementtype> *ev;
 if(eivecsfile) ev = new AuxStorage<typename LA_traits<T>::elementtype>(eivecsfile);
-if(verbose) {cout << "Davidson converged in "<<it<<" iterations.\n"; cout.flush();}
+if(verbose) {std::cout << "Davidson converged in "<<it<<" iterations.\n"; std::cout.flush();}
 for(nroot=0; nroot<nroots; ++nroot)
 	{
         eivals[nroot]=r[nroot];
@@ -254,6 +256,5 @@ else  {delete s0; delete s1;}
 if(flag) laerror("no convergence in davidson");
 }
 
-
-
+}//namespace
 #endif
diff --git a/diis.h b/diis.h
index 9d015b3..2227605 100644
--- a/diis.h
+++ b/diis.h
@@ -26,8 +26,10 @@
 #include "la_traits.h"
 #include "auxstorage.h"
 
+namespace LA {
+
 //Pulay memorial book  remarks - for numerical stabilization small addition to diagonal
-#define DIISEPS 1e-9 
+const double DIISEPS=1e-9;
 
 // Typically, T is some solution vector in form of NRVec, NRMat, or NRSMat over double or complex<double> fields
 // actually it can be anything what has operator=(const T&), clear(), dot() , axpy(), norm() and copyonwrite(), and LA_traits<T>::normtype and elementtype
@@ -150,7 +152,7 @@ rhs= (Ue)0; rhs[0]= (Ue)-1;
 NRSMat<Te> amat=bmat;
 linear_solve(amat,rhs,NULL,aktdim+1);
 }
-if(verbose) cout <<"DIIS coefficients: "<<rhs<<endl;
+if(verbose) std::cout <<"DIIS coefficients: "<<rhs<<std::endl;
 
 //build the new linear combination
 vec.clear();
@@ -170,5 +172,6 @@ else
 return norm;
 }
 
+}//namespace
 
 #endif
diff --git a/fourindex.h b/fourindex.h
index c9d727c..6da9bde 100644
--- a/fourindex.h
+++ b/fourindex.h
@@ -26,9 +26,13 @@
 #include <sys/stat.h>
 #include <unistd.h>
 #include <sys/stat.h>
-#include "la.h"
 #include "laerror.h"
+#include "vec.h"
+#include "smat.h"
+#include "mat.h"
+#include "nonclass.h"
 
+namespace LA {
 
 static unsigned int hcd0(unsigned int big,unsigned int small)
 {
@@ -591,7 +595,7 @@ return n;
 }
 
 template <class I, class T>
-ostream& operator<<(ostream &s, const fourindex_ext<I,T> &x)
+std::ostream& operator<<(std::ostream &s, const fourindex_ext<I,T> &x)
                 {
                 int n;
                 n=x.size();
@@ -609,7 +613,7 @@ ostream& operator<<(ostream &s, const fourindex_ext<I,T> &x)
 
 
 template <class I, class T>
-ostream& operator<<(ostream &s, const fourindex<I,T> &x)
+std::ostream& operator<<(std::ostream &s, const fourindex<I,T> &x)
                 {
                 int n;
                 n=x.size();
@@ -625,7 +629,7 @@ ostream& operator<<(ostream &s, const fourindex<I,T> &x)
                 }
 
 template <class I, class T>
-istream& operator>>(istream  &s, fourindex<I,T> &x)
+std::istream& operator>>(std::istream  &s, fourindex<I,T> &x)
                 {
                 typename LA_traits_io<I>::IOtype i,j,k,l;
 		typename LA_traits_io<T>::IOtype elem;
@@ -680,7 +684,7 @@ public:
 	const T& operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const;
 	void resize(const int n) {(*this).NRSMat<T>::resize(n*(n+1)/2);};
 	void putext(int f, T thr=1e-15);
-	int nbas() const {return (int)sqrt(2*(*this).nrows());};
+	int nbas() const {return (int)std::sqrt(2*(*this).nrows());};
 
 };
 
@@ -795,10 +799,10 @@ if (!NRMat<T>::v) laerror("access to unallocated fourindex_dense");
 return (*this).NRMat<T>::operator() ((j-1)*noca+i-1,(b-1)*nvra+a-1);
 }
 
-	void print(ostream &out) const
+	void print(std::ostream &out) const
 		{
 		unsigned int i,j,a,b;
-		for(i=1; i<=noca; ++i) for(j=1; j<=nocb; ++j) for(a=1; a<=nvra; ++a) for(b=1; b<=nvrb; ++b) out << i<<" "<<j<<" "<<a<<" "<<b<<" "<<(*this)(i,j,a,b)<<endl;
+		for(i=1; i<=noca; ++i) for(j=1; j<=nocb; ++j) for(a=1; a<=nvra; ++a) for(b=1; b<=nvrb; ++b) out << i<<" "<<j<<" "<<a<<" "<<b<<" "<<(*this)(i,j,a,b)<<std::endl;
 		}
 };
 
@@ -854,10 +858,10 @@ if(a<b) {minus++; unsigned int t=a; a=b; b=t;}
 }
 
 
-        void print(ostream &out) const
+        void print(std::ostream &out) const
                 {
                 unsigned int i,j,a,b;
-                for(i=1; i<=nocc; ++i) for(j=1; j<i; ++j) for(a=1; a<=nvrt; ++a) for(b=1; b<a; ++b) out << i<<" "<<j<<" "<<a<<" "<<b<<" "<<(*this)(i,j,a,b)<<endl;
+                for(i=1; i<=nocc; ++i) for(j=1; j<i; ++j) for(a=1; a<=nvrt; ++a) for(b=1; b<a; ++b) out << i<<" "<<j<<" "<<a<<" "<<b<<" "<<(*this)(i,j,a,b)<<std::endl;
                 }
 
 
@@ -886,14 +890,14 @@ public:
 	void add_unique(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem);
 	const T& operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const;
 	void resize(const int n) {nbas=n; (*this).NRSMat<T>::resize(n*(n-1)/2);};
-        void print(ostream &out) const
+        void print(std::ostream &out) const
                 {
                 unsigned int i,j,k,l;
 		for(i=1; i<=nbas; ++i)
        			 for(k=1;k<i; ++k)
                 		for(j=1; j<=i; ++j)
                         		for(l=1; l<j && (j==i ? l<=k : 1); ++l)
-                				cout << i<<" "<<k<<" "<<j<<" "<<l<<" "<<(*this)(i,k,j,l)<<endl;
+                				std::cout << i<<" "<<k<<" "<<j<<" "<<l<<" "<<(*this)(i,k,j,l)<<std::endl;
                 }
 
 
@@ -992,6 +996,6 @@ for(p= const_cast<fourindex_ext<I,T> *>(&rhs)->pbegin(); p.notend(); ++p)
 }
 
 
-
+}//namespace
 
 #endif /*_fourindex_included*/
diff --git a/gmres.h b/gmres.h
index e6eb1fb..91dbf0f 100644
--- a/gmres.h
+++ b/gmres.h
@@ -16,6 +16,8 @@
     You should have received a copy of the GNU General Public License
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
+#ifndef _GMRES_H
+#define _GMRES_H
 #include "vec.h"
 #include "smat.h"
 #include "mat.h"
@@ -24,6 +26,8 @@
 #include <iomanip>
 #include "auxstorage.h"
 
+namespace LA {
+
 //GMRES solution of a linear system
 
 //matrix can be any class which has nrows(), ncols(), diagonalof() and gemv() available
@@ -135,7 +139,7 @@ for (int l=0;l<neustart;l++)  // main loop for restarts
 	    	H(k+1,k) = tmp= z.norm();
 	    	if(tmp < 1.e-2*eps )
 			{
-	    		if(verbose) cerr <<("gmres restart performed\n");
+	    		if(verbose) std::cerr <<("gmres restart performed\n");
 	       		// Abbruchbedingung, konstruiere x_k
 	       		for (int i=0;i<k;i++) 
 				{
@@ -190,12 +194,12 @@ for (int l=0;l<neustart;l++)  // main loop for restarts
 	    	//Schritt 6: Konvergenz?
 	    	if(verbose) 
 			{
-			cout << "gmres iter "<<l<<" "<<k<<" resid "
-		<<setw(0)<<setiosflags(ios::scientific)<<setprecision(8)
-		<<abs(d[k+1])<< " thr "<<eps*beta_0<< " reduction "
-		<<setw(5)<<setprecision(2)<<resetiosflags(ios::scientific)
-		<<(d_alt - abs(d[k+1]))/d_alt*100<< "\n" <<setprecision(12);
-			cout.flush();
+			std::cout << "gmres iter "<<l<<" "<<k<<" resid "
+		<<std::setw(0)<<std::setiosflags(std::ios::scientific)<<std::setprecision(8)
+		<<std::abs(d[k+1])<< " thr "<<eps*beta_0<< " reduction "
+		<<std::setw(5)<<std::setprecision(2)<<std::resetiosflags(std::ios::scientific)
+		<<(d_alt - std::abs(d[k+1]))/d_alt*100<< "\n" <<std::setprecision(12);
+			std::cout.flush();
 			}
 		
 		d_alt = abs(d[k+1]);
@@ -258,3 +262,6 @@ if(!incore) delete st;
 return !flag;
 }
 
+}//namespace
+
+#endif
diff --git a/la.h b/la.h
index afa431c..179e066 100644
--- a/la.h
+++ b/la.h
@@ -20,10 +20,29 @@
 
 //this should be the single include file for the end user
 //
-
-#include "vec.h"
-#include "smat.h"
+#ifdef USE_TRACEBACK
+#include "traceback.h"
+#endif
+#include "la_traits.h"
+#include "laerror.h"
+#include "auxstorage.h"
+#include "bisection.h"
+#include "bitvector.h"
+#include "conjgrad.h"
+#include "davidson.h"
+#include "diis.h"
+#include "fourindex.h"
+#include "gmres.h"
 #include "mat.h"
+#include "matexp.h"
+#include "noncblas.h"
 #include "nonclass.h"
+#include "permutation.h"
+#include "qsort.h"
+#include "smat.h"
+#include "sparsemat.h"
+#include "sparsesmat.h"
+#include "vec.h"
+
 
 #endif /* _LA_H_ */
diff --git a/la_traits.h b/la_traits.h
index 2330ecb..7911498 100644
--- a/la_traits.h
+++ b/la_traits.h
@@ -27,14 +27,17 @@
 #define _LA_TRAITS_INCL
 
 
-extern bool _LA_count_check;
 
 
-using namespace std;
 #include <stdio.h>
 #include <string.h>
 #include <iostream>
 #include <complex>
+
+
+//using namespace std;
+#define complex std::complex
+
 #include "laerror.h"
 
 #ifdef NONCBLAS
@@ -53,6 +56,9 @@ extern "C" {
 }
 #endif
 
+namespace LA {
+
+extern bool _LA_count_check;
 
 //forward declarations
 template<typename C> class NRVec;
@@ -61,6 +67,7 @@ template<typename C> class NRMat_from1;
 template<typename C> class NRSMat;
 template<typename C> class NRSMat_from1;
 template<typename C> class SparseMat;
+template<typename C> class SparseSMat;
 
 
 typedef class {} Dummy_type;
@@ -95,9 +102,13 @@ SPECIALIZE_COMPLEX(complex<float>)
 SPECIALIZE_COMPLEX(char)
 SPECIALIZE_COMPLEX(unsigned char)
 SPECIALIZE_COMPLEX(short)
+SPECIALIZE_COMPLEX(unsigned short)
 SPECIALIZE_COMPLEX(int)
 SPECIALIZE_COMPLEX(unsigned int)
+SPECIALIZE_COMPLEX(long)
 SPECIALIZE_COMPLEX(unsigned long)
+SPECIALIZE_COMPLEX(long long)
+SPECIALIZE_COMPLEX(unsigned long long)
 
 
 //for general sortable classes
@@ -194,10 +205,11 @@ struct LA_traits_aux<complex<C>, scalar_true> {
 typedef complex<C> elementtype;
 typedef complex<C> producttype;
 typedef C normtype;
+static inline C sqrabs(const complex<C> x) { return x.real()*x.real()+x.imag()*x.imag();}
 static inline bool gencmp(const complex<C> *x, const complex<C> *y, int n) {return memcmp(x,y,n*sizeof(complex<C>));}
 static bool bigger(const  complex<C> &x, const complex<C> &y) {laerror("complex comparison undefined"); return false;}
 static bool smaller(const  complex<C> &x, const complex<C> &y) {laerror("complex comparison undefined"); return false;}
-static inline normtype norm (const  complex<C> &x) {return abs(x);}
+static inline normtype norm (const  complex<C> &x) {return std::abs(x);}
 static inline void axpy (complex<C> &s, const complex<C> &x, const complex<C> &c) {s+=x*c;}
 static inline void get(int fd, complex<C> &x, bool dimensions=0, bool transp=0) {if(sizeof(complex<C>)!=read(fd,&x,sizeof(complex<C>))) laerror("read error");}
 static inline void put(int fd, const complex<C> &x, bool dimensions=0, bool transp=0) {if(sizeof(complex<C>)!=write(fd,&x,sizeof(complex<C>))) laerror("write error");}
@@ -206,6 +218,7 @@ static void multiput(unsigned int n, int fd, const complex<C> *x, bool dimension
 static void copy(complex<C> *dest, complex<C> *src, unsigned int n) {memcpy(dest,src,n*sizeof(complex<C>));}
 static void clear(complex<C> *dest, unsigned int n) {memset(dest,0,n*sizeof(complex<C>));}
 static void copyonwrite(complex<C> &x) {};
+static void clearme(complex<C> &x) {x=0;};
 };
 
 //non-complex scalars
@@ -214,10 +227,11 @@ struct LA_traits_aux<C, scalar_true> {
 typedef C elementtype;
 typedef C producttype;
 typedef C normtype;
+static inline C sqrabs(const C x) { return x*x;}
 static inline bool gencmp(const C *x, const C *y, int n) {return memcmp(x,y,n*sizeof(C));}
 static inline bool bigger(const  C &x, const C &y) {return x>y;}
 static inline bool smaller(const  C &x, const C &y) {return x<y;}
-static inline normtype norm (const  C &x) {return abs(x);}
+static inline normtype norm (const  C &x) {return std::abs(x);}
 static inline void axpy (C &s, const C &x, const C &c) {s+=x*c;}
 static inline void put(int fd, const C &x, bool dimensions=0, bool transp=0) {if(sizeof(C)!=write(fd,&x,sizeof(C))) laerror("write error");}
 static inline void get(int fd, C &x, bool dimensions=0, bool transp=0) {if(sizeof(C)!=read(fd,&x,sizeof(C))) laerror("read error");}
@@ -226,6 +240,7 @@ static void multiget(unsigned int n, int fd, C *x, bool dimensions=0) {if((ssize
 static void copy(C *dest, C *src, unsigned int n) {memcpy(dest,src,n*sizeof(C));}
 static void clear(C *dest, unsigned int n) {memset(dest,0,n*sizeof(C));}
 static void copyonwrite(C &x) {};
+static void clearme(complex<C> &x) {x=0;};
 };
 
 
@@ -252,6 +267,7 @@ static void multiget(unsigned int n,int fd, C *x, bool dimensions=1) {for(unsign
 static void copy(C *dest, C *src, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i]=src[i];} \
 static void clear(C *dest, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i].clear();}\
 static void copyonwrite(X<C> &x) {x.copyonwrite();}\
+static void clearme(X<C> &x) {x.clear();}\
 };
 
 
@@ -260,6 +276,7 @@ generate_traits(NRMat)
 generate_traits(NRMat_from1)
 generate_traits(NRVec)
 generate_traits(SparseMat)
+generate_traits(SparseSMat) //product leading to non-symmetric result not implemented
 
 #undef generate_traits
 
@@ -282,6 +299,7 @@ static void multiget(unsigned int n,int fd, C *x, bool dimensions=1) {for(unsign
 static void copy(C *dest, C *src, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i]=src[i];}  \
 static void clear(C *dest, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i].clear();} \
 static void copyonwrite(X<C> &x) {x.copyonwrite();} \
+static void clearme(X<C> &x) {x.clear();} \
 };
 
 generate_traits_smat(NRSMat)
@@ -292,4 +310,6 @@ generate_traits_smat(NRSMat_from1)
 template<typename C>
 struct LA_traits : LA_traits_aux<C, typename isscalar<C>::scalar_type> {};
 
+}//namespace
+
 #endif
diff --git a/laerror.cc b/laerror.cc
index ba130c1..219ecba 100644
--- a/laerror.cc
+++ b/laerror.cc
@@ -30,6 +30,8 @@
 #include "traceback.h"
 #endif
 
+namespace LA {
+
 bool _LA_count_check=true;
 
 extern "C" void _findme(void) {}; //for autoconf test we need a function with C linkage
@@ -86,3 +88,5 @@ int cblas_errprn(int ierr, int info, char *form, ...)
    laerror(msg0);
 return 0;
 }
+
+}//namespace
diff --git a/laerror.h b/laerror.h
index 208f2c1..1dbefea 100644
--- a/laerror.h
+++ b/laerror.h
@@ -17,9 +17,10 @@
 */
 #ifndef _LAERROR_H_
 #define _LAERROR_H_
-
 #include <iostream>
 
+namespace LA {
+
 //exception class for laerror
 class LAerror
 	{
@@ -36,5 +37,6 @@ s << x.msg;
 return s;
 }
 
+}//namespace
 
 #endif
diff --git a/mat.cc b/mat.cc
index 013f89c..9375348 100644
--- a/mat.cc
+++ b/mat.cc
@@ -31,6 +31,7 @@ extern ssize_t write(int, const void *, size_t);
 // TODO :
 //
 
+namespace LA {
 
 /*
  *  Templates first, specializations for BLAS next
@@ -439,9 +440,9 @@ NRSMat<complex<double> > r(mm,mm);
 int i,j;
 for(i=0; i<mm; ++i) for(j=0; j<=i; ++j)
 #ifdef MATPTR
-	cblas_zdotc_sub(nn, v[0]+i , mm,v[0]+j, mm, (void *)(&r(i,j)));
+	cblas_zdotc_sub(nn, v[0]+i , mm,v[0]+j, mm, &r(i,j));
 #else
-	cblas_zdotc_sub(nn, v+i , mm,v+j, mm, (void *)(&r(i,j)));
+	cblas_zdotc_sub(nn, v+i , mm,v+j, mm, &r(i,j));
 #endif
 return r;
 }
@@ -486,9 +487,9 @@ NRSMat<complex<double> > r(nn,nn);
 int i,j;
 for(i=0; i<nn; ++i) for(j=0; j<=i; ++j)
 #ifdef MATPTR
-	cblas_zdotc_sub(mm, v[i] , 1,v[j], 1, (void *)(&r(i,j)));
+	cblas_zdotc_sub(mm, v[i] , 1,v[j], 1, &r(i,j));
 #else
-	cblas_zdotc_sub(mm, v+i*mm , 1,v+j*mm, 1, (void *)(&r(i,j)));
+	cblas_zdotc_sub(mm, v+i*mm , 1,v+j*mm, 1, &r(i,j));
 #endif
 return r;
 }
@@ -548,7 +549,7 @@ NRMat< complex<double> > &
 NRMat< complex<double> >::operator*=(const complex<double> &a)
 {
 	copyonwrite();
-	cblas_zscal(nn*mm, &a, (void *)(*this)[0], 1);
+	cblas_zscal(nn*mm, &a, (*this)[0], 1);
 	return *this;
 }
 
@@ -591,7 +592,7 @@ NRMat< complex<double> >::operator+=(const NRMat< complex<double> >  &rhs)
 		laerror("Mat += Mat of incompatible matrices");
 #endif
 	 copyonwrite();
-	 cblas_zaxpy(nn*mm, &CONE, (void *)rhs[0], 1, (void *)(*this)[0], 1);
+	 cblas_zaxpy(nn*mm, &CONE, rhs[0], 1, (*this)[0], 1);
 	 return *this;
 }
 
@@ -639,7 +640,7 @@ NRMat< complex<double> >::operator-=(const NRMat< complex<double> >  &rhs)
 		laerror("Mat -= Mat of incompatible matrices");
 #endif
 	 copyonwrite();
-	 cblas_zaxpy(nn*mm, &CMONE, (void *)rhs[0], 1, (void *)(*this)[0], 1);
+	 cblas_zaxpy(nn*mm, &CMONE, rhs[0], 1, (*this)[0], 1);
 	 return *this;
 }
 
@@ -696,12 +697,12 @@ NRMat< complex<double> >::operator+=(const NRSMat< complex<double> > &rhs)
 	const complex<double> *p = rhs;
 	copyonwrite();
 	for (int i=0; i<nn; i++) {
-		cblas_zaxpy(i+1, (void *)&CONE, (void *)p, 1, (void *)(*this)[i], 1); 
+		cblas_zaxpy(i+1, &CONE, p, 1, (*this)[i], 1); 
 		p += i+1;
 	}
 	p = rhs; p++;
 	for (int i=1; i<nn; i++) {
-		cblas_zaxpy(i, (void *)&CONE, (void *)p, 1, (void *)((*this)[i]+i), nn);
+		cblas_zaxpy(i, &CONE, p, 1, (*this)[0]+i, nn);
 		p += i+1;
 	}
 	return *this;
@@ -725,7 +726,7 @@ NRMat<T> & NRMat<T>::operator+=(const NRSMat<T> &rhs)
         }
         p = rhs; p++;
         for (int i=1; i<nn; i++) {
-		for(int j=0; j<i; ++j) *((*this)[i]+i+nn*j) += p[j];
+		for(int j=0; j<i; ++j) *((*this)[j]+i) += p[j];
                 p += i+1;
         }
         return *this;
@@ -765,12 +766,12 @@ NRMat< complex<double> >::operator-=(const NRSMat< complex<double> > &rhs)
 	const complex<double> *p = rhs;
 	copyonwrite();
 	for (int i=0; i<nn; i++) {
-		cblas_zaxpy(i+1, (void *)&CMONE, (void *)p, 1, (void *)(*this)[i], 1); 
+		cblas_zaxpy(i+1, &CMONE, p, 1, (*this)[i], 1); 
 		p += i+1;
 	}
 	p = rhs; p++;
 	for (int i=1; i<nn; i++) {
-		cblas_zaxpy(i, (void *)&CMONE, (void *)p, 1, (void *)((*this)[i]+i), nn);
+		cblas_zaxpy(i, &CMONE, p, 1, (*this)[0]+i, nn);
 		p += i+1;
 	}
 	return *this;
@@ -792,7 +793,7 @@ NRMat<T> & NRMat<T>::operator-=(const NRSMat<T> &rhs)
         }
         p = rhs; p++;
         for (int i=1; i<nn; i++) {
-		for(int j=0; j<i; ++j) *((*this)[i]+i+nn*j) -= p[j];
+		for(int j=0; j<i; ++j) *((*this)[j]+i) -= p[j];
                 p += i+1;
         }
         return *this;
@@ -821,8 +822,8 @@ NRMat< complex<double> >::dot(const NRMat< complex<double> > &rhs) const
 	if(nn!=rhs.nn || mm!= rhs.mm) laerror("Mat.Mat incompatible matrices");
 #endif
 	complex<double> dot;
-	cblas_zdotc_sub(nn*mm, (void *)(*this)[0], 1, (void *)rhs[0], 1, 
-			(void *)(&dot));
+	cblas_zdotc_sub(nn*mm, (*this)[0], 1, rhs[0], 1, 
+			&dot);
 	return dot;
 }
 
@@ -851,8 +852,8 @@ NRMat< complex<double> >::operator*(const NRMat< complex<double> > &rhs) const
 #endif
 	NRMat< complex<double> > result(nn, rhs.mm);
 	cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nn, rhs.mm, mm,
-			(const void *)(&CONE),(const void *)(*this)[0], mm, (const void *)rhs[0],
-			rhs.mm, (const void *)(&CZERO), (void *)result[0], rhs.mm);
+			&CONE,(*this)[0], mm, rhs[0],
+			rhs.mm, &CZERO, result[0], rhs.mm);
 	return result;
 }
 
@@ -936,8 +937,8 @@ NRMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
 #endif
 	NRMat< complex<double> > result(nn, rhs.ncols());
 	for (int i=0; i<nn; i++)
-		cblas_zhpmv(CblasRowMajor, CblasLower, mm, (void *)&CONE, (void *)&rhs[0],
-				(void *)(*this)[i], 1, (void *)&CZERO, (void *)result[i], 1);
+		cblas_zhpmv(CblasRowMajor, CblasLower, mm, &CONE, &rhs[0],
+				(*this)[i], 1, &CZERO, result[i], 1);
 	return result;
 }
 
@@ -969,7 +970,7 @@ NRMat< complex<double> >::transpose(bool conj) const
 {
 	NRMat< complex<double> > result(mm,nn);
 	for (int i=0; i<nn; i++) 
-		cblas_zcopy(mm, (void *)(*this)[i], 1, (void *)(result[0]+i), nn);
+		cblas_zcopy(mm, (*this)[i], 1, (result[0]+i), nn);
 	if (conj) cblas_dscal(mm*nn, -1.0, (double *)(result[0])+1, 2);
 	return result;
 }
@@ -1039,7 +1040,7 @@ const double  NRMat<double>::norm(const double scalar) const
 			if (i==j) tmp -= scalar;
 			sum += tmp*tmp;
 		}
-	return sqrt(sum);
+	return std::sqrt(sum);
 }
 
 
@@ -1060,7 +1061,7 @@ const double  NRMat< complex<double> >::norm(const complex<double> scalar) const
 			if (i==j) tmp -= scalar;
 			sum += tmp.real()*tmp.real()+tmp.imag()*tmp.imag();
 		}
-	return sqrt(sum);
+	return std::sqrt(sum);
 }
 
 
@@ -1087,7 +1088,7 @@ void NRMat< complex<double> >::axpy(const complex<double> alpha,
 	if (nn!=mat.nn || mm!=mat.mm) laerror("zaxpy of incompatible matrices");
 #endif
 	copyonwrite();
-	cblas_zaxpy(nn*mm, (void *)&alpha, mat, 1, (void *)(*this)[0], 1);
+	cblas_zaxpy(nn*mm, &alpha, mat, 1, (*this)[0], 1);
 }
 
 
@@ -1186,7 +1187,7 @@ if(rowcol) //vectors are rows
 		NRVec<double> tmp = *metric * (*this).row(j);
                 double norm = cblas_ddot(mm,(*this)[j],1,tmp,1);
                 if(norm<=0.) laerror("zero vector in orthonormalize or nonpositive metric");
-                cblas_dscal(mm,1./sqrt(norm),(*this)[j],1);
+                cblas_dscal(mm,1./std::sqrt(norm),(*this)[j],1);
                 }
 	}
 else	//vectors are columns
@@ -1203,7 +1204,7 @@ else	//vectors are columns
 		NRVec<double> tmp = *metric * (*this).column(j);
 		double norm = cblas_ddot(nn,&(*this)[0][j],mm,tmp,1);
                 if(norm<=0.) laerror("zero vector in orthonormalize or nonpositive metric");
-                cblas_dscal(nn,1./sqrt(norm),&(*this)[0][j],mm);
+                cblas_dscal(nn,1./std::sqrt(norm),&(*this)[0][j],mm);
                 }
 	}
 }
@@ -1261,3 +1262,4 @@ template class NRMat<unsigned int>;
 template class NRMat<unsigned long>;
 template class NRMat<unsigned long long>;
 
+}//namespace
diff --git a/mat.h b/mat.h
index a1ce209..9ed2abd 100644
--- a/mat.h
+++ b/mat.h
@@ -19,9 +19,9 @@
 */
 #ifndef _LA_MAT_H_
 #define _LA_MAT_H_
-
 #include "la_traits.h"
 
+namespace LA {
 template <typename T>
 class NRMat {
 protected:
@@ -91,6 +91,7 @@ public:
 	const NRVec<complex<T> > operator*(const NRVec<complex<T> > &rhs) const {NRVec<complex<T> > result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec
 	const NRVec<T> rsum() const; //sum of rows
 	const NRVec<T> csum() const; //sum of columns
+	void orthonormalize(const bool rowcol, const NRSMat<T> *metric=NULL);//orthonormalize (true - vectors are rows)
 	const NRVec<T> row(const int i, int l= -1) const; //row of, efficient
 	const NRVec<T> column(const int j, int l= -1) const {if(l<0) l=nn; NRVec<T> r(l); for(int i=0; i<l; ++i) r[i]= (*this)(i,j); return r;}; //column of, general but not very efficient
 	const T* diagonalof(NRVec<T> &, const bool divide=0, bool cache=false) const; //get diagonal 
@@ -140,11 +141,13 @@ public:
 #endif
 };
 
+}//namespace
 //due to mutual includes this has to be after full class declaration
 #include "vec.h"
 #include "smat.h"
 #include "sparsemat.h"
 
+namespace LA {
 // ctors
 template <typename T>
 NRMat<T>::NRMat(const int n, const int m) : nn(n), mm(m), count(new int)
@@ -374,9 +377,9 @@ template<>
 inline const complex<double>  NRMat< complex<double> >::amax() const
 {
 #ifdef MATPTR
-	return v[0][cblas_izamax(nn*mm, (void *)v[0], 1)];
+	return v[0][cblas_izamax(nn*mm, v[0], 1)];
 #else
-	return v[cblas_izamax(nn*mm, (void *)v, 1)];
+	return v[cblas_izamax(nn*mm, v, 1)];
 #endif
 }
 
@@ -576,7 +579,7 @@ return r;
 
 // I/O
 template <typename T>
-ostream& operator<<(ostream &s, const NRMat<T> &x)
+std::ostream& operator<<(std::ostream &s, const NRMat<T> &x)
                 {
                 int i,j,n,m;
                 n=x.nrows();
@@ -590,7 +593,7 @@ ostream& operator<<(ostream &s, const NRMat<T> &x)
                 }
 
 template <typename T>
-istream& operator>>(istream  &s, NRMat<T> &x)
+std::istream& operator>>(std::istream  &s, NRMat<T> &x)
                 {
                 int i,j,n,m;
                 s >> n >> m;
@@ -655,4 +658,5 @@ NRVECMAT_OPER(Mat,*)
 NRVECMAT_OPER2(Mat,+)
 NRVECMAT_OPER2(Mat,-)
 
+}//namespace
 #endif /* _LA_MAT_H_ */
diff --git a/matexp.h b/matexp.h
index d7d3c51..9910cb5 100644
--- a/matexp.h
+++ b/matexp.h
@@ -21,10 +21,11 @@
 //of matrix-matrix multiplications on cost of additions and memory
 // the polynom and exp routines will work on any type, for which traits class
 // is defined containing definition of an element type, norm and axpy operation
-
 #include "la_traits.h"
 #include "laerror.h"
+#include <math.h>
 
+namespace LA {
 
 template<class T,class R>
 const T polynom0(const T &x, const NRVec<R> &c)
@@ -40,6 +41,7 @@ else
 	z=x*c[order];
 	for(i=order-1; i>=0; i--)
 		{
+		//std::cerr<<"TEST polynom0 "<<i<<'\n';
 		if(i<order-1) z=y*x;
 		y=z+c[i];
 		}
@@ -138,7 +140,7 @@ for(int i=1; i<=n; ++i)
 	{
 	factor/=i;
 	z= z*t-t*z;
-	if(verbose) cerr << "BCH contribution at order "<<i<<" : "<<z.norm()*factor<<endl;
+	if(verbose) std::cerr << "BCH contribution at order "<<i<<" : "<<z.norm()*factor<<std::endl;
 	result+= z*factor; 
 	}
 return result;
@@ -169,10 +171,10 @@ return y;
 
 inline int nextpow2(const double n)
 {
-const double log2=log(2.);
+const double log2=std::log(2.);
 if(n<=.75) return 0; //try to keep the taylor expansion short
 if(n<=1.) return 1;
-return int(ceil(log(n)/log2-log(.75)));
+return int(std::ceil(std::log(n)/log2-std::log(.75)));
 }
 
 //should better be computed by mathematica to have accurate last digits, perhaps chebyshev instead, see exp in glibc
@@ -208,10 +210,10 @@ template<class T, class C, class S>
 NRVec<C> exp_aux(const T &x, int &power, int maxpower, int maxtaylor, S prescale)
 {
 
-double mnorm= x.norm() * abs(prescale);
+double mnorm= x.norm() * std::abs(prescale);
 power=nextpow2(mnorm); 
 if(maxpower>=0 && power>maxpower) power=maxpower;
-double scale=exp(-log(2.)*power);
+double scale=std::exp(-std::log(2.)*power);
 
 
 //find how long taylor expansion will be necessary
@@ -236,8 +238,8 @@ for(i=0,t=1.;i<=n;i++)
 	taylor2[i]=exptaylor[i]*t;
 	t*=scale;
 	}
-//cout <<"TEST power, scale "<<power<<" "<<scale<<endl; 
-//cout <<"TEST taylor2 "<<taylor2<<endl;
+//std::cout <<"TEST power, scale "<<power<<" "<<scale<<std::endl; 
+//std::cout <<"TEST taylor2 "<<taylor2<<std::endl;
 return taylor2;
 }
 
@@ -246,10 +248,10 @@ return taylor2;
 template<class T, class C, class S>
 void sincos_aux(NRVec<C> &si, NRVec<C> &co, const T &x, int &power,int maxpower, int maxtaylor, const S prescale)
 {
-double mnorm= x.norm() * abs(prescale);
+double mnorm= x.norm() * std::abs(prescale);
 power=nextpow2(mnorm); 
 if(maxpower>=0 && power>maxpower) power=maxpower;
-double scale=exp(-log(2.)*power);
+double scale=std::exp(-std::log(2.)*power);
 
 //find how long taylor expansion will be necessary
 const double precision=1e-14; //further decreasing brings nothing
@@ -275,8 +277,8 @@ for(i=0,t=1.;i<=n;i++)
 	else 	co[i>>1] = exptaylor[i]* (i&2?-t:t);
 	t*=scale;
 	}
-//cout <<"TEST sin "<<si<<endl;
-//cout <<"TEST cos "<<co<<endl;
+//std::cout <<"TEST sin "<<si<<std::endl;
+//std::cout <<"TEST cos "<<co<<std::endl;
 }
 
 
@@ -290,6 +292,7 @@ int power;
 //prepare the polynom of and effectively scale T
 NRVec<typename LA_traits<T>::normtype> taylor2=exp_aux<T,typename LA_traits<T>::normtype,double>(x,power,maxpower,maxtaylor,1.);
 
+//std::cerr <<"TEST power "<<power<<std::endl;
 
 T r= horner?polynom0(x,taylor2):polynom(x,taylor2); 
 //for accuracy summing from the smallest terms up would be better, but this is more efficient for matrices
@@ -309,6 +312,7 @@ int power;
 NRVec<typename LA_traits<T>::normtype> taylors,taylorc;
 sincos_aux<T,typename LA_traits<T>::normtype>(taylors,taylorc,x,power,maxpower,maxtaylor,1.);
 
+
 //could we save something by computing both polynoms simultaneously?
 {
 T x2 = x*x;
@@ -446,16 +450,14 @@ NRVec<typename LA_traits<V>::normtype> taylors,taylorc;
 sincos_aux<M,typename LA_traits<V>::normtype>(taylors,taylorc,mat,power,maxpower,maxtaylor,scale);
 if(taylors.size()!=taylorc.size()) laerror("internal error - same size of sin and cos expansions assumed");
 //the actual computation and resursive "squaring"
-cout <<"TEST power "<<power<<endl;
+//std::cout <<"TEST power "<<power<<std::endl;
 sincostimes_aux(mat,si,co,rhs,taylors,taylorc,transpose,scale,power);
 
 return;
 }
 
 
-
-
-
 //@@@ power series matrix logarithm?
 
+}//namespace
 #endif
diff --git a/noncblas.cc b/noncblas.cc
index 6fcb9d0..00be7d1 100644
--- a/noncblas.cc
+++ b/noncblas.cc
@@ -21,6 +21,8 @@
 #include "laerror.h"
 #include "mat.h"
 
+namespace LA {
+
 #ifdef FORTRAN_
 #define FORNAME(x) x##_
 #else
@@ -255,3 +257,5 @@ return INFO;
 
 #endif
 
+}//namespace
+
diff --git a/nonclass.cc b/nonclass.cc
index 3a5aee3..0dd38e4 100644
--- a/nonclass.cc
+++ b/nonclass.cc
@@ -24,6 +24,7 @@
 #include "nonclass.h"
 #include "qsort.h"
 
+namespace LA {
 
 #ifdef FORTRAN_
 #define FORNAME(x) x##_
@@ -39,9 +40,13 @@ INSTANTIZE(complex<double>)
 INSTANTIZE(int)
 INSTANTIZE(short)
 INSTANTIZE(char)
+INSTANTIZE(long)
+INSTANTIZE(long long)
 INSTANTIZE(unsigned char)
-INSTANTIZE(unsigned long)
+INSTANTIZE(unsigned short)
 INSTANTIZE(unsigned int)
+INSTANTIZE(unsigned long)
+INSTANTIZE(unsigned long long)
 
 #define EPSDET 1e-300
 
@@ -146,7 +151,7 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
 	if (det && r==0) {
 		*det = 1.;
 		//take into account some numerical instabilities in dgesv for singular matrices
-		for (int i=0; i<n; ++i)  {double t=A[i][i]; if(!finite(t) || abs(t) < EPSDET ) {*det=0.; break;} else *det *=t;}
+		for (int i=0; i<n; ++i)  {double t=A[i][i]; if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *=t;}
 		//change sign of det by parity of ipiv permutation
 		if(*det) for (int i=0; i<n; ++i) if(
 #ifdef NONCBLAS
@@ -158,9 +163,9 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
 	}
 	if(det && r>0) *det = 0;
 /*
-	cout <<"ipiv = ";
-	for (int i=0; i<n; ++i) cout <<ipiv[i]<<" ";
-	cout <<endl;
+	std::cout <<"ipiv = ";
+	for (int i=0; i<n; ++i) std::cout <<ipiv[i]<<" ";
+	std::cout <<std::endl;
 */
 	delete [] ipiv;
 	if (r>0 && B) laerror("singular matrix in lapack_gesv");
@@ -202,7 +207,7 @@ static void linear_solve_do(NRSMat<double> &a, double *b, const int nrhs, const
 	}
 	if (det && r == 0) {
 		*det = 1.;
-		for (int i=1; i<n; i++) {double t=a(i,i); if(!finite(t) || abs(t) < EPSDET ) {*det=0.; break;} else *det *= t;}
+		for (int i=1; i<n; i++) {double t=a(i,i); if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *= t;}
 		//do not use ipiv, since the permutation matrix occurs twice in the decomposition and signs thus cancel (man dspsv)
 	}
         if (det && r>0)  *det = 0;
@@ -452,6 +457,9 @@ void singular_decomposition(NRMat<double> &a, NRMat<double> *u, NRVec<double> &s
 	if (r > 0) laerror("convergence problem in gesvd() of ingular_decomposition()");
 }
 
+//QR decomposition
+//extern "C" void FORNAME(dgeqrf)(const int *M, const int *N, double *A, const int *LDA, double *TAU, double *WORK, int *LWORK, int *INFO);
+
 
 extern "C" void FORNAME(dgeev)(const char *JOBVL, const char *JOBVR, const int *N,
 		double *A, const int *LDA, double *WR, double *WI, double *VL, const int *LDVL,
@@ -555,7 +563,7 @@ void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
 			&ldvr, vl?vl[0]:(double *)0, &ldvl, work, &lwork, &r);
 	delete[] work;
 
-//cout <<"TEST dgeev\n"<<wr<<wi<<*vr<<*vl<<endl;
+//std::cout <<"TEST dgeev\n"<<wr<<wi<<*vr<<*vl<<std::endl;
 
 	if (r < 0) laerror("illegal argument in ggev/geev in gdiagonalize()");
 	if (r > 0) laerror("convergence problem in ggev/geev in gdiagonalize()");
@@ -589,9 +597,9 @@ void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
 				x=.5*t*s11;
 				y=.5*t*s12;
 				double alp,bet;
-				t=.5*sqrt(t);
-				alp=sqrt(.5*(t+x));
-				bet=sqrt(.5*(t-x));
+				t=.5*std::sqrt(t);
+				alp=std::sqrt(.5*(t+x));
+				bet=std::sqrt(.5*(t-x));
 				if(y<0.) bet= -bet;
 
 				//rotate left ev
@@ -722,6 +730,16 @@ const NRMat< complex<double> > imagmatrix (const NRMat<double> &a)
 	return result;
 }
 
+const NRMat< complex<double> > complexmatrix (const NRMat<double> &re, const NRMat<double> &im)
+{
+	if(re.nrows()!=im.nrows() || re.ncols() != im.ncols()) laerror("incompatible sizes of real and imaginary parts");
+        NRMat< complex<double> > result(re.nrows(), re.ncols());
+        cblas_dcopy(re.nrows()*re.ncols(), re, 1, (double *)result[0], 2);
+        cblas_dcopy(re.nrows()*re.ncols(), im, 1, (double *)result[0]+1, 2);
+        return result;
+}
+
+
 
 NRMat<double> matrixfunction(NRMat<double> a, complex<double>
 		(*f)(const complex<double> &), const bool adjust)
@@ -735,13 +753,13 @@ NRMat<complex<double> > a0=complexify(a);
 	gdiagonalize(a, w, &u, &v);//a gets destroyed, eigenvectors are rows
 	NRVec< complex<double> > z = diagofproduct(u, v, 1, 1);
 /*
-cout <<"TEST matrixfunction\n"<<w<<u<<v<<z;
-cout <<"TEST matrixfunction1 "<< u*a0 - diagonalmatrix(w)*u<<endl;
-cout <<"TEST matrixfunction2 "<< a0*v.transpose(1) - v.transpose(1)*diagonalmatrix(w)<<endl;
-cout <<"TEST matrixfunction3 "<< u*v.transpose(1)<<diagonalmatrix(z)<<endl;
+std::cout <<"TEST matrixfunction\n"<<w<<u<<v<<z;
+std::cout <<"TEST matrixfunction1 "<< u*a0 - diagonalmatrix(w)*u<<std::endl;
+std::cout <<"TEST matrixfunction2 "<< a0*v.transpose(1) - v.transpose(1)*diagonalmatrix(w)<<std::endl;
+std::cout <<"TEST matrixfunction3 "<< u*v.transpose(1)<<diagonalmatrix(z)<<std::endl;
 NRVec< complex<double> > wz(n);
 for (int i=0; i<a.nrows(); i++) wz[i] = w[i]/z[i];
-cout <<"TEST matrixfunction4 "<< a0<< v.transpose(true)*diagonalmatrix(wz)*u<<endl;
+std::cout <<"TEST matrixfunction4 "<< a0<< v.transpose(true)*diagonalmatrix(wz)*u<<std::endl;
 */
 
 	for (int i=0; i<a.nrows(); i++) w[i] = (*f)(w[i])/z[i];
@@ -751,7 +769,7 @@ cout <<"TEST matrixfunction4 "<< a0<< v.transpose(true)*diagonalmatrix(wz)*u<<en
 	r.gemm(0.0, v, 'c', u, 'n', 1.0);
 	double inorm = cblas_dnrm2(n*n, (double *)r[0]+1, 2);
 	if (inorm > 1e-10) {
-		cout << "norm = " << inorm << endl;
+		std::cout << "norm = " << inorm << std::endl;
 		laerror("nonzero norm of imaginary part of real matrixfunction");
 	}
 	return realpart(r);
@@ -817,18 +835,18 @@ complex<double> myclog (const complex<double> &x)
 
 complex<double> mycexp (const complex<double> &x)
 {
-        return exp(x);
+        return std::exp(x);
 }
 
 
 complex<double> sqrtinv (const complex<double> &x)
 {
-        return 1./sqrt(x);
+        return 1./std::sqrt(x);
 }
 
 double sqrtinv (const double x)
 {
-        return 1./sqrt(x);
+        return 1./std::sqrt(x);
 }
 
 
@@ -959,7 +977,7 @@ for(j=i; j<n; ++j)
 	if(i==j)
 		{
 		if(x<=0) laerror("not positive definite metric in gendiagonalize");
-       	 	dl[i] = sqrt(x);
+       	 	dl[i] = std::sqrt(x);
 		}
 	else    
 		b(j,i) = x / dl[i];
@@ -1023,12 +1041,13 @@ if(nchange&1)//still adjust to get determinant=1
         {
         int imin=-1; double min=1e200;
         for(int i=0; i<n;++i)
-                if(abs(v[i][i])<min)
+                if(std::abs(v[i][i])<min)
                         {
                         imin=i;
-                        min=abs(v[i][i]);
+                        min=std::abs(v[i][i]);
                         } 
         cblas_dscal(n,-1.,v[imin],1);
         }
 }
 
+}//namespace
diff --git a/nonclass.h b/nonclass.h
index f7caa32..9f7a3ad 100644
--- a/nonclass.h
+++ b/nonclass.h
@@ -17,12 +17,12 @@
 */
 #ifndef _LA_NONCLASS_H_
 #define _LA_NONCLASS_H_
-
 #include "vec.h"
 #include "smat.h"
 #include "mat.h"
 #include "la_traits.h"
 
+namespace LA {
 
 //MISC
 template <class T>
@@ -120,11 +120,11 @@ extern complex<double> sqrtinv(const complex<double> &);
 extern double sqrtinv(const double);
 
 //functions on matrices
-inline NRMat<double>  sqrt(const NRSMat<double> &a) { return matrixfunction(a,&sqrt); }
+inline NRMat<double>  sqrt(const NRSMat<double> &a) { return matrixfunction(a,&std::sqrt); }
 inline NRMat<double>  sqrtinv(const NRSMat<double> &a) { return matrixfunction(a,&sqrtinv); }
-inline NRMat<double>  realsqrt(const NRMat<double> &a) { return realmatrixfunction(a,&sqrt); }
+inline NRMat<double>  realsqrt(const NRMat<double> &a) { return realmatrixfunction(a,&std::sqrt); }
 inline NRMat<double>  realsqrtinv(const NRMat<double> &a) { return realmatrixfunction(a,&sqrtinv); }
-inline NRMat<double>  log(const NRSMat<double> &a) { return matrixfunction(a,&log); }
+inline NRMat<double>  log(const NRSMat<double> &a) { return matrixfunction(a,&std::log); }
 extern NRMat<double> log(const NRMat<double> &a);
 extern NRMat<double> exp0(const NRMat<double> &a);
 
@@ -133,6 +133,7 @@ extern const NRMat<double> realpart(const NRMat< complex<double> >&);
 extern const NRMat<double> imagpart(const NRMat< complex<double> >&);
 extern const NRMat< complex<double> > realmatrix (const NRMat<double>&);
 extern const NRMat< complex<double> > imagmatrix (const NRMat<double>&);
+extern const NRMat< complex<double> > complexmatrix (const NRMat<double>&, const NRMat<double>&);
 
 //inverse by means of linear solve, preserving rhs intact
 template<typename T>
@@ -186,7 +187,7 @@ if(ignoresign)
 {
 for(int i=0; i<nrows; ++i)
         for(int j=0; j<ncols; ++j)
-                r(i,j) = rows[i]*cols[j]<0?0.:a(abs(rows[i])+indexshift,abs(cols[j])+indexshift);
+                r(i,j) = rows[i]*cols[j]<0?0.:a(std::abs(rows[i])+indexshift,std::abs(cols[j])+indexshift);
 }
 else
 {
@@ -201,7 +202,7 @@ if(ignoresign)
 {
 for(int i=0; i<nrows; ++i)
         for(int j=0; j<ncols; ++j)
-                r(i,j) = a(abs(rows[i])+indexshift,abs(cols[j])+indexshift);
+                r(i,j) = a(std::abs(rows[i])+indexshift,std::abs(cols[j])+indexshift);
 }
 else
 {
@@ -218,11 +219,6 @@ return r;
 //auxiliary routine to adjust eigenvectors to guarantee real logarithm
 extern void adjustphases(NRMat<double> &v);
 
-template<class T>
-T abssqr(const complex<T> &x)
-{
-return x.real()*x.real()+x.imag()*x.imag();
-}
 
 //declaration of template interface to cblas routines with full options available
 //just to facilitate easy change between float, double, complex in a user code
@@ -272,6 +268,5 @@ return r;
 }
 
 
-
-
+}//namespace
 #endif
diff --git a/permutation.h b/permutation.h
index 9f86375..189ad94 100644
--- a/permutation.h
+++ b/permutation.h
@@ -15,6 +15,9 @@
     You should have received a copy of the GNU General Public License
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
+#ifndef _PERMUTATION_H
+#define _PERMUTATION_H
+namespace LA {
 template<typename T>
 const NRVec<T> inversepermutation(const NRVec<T> &p, const T offset=0)
 {
@@ -24,3 +27,5 @@ if(!offset) for(int i=0; i<n; ++i) q[p[i]]=i;
 else for(int i=0; i<n; ++i) q[p[i]-offset]=i+offset;
 return q;
 }
+}//namespace
+#endif
diff --git a/qsort.h b/qsort.h
index 55063d6..c7f4cc1 100644
--- a/qsort.h
+++ b/qsort.h
@@ -18,6 +18,8 @@
 #ifndef _QSORT_H
 #define _QSORT_H
 //quicksort, returns parity of the permutation
+//
+namespace LA {
 
 template<typename INDEX, typename COMPAR>
 int genqsort(INDEX l, INDEX r,COMPAR (*cmp)(const INDEX, const INDEX), void (*swap)(const INDEX,const INDEX)) 
@@ -181,4 +183,5 @@ else
 return parity;
 }
 
+}//namespace
 #endif
diff --git a/smat.cc b/smat.cc
index bdc6788..75ddccd 100644
--- a/smat.cc
+++ b/smat.cc
@@ -31,7 +31,7 @@ extern ssize_t write(int, const void *, size_t);
 // TODO
 // specialize unary minus
 
-
+namespace LA {
 
 
 /*
@@ -298,7 +298,7 @@ NRSMat< complex<double> >::dot(const NRSMat< complex<double> > &rhs) const
 	if (nn != rhs.nn) laerror("dot of incompatible SMat's");
 #endif
 	complex<double> dot;
-	cblas_zdotc_sub(NN2, (void *)v, 1, (void *)rhs.v, 1, (void *)(&dot));
+	cblas_zdotc_sub(NN2, v, 1, rhs.v, 1, &dot);
 	return dot;
 }
 
@@ -322,7 +322,7 @@ NRSMat< complex<double> >::dot(const NRVec< complex<double> > &rhs) const
 	if (NN2 != rhs.nn) laerror("dot of incompatible SMat's");
 #endif
 	complex<double> dot;
-	cblas_zdotc_sub(NN2, (void *)v, 1, (void *)rhs.v, 1, (void *)(&dot));
+	cblas_zdotc_sub(NN2, v, 1, rhs.v, 1, &dot);
 	return dot;
 }
 
@@ -341,7 +341,7 @@ const double  NRSMat<double>::norm(const double scalar) const
 			if (i == j) tmp -= scalar;
 			sum += tmp*tmp;
 		}
-	return sqrt(sum);
+	return std::sqrt(sum);
 }
 
 
@@ -350,7 +350,7 @@ template<>
 const double NRSMat< complex<double> >::norm(const complex<double> scalar) const
 {
 	if (!(scalar.real()) && !(scalar.imag()))
-		return cblas_dznrm2(NN2, (void *)v, 1);
+		return cblas_dznrm2(NN2, v, 1);
 	double sum = 0;
 	complex<double> tmp;
 	int k = 0;
@@ -360,7 +360,7 @@ const double NRSMat< complex<double> >::norm(const complex<double> scalar) const
 			if (i == j) tmp -= scalar;
 			sum += tmp.real()*tmp.real() + tmp.imag()*tmp.imag();
 		}
-	return sqrt(sum);
+	return std::sqrt(sum);
 }
 
 
@@ -388,7 +388,7 @@ void NRSMat< complex<double> >::axpy(const complex<double> alpha,
 	if (nn != x.nn) laerror("axpy of incompatible SMats");
 #endif
 	copyonwrite();
-	cblas_zaxpy(nn, (void *)(&alpha), (void *)x.v, 1, (void *)v, 1);
+	cblas_zaxpy(nn, &alpha, x.v, 1, v, 1);
 }
 
 //complex from real
@@ -407,9 +407,16 @@ cblas_dcopy(nn*(nn+1)/2,&rhs(0,0),1,((double *)v) + (imagpart?1:0),2);
 ////// forced instantization in the corresponding object file
 template class NRSMat<double>;
 template class NRSMat< complex<double> >;
+
+template class NRSMat<long long>;
+template class NRSMat<long>;
 template class NRSMat<int>;
 template class NRSMat<short>;
 template class NRSMat<char>;
+template class NRSMat<unsigned char>;
+template class NRSMat<unsigned short>;
 template class NRSMat<unsigned int>;
 template class NRSMat<unsigned long>;
+template class NRSMat<unsigned long long>;
 
+}//namespace
diff --git a/smat.h b/smat.h
index 4520303..1b9b93e 100644
--- a/smat.h
+++ b/smat.h
@@ -19,9 +19,9 @@
 */
 #ifndef _LA_SMAT_H_
 #define _LA_SMAT_H_
-
 #include "la_traits.h"
 
+namespace LA {
 #define NN2 (nn*(nn+1)/2)
 template <class T>
 class NRSMat { // symmetric or complex hermitean matrix in packed form
@@ -94,15 +94,17 @@ public:
 	void fscanf(FILE *f, const char *format); 
 //members concerning sparse matrix
 	explicit NRSMat(const SparseMat<T> &rhs);               // dense from sparse
+	explicit NRSMat(const SparseSMat<T> &rhs);               // dense from sparse
 	inline void simplify() {}; //just for compatibility with sparse ones
 	bool issymmetric() const {return 1;}
 };
 
+}//namespace
 //due to mutual includes this has to be after full class declaration
 #include "vec.h"
 #include "mat.h"
-#include "sparsemat.h"
 
+namespace LA {
 
 // ctors
 template <typename T>
@@ -162,7 +164,7 @@ inline NRSMat< complex<double> > &
 NRSMat< complex<double> >::operator*=(const complex<double> & a)
 {
 	copyonwrite();
-	cblas_zscal(NN2, (void *)(&a), (void *)v, 1);
+	cblas_zscal(NN2, &a, v, 1);
 	return *this;
 }
 template <typename T>
@@ -214,7 +216,7 @@ NRSMat< complex<double> >::operator+=(const NRSMat< complex<double> > & rhs)
 	if (nn != rhs.nn) laerror("incompatible SMats in SMat::operator+=");
 #endif
 	copyonwrite();
-	cblas_zaxpy(NN2, (void *)(&CONE), (void *)(&rhs.v), 1, (void *)(&v), 1);
+	cblas_zaxpy(NN2, &CONE, rhs.v, 1, v, 1);
 	return *this;
 }
 
@@ -251,7 +253,7 @@ NRSMat< complex<double> >::operator-=(const NRSMat< complex<double> > & rhs)
 	if (nn != rhs.nn) laerror("incompatible SMats in SMat::operator-=");
 #endif
 	copyonwrite();
-	cblas_zaxpy(NN2, (void *)(&CMONE), (void *)(&rhs.v), 1, (void *)(&v), 1);
+	cblas_zaxpy(NN2, &CMONE, rhs.v, 1, v, 1);
 	return *this;
 }
 
@@ -406,7 +408,7 @@ inline const double NRSMat<double>::amax() const
 template<>
 inline const complex<double> NRSMat< complex<double> >::amax() const
 {
-	return v[cblas_izamax(NN2, (void *)v, 1)];
+	return v[cblas_izamax(NN2, v, 1)];
 }
 
 // reference pointer to Smat
@@ -554,7 +556,7 @@ return r;
 
 // I/O
 template <typename T>
-ostream& operator<<(ostream &s, const NRSMat<T> &x)
+std::ostream& operator<<(std::ostream &s, const NRSMat<T> &x)
                 {
                 int i,j,n;
                 n=x.nrows();
@@ -568,7 +570,7 @@ ostream& operator<<(ostream &s, const NRSMat<T> &x)
 
 
 template <typename T>
-istream& operator>>(istream  &s, NRSMat<T> &x)
+std::istream& operator>>(std::istream  &s, NRSMat<T> &x)
                 {
                 int i,j,n,m;
                 s >> n >> m;
@@ -618,5 +620,5 @@ public:
 		}
 };
 
-
+}//namespace
 #endif /* _LA_SMAT_H_ */
diff --git a/sparsemat.cc b/sparsemat.cc
index 57c1ae2..c103f13 100644
--- a/sparsemat.cc
+++ b/sparsemat.cc
@@ -25,6 +25,8 @@
 #include <errno.h>
 #include "sparsemat.h"
 
+namespace LA {
+
 template<typename T>
 static inline const T MAX(const T &a, const T &b)
         {return b > a ? (b) : (a);}
@@ -333,11 +335,7 @@ for(i=1; i<n;i++)
 
 //check if  summed to zero
 for(i=0; i<n;i++) if(rowsorted[i] &&
-#ifdef SPARSEEPSILON
-	abs(rowsorted[i]->elem)<SPARSEEPSILON
-#else
-	! rowsorted[i]->elem 
-#endif
+	std::abs(rowsorted[i]->elem)<=SPARSEEPSILON
 	)	 {delete rowsorted[i]; rowsorted[i]=NULL;}
 
 //restore connectivity
@@ -397,11 +395,7 @@ void SparseMat<T>::addsafe(const SPMatindex n, const SPMatindex m, const T elem)
 #ifdef debug
 if(n<0||n>=nn||m<0||m>=mm) laerror("SparseMat out of range");
 #endif
-#ifdef SPARSEEPSILON
-if(abs(elem)<SPARSEEPSILON) return;
-#else
-if(!elem) return;
-#endif
+if(std::abs(elem)<=SPARSEEPSILON) return;
 if(!count) {count=new int;  *count=1; list=NULL;} //blank new matrix
 else copyonwrite(); //makes also unsort
 add(n,m,elem);
@@ -452,9 +446,7 @@ register matel<T> *l=rhs.list;
 while(l)
 	{
 	if(rhs.symmetric || lower && l->row <=l->col || !lower && l->row >=l->col)
-#ifdef SPARSEEPSILON
-	    if(abs(l->elem)>SPARSEEPSILON) 
-#endif
+	    if(std::abs(l->elem)>SPARSEEPSILON) 
 		add( l->row,l->col,sign=='+'?l->elem:- l->elem);
 	l=l->next;
 	}
@@ -473,21 +465,18 @@ register matel<T> *l=rhs.list;
 if(symmetrize)
 while(l)
         {
-#ifdef SPARSEEPSILON
-        if(abs(l->elem)>SPARSEEPSILON) 
-#endif
+        if(std::abs(l->elem)>SPARSEEPSILON) 
 		{add( l->row,l->col,l->elem); if( l->row!=l->col) add( l->col,l->row,l->elem);}
         l=l->next;
         }
 else
 while(l)
 	{
-#ifdef SPARSEEPSILON
-	if(abs(l->elem)>SPARSEEPSILON) 
-#endif
+	if(std::abs(l->elem)>SPARSEEPSILON) 
 		add( l->row,l->col,l->elem);
 	l=l->next;
 	}
+simplify(); //maybe leave up to the user
 return *this;
 }
 
@@ -503,21 +492,18 @@ register matel<T> *l=rhs.list;
 if(symmetrize)
 while(l)
         {
-#ifdef SPARSEEPSILON
-        if(abs(l->elem)>SPARSEEPSILON) 
-#endif
+        if(std::abs(l->elem)>SPARSEEPSILON) 
 		{add( l->row,l->col,- l->elem); if( l->row!=l->col) add( l->col,l->row,- l->elem);}
         l=l->next;
         }
 else
 while(l)
         {
-#ifdef SPARSEEPSILON
-        if(abs(l->elem)>SPARSEEPSILON) 
-#endif
+        if(std::abs(l->elem)>SPARSEEPSILON) 
 		add( l->row,l->col,- l->elem);
         l=l->next;
         }
+simplify(); //maybe leave up to the user
 return *this;
 }
 
@@ -537,11 +523,7 @@ SPMatindex i,j;
 for(i=0;i<nn;++i)
 	for(j=0; j<mm;++j)
 		{register T t(rhs(i,j));
-#ifdef SPARSEEPSILON
-		if( abs(t)>SPARSEEPSILON)
-#else
-		if(t)
-#endif
+		if( std::abs(t)>SPARSEEPSILON)
 		 add(i,j,t);
 		}
 }
@@ -660,40 +642,38 @@ else while(l)
 
 //assignment of a scalar matrix
 template <class T>
-SparseMat<T> & SparseMat<T>::operator=(const T a)
+SparseMat<T> & SparseMat<T>::operator=(const T &a)
 {
 if(!count ||nn<=0||mm<=0) laerror("assignment of scalar to undefined sparse matrix");
 if(nn!=mm) laerror("assignment of scalar to non-square sparse matrix");
 resize(nn,mm);//clear 
-#ifdef SPARSEEPSILON
-if(abs(a)<SPARSEEPSILON) return *this;
-#else
-if(a==(T)0) return *this;
-#endif
+if(std::abs(a)<=SPARSEEPSILON) return *this;
 SPMatindex i;
 for(i=0;i<nn;++i) add(i,i,a);
 return *this;
 }
 
 template <class T>
-SparseMat<T> & SparseMat<T>::operator+=(const T a)
+SparseMat<T> & SparseMat<T>::operator+=(const T &a)
 {
 if(!count ||nn<=0||mm<=0) laerror("assignment of scalar to undefined sparse matrix");
 if(nn!=mm) laerror("assignment of scalar to non-square sparse matrix");
 if(a==(T)0) return *this;
 SPMatindex i;
 for(i=0;i<nn;++i) add(i,i,a);
+simplify(); //maybe leave up to the user
 return *this;
 }
 
 template <class T>
-SparseMat<T> & SparseMat<T>::operator-=(const T a)
+SparseMat<T> & SparseMat<T>::operator-=(const T &a)
 {
 if(!count ||nn<=0||mm<=0) laerror("assignment of scalar to undefined sparse matrix");
 if(nn!=mm) laerror("assignment of scalar to non-square sparse matrix");
 if(a==(T)0) return *this;
 SPMatindex i;
 for(i=0;i<nn;++i) add(i,i,-a);
+simplify(); //maybe leave up to the user
 return *this;
 }
 
@@ -715,11 +695,7 @@ for(i=0;i<nn;++i)
         for(j=0; j<=i;++j)
 		{register T t;
                 if(
-#ifdef SPARSEEPSILON
-			abs(t=rhs(i,j))>SPARSEEPSILON
-#else
-			t=rhs(i,j)
-#endif
+			std::abs(t=rhs(i,j))>SPARSEEPSILON
 			) add(i,j,t);
 		}
 }
@@ -751,9 +727,7 @@ matel<T> *l=list;
 while(l) //include the mirror picture of elements into the list
 	{
 	if(
-#ifdef SPARSEEPSILON
-		abs(l->elem)>SPARSEEPSILON && 
-#endif
+		std::abs(l->elem)>SPARSEEPSILON && 
 	   l->row!=l->col) add(l->col,l->row,l->elem); //not OK for complex-hermitean
 	l=l->next;
 	}
@@ -761,12 +735,12 @@ while(l) //include the mirror picture of elements into the list
 
 
 template <class T>
-SparseMat<T> & SparseMat<T>::operator*=(const T a)
+SparseMat<T> & SparseMat<T>::operator*=(const T &a)
 {
 if(!count) laerror("operator*= on undefined lhs");
 if(!list||a==(T)1) return *this;
-if(a==(T)0) resize(nn,mm);
-else copyonwrite();
+if(a==(T)0) {clear(); return *this;}
+copyonwrite();
 
 register matel<T> *l=list;
 while(l) 
@@ -936,15 +910,14 @@ return sum;
 }
 
 
-//not OK for complex hermitean
 template<class T>
-const T SparseMat<T>::norm(const T scalar) const
+const typename LA_traits<T>::normtype SparseMat<T>::norm(const T scalar) const
 {
-if(!list) return T(0);
+if(!list) return typename LA_traits<T>::normtype(0);
 const_cast<SparseMat<T> *>(this)->simplify(); 
 
 matel<T> *l=list;
-T sum(0);
+typename LA_traits<T>::normtype sum(0);
 if(scalar!=(T)0)
 	{
 	if(symmetric) 
@@ -953,7 +926,7 @@ if(scalar!=(T)0)
 		T hlp=l->elem;
 		bool b=l->row==l->col;
 		if(b) hlp-=scalar;
-		T tmp=hlp*hlp;
+		typename LA_traits<T>::normtype tmp=LA_traits<T>::sqrabs(hlp);
        	 sum+= tmp;
 		if(!b) sum+=tmp;
        	 l=l->next;
@@ -963,7 +936,7 @@ if(scalar!=(T)0)
 	       	 {
 		T hlp=l->elem;
 		if(l->row==l->col) hlp-=scalar;
-	        sum+= hlp*hlp; 
+	        sum+= LA_traits<T>::sqrabs(hlp);
 	        l=l->next;
 	        }
 	}
@@ -972,7 +945,7 @@ else
         if(symmetric)
         while(l)
                 {
-                T tmp=l->elem*l->elem;
+                typename LA_traits<T>::normtype tmp=LA_traits<T>::sqrabs(l->elem);
          sum+= tmp;
                 if(l->row!=l->col) sum+=tmp;
          l=l->next;
@@ -980,11 +953,11 @@ else
         else
         while(l)
                  {
-                sum+= l->elem*l->elem;
+                sum+= LA_traits<T>::sqrabs(l->elem);
                 l=l->next;
                 }
 	}
-return sqrt(sum); //not OK for int, would need traits technique
+return (typename LA_traits<T>::normtype) std::sqrt(sum); //not OK for int, would need traits technique
 }
 
 
@@ -1009,9 +982,7 @@ if(transp)
 while(l)
         {
 	register T t=alpha*l->elem;
-#ifdef SPARSEEPSILON
-	if(abs(t)>SPARSEEPSILON) 
-#endif
+	if(std::abs(t)>SPARSEEPSILON) 
 				add( l->col,l->row,t);
         l=l->next;
         }
@@ -1019,9 +990,7 @@ else
 while(l)
 	{
 	register T t=alpha*l->elem;
-#ifdef SPARSEEPSILON
-	if(abs(t)>SPARSEEPSILON) 
-#endif
+	if(std::abs(t)>SPARSEEPSILON) 
 				add( l->row,l->col,t);
 	l=l->next;
 	}
@@ -1227,11 +1196,7 @@ incsize(rhs.nrows(), rhs.ncols());
 T tmp;
 for(SPMatindex i=0; i<(SPMatindex)rhs.nrows(); ++i)
 	for(SPMatindex j=0; j<(SPMatindex)rhs.ncols(); ++j)
-#ifdef SPARSEEPSILON
-		if(abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
-#else
-		if((tmp=rhs(i,j))!=(T)0) add(n0+i,m0+j,tmp);
-#endif
+		if(std::abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
 return *this;
 }
 
@@ -1247,11 +1212,7 @@ T tmp;
 incsize(rhs.nrows(), rhs.ncols());
 for(SPMatindex i=0; i<(SPMatindex)rhs.nrows(); ++i)
         for(SPMatindex j=0; j<(SPMatindex)rhs.ncols(); ++j)
-#ifdef SPARSEEPSILON
-                if(abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
-#else
-                if((tmp=rhs(i,j))!=(T)0) add(n0+i,m0+j,tmp);
-#endif
+                if(std::abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
 return *this;
 }
 
@@ -1270,9 +1231,7 @@ incsize(rhs.nrows(), rhs.ncols());
 register matel<T> *l=rhs.list;
 while(l)
 	{
-#ifdef SPARSEEPSILON
-	if(abs(l->elem)>SPARSEEPSILON) 
-#endif
+	if(std::abs(l->elem)>SPARSEEPSILON) 
 		add(n0+l->row,m0+l->col,l->elem);
 	l=l->next;
 	}
@@ -1303,18 +1262,18 @@ template SparseMat<T>::SparseMat(const NRMat<T> &rhs); \
 template SparseMat<T>::SparseMat(const NRSMat<T> &rhs); \
 template void SparseMat<T>::transposeme(); \
 template const T* SparseMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const; \
-template SparseMat<T> & SparseMat<T>::operator*=(const T a); \
+template SparseMat<T> & SparseMat<T>::operator*=(const T &a); \
 template void SparseMat<T>::setunsymmetric(); \
-template SparseMat<T> & SparseMat<T>::operator=(const T a); \
-template SparseMat<T> & SparseMat<T>::operator+=(const T a); \
-template SparseMat<T> & SparseMat<T>::operator-=(const T a); \
+template SparseMat<T> & SparseMat<T>::operator=(const T &a); \
+template SparseMat<T> & SparseMat<T>::operator+=(const T &a); \
+template SparseMat<T> & SparseMat<T>::operator-=(const T &a); \
 template NRMat<T>::NRMat(const SparseMat<T> &rhs); \
 template NRSMat<T>::NRSMat(const SparseMat<T> &rhs); \
 template NRVec<T>::NRVec(const SparseMat<T> &rhs); \
 template const NRVec<T> NRVec<T>::operator*(const SparseMat<T> &mat) const; \
 template SparseMat<T> & SparseMat<T>::join(SparseMat<T> &rhs); \
 template const T SparseMat<T>::trace() const; \
-template const T SparseMat<T>::norm(const T scalar) const; \
+template const LA_traits<T>::normtype SparseMat<T>::norm(const T scalar) const; \
 template void SparseMat<T>::axpy(const T alpha, const SparseMat<T> &x, const bool transp); \
 template const SparseMat<T> SparseMat<T>::operator*(const SparseMat<T> &rhs) const; \
 template const T SparseMat<T>::dot(const SparseMat<T> &rhs) const; \
@@ -1336,3 +1295,4 @@ INSTANTIZE(complex<double>) //some functions are not OK for hermitean matrices,
 template class SparseMat<double>;
 template class SparseMat<complex<double> >;
 
+}//namespace
diff --git a/sparsemat.h b/sparsemat.h
index 6633827..ece4a04 100644
--- a/sparsemat.h
+++ b/sparsemat.h
@@ -17,13 +17,13 @@
 */
 #ifndef _SPARSEMAT_H_
 #define _SPARSEMAT_H_
-
 #include "la_traits.h"
 
+namespace LA {
 
 //threshold for neglecting elements, if not defined, no tests are done except exact zero test in simplify - might be even faster
 //seems to perform better with a threshold, in spite of abs() tests
-#define  SPARSEEPSILON 1e-13 
+const double SPARSEEPSILON=1e-35;
 
 typedef unsigned int SPMatindex;
 typedef int SPMatindexdiff; //more clear would be to use traits
@@ -81,10 +81,10 @@ public:
 	explicit SparseMat(const NRMat<T> &rhs); //construct from a dense one
 	explicit SparseMat(const NRSMat<T> &rhs); //construct from a dense symmetric one
 	SparseMat & operator=(const SparseMat &rhs);
-	SparseMat & operator=(const T a);          //assign a to diagonal
-    	SparseMat & operator+=(const T a);         //assign a to diagonal
-	SparseMat & operator-=(const T a);         //assign a to diagonal
-        SparseMat & operator*=(const T a);         //multiply by a scalar
+	SparseMat & operator=(const T &a);          //assign a to diagonal
+    	SparseMat & operator+=(const T &a);         //assign a to diagonal
+	SparseMat & operator-=(const T &a);         //assign a to diagonal
+        SparseMat & operator*=(const T &a);         //multiply by a scalar
         SparseMat & operator+=(const SparseMat &rhs);
 	SparseMat & addtriangle(const SparseMat &rhs, const bool lower, const char sign);
         SparseMat & join(SparseMat &rhs); //more efficient +=, rhs will be emptied
@@ -135,17 +135,21 @@ public:
 	void clear() {copyonwrite();unsort();deletelist();}
 	void simplify();
 	const T trace() const;
-	const T norm(const T scalar=(T)0) const; //is const only mathematically, not in internal implementation - we have to simplify first
+	const typename LA_traits<T>::normtype norm(const T scalar=(T)0) const; //is const only mathematically, not in internal implementation - we have to simplify first
 	unsigned int sort(int type) const;
 	inline void add(const SPMatindex n, const SPMatindex m, const T elem) {matel<T> *ltmp= new matel<T>; ltmp->next=list; list=ltmp; list->row=n; list->col=m; list->elem=elem;}
 	void addsafe(const SPMatindex n, const SPMatindex m, const T elem);
 };
 
+}//namespace
+
 //due to mutual includes this has to be after full class declaration
 #include "vec.h"
 #include "smat.h"
 #include "mat.h"
 
+namespace LA {
+
 template <typename T>
 inline const NRVec<T> SparseMat<T>::operator*(const NRVec<T> &rhs) const
 {NRVec<T> result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;};
@@ -155,7 +159,7 @@ 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 <class T>
-ostream& operator<<(ostream &s, const SparseMat<T> &x)
+std::ostream& operator<<(std::ostream &s, const SparseMat<T> &x)
                 {
                 SPMatindex n,m;
                 n=x.nrows();
@@ -172,7 +176,7 @@ ostream& operator<<(ostream &s, const SparseMat<T> &x)
                 }
 
 template <class T>
-istream& operator>>(istream  &s, SparseMat<T> &x)
+std::istream& operator>>(std::istream  &s, SparseMat<T> &x)
                 {
                 int i,j;
 		int n,m;
@@ -298,4 +302,5 @@ while(l)
 return *this;
 }
 
+}//namespace
 #endif
diff --git a/sparsesmat.cc b/sparsesmat.cc
index 248bc9e..0ecaa39 100644
--- a/sparsesmat.cc
+++ b/sparsesmat.cc
@@ -25,10 +25,11 @@
 #include <errno.h>
 #include "sparsesmat.h"
 
+namespace LA {
+
 template <typename T>
 void SparseSMat<T>::gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha)
 {
-std::cerr << "enter gemm\n";
 (*this) *= beta;
 if(alpha==(T)0) return;
 if(a.nn!=b.nn || a.nn!=nn) laerror("incompatible sizes in SparseSMat::gemm");
@@ -64,11 +65,7 @@ for(SPMatindex k=0; k<nn; ++k) //summation loop
 		add(ai[i],bi[j],prod(i,j),false);
 
 	}
-std::cerr << "before simplify in gemm\n";
-std::cerr << " nterms = "<<
-simplify() 
-<<std::endl; //erase elements below threshold
-std::cerr << "regular exit gemm\n";
+simplify();
 }
 
 
@@ -183,7 +180,7 @@ for(SPMatindex i=0; i<nn; ++i)
 		}
 	else sum += LA_traits<T>::sqrabs(scalar); //missing diagonal element
 
-return sqrt(sum);
+return std::sqrt(sum);
 }
 
 
@@ -208,4 +205,4 @@ INSTANTIZE(complex<double>)
 template class SparseSMat<double>;
 template class SparseSMat<complex<double> >;
 
-
+}//namespace
diff --git a/sparsesmat.h b/sparsesmat.h
index 298b8b3..58bd828 100644
--- a/sparsesmat.h
+++ b/sparsesmat.h
@@ -19,7 +19,6 @@
 #ifndef _SPARSESMAT_H_
 #define _SPARSESMAT_H_
 
-
 #include <string>
 #include <cmath>
 #include <stdlib.h>
@@ -36,6 +35,7 @@
 #include <map>
 #include <list>
 
+namespace LA {
 
 //symmetric sparse matrix class with a representation for efficient exponentiatiation
 //in particular we need a unitary symmetric complex matrix as exp(iH) with H real symmetric
@@ -360,5 +360,5 @@ std::istream& operator>>(std::istream  &s, SparseSMat<T> &x)
                 }
 
 
-
+}//namespace
 #endif //_SPARSESMAT_H_
diff --git a/strassen.cc b/strassen.cc
index fa55468..9a350ea 100644
--- a/strassen.cc
+++ b/strassen.cc
@@ -20,6 +20,7 @@
 #include "mat.h"
 
 #ifndef NO_STRASSEN
+namespace LA {
 /*Strassen algorithm*/
 // called routine is fortran-compatible
 extern "C" void fmm(const char c_transa,const char c_transb,const int m,const int n,const int k,const double alpha,
@@ -45,4 +46,5 @@ copyonwrite();
 //swap transpositions and order of matrices
 fmm(transb,transa,mm,nn,k,alpha,b,b.mm, a, a.mm, beta,*this, mm,NULL,0);
 }
+}//namespace
 #endif
diff --git a/t.cc b/t.cc
index 1cb8363..183f900 100644
--- a/t.cc
+++ b/t.cc
@@ -20,18 +20,9 @@
 
 #include <time.h>
 #include "la.h"
-#include "sparsemat.h"
-#include "matexp.h"
-#include "fourindex.h"
-#include "davidson.h"
-#include "gmres.h"
-#include "conjgrad.h"
-#include "diis.h"
-#include "bitvector.h"
-#ifdef USE_TRACEBACK
-#include "traceback.h"
-#endif
 
+using namespace std;
+using namespace LA;
 
 extern void test(const NRVec<double> &x);
 
@@ -1335,7 +1326,7 @@ cout << "Unitary\n"<<U;
 cout <<"Error = "<<(U*U.transpose() -1.).norm()<<endl;
 }
 
-if(1)
+if(0)
 {
 double t;
 NRMat<complex<double> > ham;
@@ -1397,7 +1388,58 @@ cout <<"sparse sincostimes "<<clock()/((double) (CLOCKS_PER_SEC))-t <<"\n";
 NRVec<complex<double> > r4 = cov2 + siv2*complex<double>(0,1);
 cout <<"sincostimes error 2 = "<<(r1-r4).norm()<<endl;
 */
+}
 
+if(0)
+{
+double t;
+SparseMat<complex<double> > h;
+cin >> h;
+h *= complex<double>(0,1);
+h.simplify();
+cout <<"length of hamiltonian "<<h.length()<<endl;
+t=clock()/((double) (CLOCKS_PER_SEC));
+SparseMat<complex<double> > hexp = exp(h);
+cout <<"sparse exp time "<<clock()/((double) (CLOCKS_PER_SEC))-t <<"\n";
+cout <<"length of exp hamiltonian "<<hexp.length()<<endl;
+NRMat<complex<double> > he(hexp);
+NRMat<complex<double> > hec(he);
+NRMat<complex<double> > het(he);
+hec.conjugateme();
+het.transposeme();
+//cout <<he;
+cout <<"unitarity error "<<(hec*he).norm(1)<<endl;
+cout <<"symmetry error "<<(het-he).norm()<<endl;
+}
+
+if(0)
+{
+NRSMat<double> hd(100);
+hd.randomize(1);
+SparseSMat<double> h(hd);
+NRMat<double> rd = hd*hd;
+SparseSMat<double> r = h*h;
+NRSMat<double>rx(r);
+NRMat<double> r2(rx);
+cout <<"Error = "<<(r2-rd).norm()<<endl;
+}
+
+if(1)
+{
+SparseSMat<complex<double> > h;
+cin>>h;
+h *= complex<double>(0,1);
+double t=clock()/((double) (CLOCKS_PER_SEC));
+SparseSMat<complex<double> > r = exp(h);
+cout <<"SparseSMat time "<<clock()/((double) (CLOCKS_PER_SEC))-t <<"\n";
+if(h.nrows()<=1024)
+{
+NRSMat<complex<double> > h3(h);
+NRMat<complex<double> > h2(h3);
+NRMat<complex<double> >r2 = exp(h2);
+cout <<"error = "<<(r2-NRMat<complex<double> >(NRSMat<complex<double> >(r))).norm()<<endl;
+cout <<"errorx = "<<(r2-NRSMat<complex<double> >(r)).norm()<<endl;
+}
 }
 
 
diff --git a/t2.cc b/t2.cc
index 7253842..e845056 100644
--- a/t2.cc
+++ b/t2.cc
@@ -28,13 +28,15 @@
 #include "traceback.h"
 #endif
 
+using namespace LA;
+
 void test(const NRVec<double> &x)
 {
 NRMat<double> aa(0.,2,2);
 NRMat<double> cc(aa);
 cc.copyonwrite(); cc[0][0]=1.;
 
-cout << aa << cc <<"\n";
-cout <<"test x= "<<x<<"\n";
+std::cout << aa << cc <<"\n";
+std::cout <<"test x= "<<x<<"\n";
 }
 
diff --git a/test.cc b/test.cc
index 6649a46..f660e8e 100644
--- a/test.cc
+++ b/test.cc
@@ -25,6 +25,8 @@
 #include "traceback.h"
 #endif
 
+using namespace std;
+
 void test2(char *msg)
 {
 laerror(msg);
diff --git a/vec.cc b/vec.cc
index 41708ff..03adfe0 100644
--- a/vec.cc
+++ b/vec.cc
@@ -18,19 +18,22 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+
 #include <iostream>
-#include "vec.h"
 #include <stdlib.h>
 #include <sys/types.h>
 #include <sys/stat.h>
 #include <fcntl.h>
 #include <errno.h>
+#include "vec.h"
 #include "qsort.h"
 extern "C" {
 extern ssize_t read(int, void *, size_t);
 extern ssize_t write(int, const void *, size_t);
 }
 
+namespace LA {
+
 
 //////////////////////////////////////////////////////////////////////////////
 //// forced instantization in the corespoding object file
@@ -42,13 +45,16 @@ template void NRVec<T>::get(int fd, bool dim, bool transp); \
 
 INSTANTIZE(double)
 INSTANTIZE(complex<double>)
-INSTANTIZE(int)
-INSTANTIZE(short)
-INSTANTIZE(unsigned short)
 INSTANTIZE(char)
+INSTANTIZE(short)
+INSTANTIZE(int)
+INSTANTIZE(long)
+INSTANTIZE(long long)
 INSTANTIZE(unsigned char)
+INSTANTIZE(unsigned short)
 INSTANTIZE(unsigned int)
 INSTANTIZE(unsigned long)
+INSTANTIZE(unsigned long long)
 
 
 
@@ -206,7 +212,7 @@ void NRVec< complex<double> >::axpy(const complex<double> alpha,
 	if (nn != x.nn) laerror("axpy of incompatible vectors");
 #endif
 	copyonwrite();
-	cblas_zaxpy(nn, (void *)(&alpha), (void *)(x.v), 1, (void *)v, 1);
+	cblas_zaxpy(nn, &alpha, x.v, 1, v, 1);
 }
 
 // axpy call for T = double (strided)
@@ -223,7 +229,7 @@ void NRVec< complex<double> >::axpy(const complex<double> alpha,
 			const complex<double> *x, const int stride)
 {
 	copyonwrite();
-	cblas_zaxpy(nn, (void *)(&alpha), (void *)x, stride, v, 1);
+	cblas_zaxpy(nn, &alpha, x, stride, v, 1);
 }
 
 // unary minus
@@ -242,7 +248,7 @@ NRVec< complex<double> >::operator-() const
 {
 	NRVec< complex<double> > result(*this);
 	result.copyonwrite();
-	cblas_zdscal(nn, -1.0, (void *)(result.v), 1);
+	cblas_zdscal(nn, -1.0, result.v, 1);
 	return result;
 }
 
@@ -279,37 +285,26 @@ template<>
 NRVec< complex<double> > & NRVec< complex<double> >::normalize()
 {
 	complex<double> tmp;
-	tmp = cblas_dznrm2(nn, (void *)v, 1);
+	tmp = cblas_dznrm2(nn, v, 1);
 #ifdef DEBUG
 	if(!(tmp.real()) && !(tmp.imag())) laerror("normalization of zero vector");
 #endif
 	copyonwrite();
 	tmp = 1.0/tmp;
-	cblas_zscal(nn, (void *)(&tmp), (void *)v, 1);
+	cblas_zscal(nn, &tmp, v, 1);
 	return *this;
 }
 
 //stubs for linkage
-template<>
-NRVec<unsigned char> &  NRVec<unsigned char>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
-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<>
-NRVec<unsigned long> &  NRVec<unsigned long>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
-template<>
-NRVec<unsigned int> &  NRVec<unsigned int>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
-
 
 #define INSTANTIZE_DUMMY(T) \
 template<> void NRVec<T>::gemv(const T beta, const NRMat<T> &a, const char trans,  const T alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } \
 template<> void NRVec<T>::gemv(const T beta, const NRSMat<T> &a, const char trans,  const T alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } \
 template<> void NRVec<T>::gemv(const T beta, const SparseMat<T> &a, const char trans,  const T alpha, const NRVec<T> &x, bool s) { laerror("gemv on unsupported types"); } \
 template<> void NRVec<T>::gemv(const LA_traits_complex<T>::Component_type beta, const  LA_traits_complex<T>::NRMat_Noncomplex_type  &a, const char trans,  const  LA_traits_complex<T>::Component_type alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } \
-template<> void NRVec<T>::gemv(const  LA_traits_complex<T>::Component_type beta, const  LA_traits_complex<T>::NRSMat_Noncomplex_type  &a, const char trans,  const  LA_traits_complex<T>::Component_type alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } 
+template<> void NRVec<T>::gemv(const  LA_traits_complex<T>::Component_type beta, const  LA_traits_complex<T>::NRSMat_Noncomplex_type  &a, const char trans,  const  LA_traits_complex<T>::Component_type alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } \
+template<> NRVec<T> & NRVec<T>::normalize() {laerror("normalize() impossible for integer types"); return *this;} \
+template<> const NRMat<T> NRVec<T>::otimes(const NRVec<T> &b,const bool conj, const T &scale) const {laerror("otimes presently implemented only for double and complex double"); return NRMat<T> ();}
 
 
 
@@ -401,20 +396,22 @@ void NRVec< complex<double> >::gemv(const complex<double> beta,
 
 
 
-// Direc product Mat = Vec | Vec
+// Direct product Mat = Vec | Vec
 template<>
-const NRMat<double> NRVec<double>::operator|(const NRVec<double> &b) const
+const NRMat<double> NRVec<double>::otimes(const NRVec<double> &b,const bool conj, const double &scale) const
 {
 	NRMat<double> result(0.,nn,b.nn);
-	cblas_dger(CblasRowMajor, nn, b.nn, 1., v, 1, b.v, 1, result, b.nn);
+	cblas_dger(CblasRowMajor, nn, b.nn, scale, 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
+NRVec<complex<double> >::otimes(const NRVec< complex<double> > &b, const bool conj, const  complex<double> &scale) const
 {
 	NRMat< complex<double> > result(0.,nn,b.nn);
-	cblas_zgerc(CblasRowMajor, nn, b.nn, &CONE, v, 1, b.v, 1, result, b.nn);
+	if(conj) cblas_zgerc(CblasRowMajor, nn, b.nn, &scale, v, 1, b.v, 1, result, b.nn);
+	else cblas_zgeru(CblasRowMajor, nn, b.nn, &scale, v, 1, b.v, 1, result, b.nn);
 	return result;
 }
 
@@ -436,24 +433,39 @@ else return memqsort<0,NRVec<T>,int,int>(*this,perm,from,to);
 template class NRVec<double>;
 template class NRVec<complex<double> >;
 template class NRVec<char>;
-template class NRVec<unsigned char>;
 template class NRVec<short>;
 template class NRVec<int>;
+template class NRVec<long>;
+template class NRVec<long long>;
+template class NRVec<unsigned char>;
+template class NRVec<unsigned short>;
 template class NRVec<unsigned int>;
 template class NRVec<unsigned long>;
+template class NRVec<unsigned long long>;
 
 
 INSTANTIZE_DUMMY(char)
-INSTANTIZE_DUMMY(unsigned char)
 INSTANTIZE_DUMMY(short)
 INSTANTIZE_DUMMY(int)
+INSTANTIZE_DUMMY(long)
+INSTANTIZE_DUMMY(long long)
+INSTANTIZE_DUMMY(unsigned char)
+INSTANTIZE_DUMMY(unsigned short)
 INSTANTIZE_DUMMY(unsigned int)
 INSTANTIZE_DUMMY(unsigned long)
+INSTANTIZE_DUMMY(unsigned long long)
 INSTANTIZE_DUMMY(complex<char>)
-INSTANTIZE_DUMMY(complex<unsigned char>)
 INSTANTIZE_DUMMY(complex<short>)
 INSTANTIZE_DUMMY(complex<int>)
+INSTANTIZE_DUMMY(complex<long>)
+INSTANTIZE_DUMMY(complex<long long>)
+INSTANTIZE_DUMMY(complex<unsigned char>)
+INSTANTIZE_DUMMY(complex<unsigned short>)
 INSTANTIZE_DUMMY(complex<unsigned int>)
+INSTANTIZE_DUMMY(complex<unsigned long>)
+INSTANTIZE_DUMMY(complex<unsigned long long>)
+
 INSTANTIZE_DUMMY(complex<complex<double> >)
 INSTANTIZE_DUMMY(complex<complex<float> >)
-INSTANTIZE_DUMMY(complex<unsigned long>)
+
+}//namespace
diff --git a/vec.h b/vec.h
index 6ce250a..75d5b7b 100644
--- a/vec.h
+++ b/vec.h
@@ -21,6 +21,8 @@
 
 #include "la_traits.h"
 
+namespace LA {
+
 //////////////////////////////////////////////////////////////////////////////
 // Forward declarations
 template <typename T> void lawritemat(FILE *file,const T *a,int r,int c,
@@ -104,7 +106,8 @@ public:
 	const NRVec operator*(const NRMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
 	const NRVec operator*(const NRSMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
 	const NRVec operator*(const SparseMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
-	const NRMat<T> operator|(const NRVec<T> &rhs) const;
+	const NRMat<T> otimes(const NRVec<T> &rhs, const bool conjugate=false, const T &scale=1) const; //outer product
+	inline const NRMat<T> operator|(const NRVec<T> &rhs) const {return otimes(rhs,true);};
 	inline const T sum() const {T sum=0; for(int i=0; i<nn; i++) sum += v[i]; return sum;}; //sum of its elements
 	inline const T asum() const; //sum of its elements absolute values
 	inline const T dot(const T *a, const int stride=1) const; // ddot with a stride-vector
@@ -136,26 +139,28 @@ public:
 	int sort(int direction=0, int from=0, int to= -1, int *perm=NULL); //sort, ascending by default, returns parity of permutation
 };
 
+}//namespace
 
 //due to mutual includes this has to be after full class declaration
 #include "mat.h"
 #include "smat.h"
 #include "sparsemat.h"
 
+namespace LA {
 // formatted I/O
 template <typename T>
-ostream & operator<<(ostream &s, const NRVec<T> &x)
+std::ostream & operator<<(std::ostream &s, const NRVec<T> &x)
 {
   int i, n;
 
   n = x.size();
-  s << n << endl;
+  s << n << std::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)
+std::istream & operator>>(std::istream &s, NRVec<T> &x)
 {
   int i,n;
 
@@ -238,7 +243,7 @@ inline NRVec< complex<double> > &
 NRVec< complex<double> >::operator+=(const complex<double> &a)
 {
 	copyonwrite();
-	cblas_zaxpy(nn, (void *)(&CONE), (void *)(&a), 0, (void *)v, 1);
+	cblas_zaxpy(nn, &CONE, &a, 0, v, 1);
 	return *this;
 }
 
@@ -267,7 +272,7 @@ inline NRVec< complex<double> > &
 NRVec< complex<double> >::operator-=(const complex<double> &a)
 {
 	copyonwrite();
-	cblas_zaxpy(nn, (void *)(&CMONE), (void *)(&a), 0, (void *)v, 1);
+	cblas_zaxpy(nn, &CMONE, &a, 0, v, 1);
 	return *this;
 }
 
@@ -302,7 +307,7 @@ NRVec< complex<double> >::operator+=(const NRVec< complex<double> > &rhs)
 	if (nn != rhs.nn) laerror("daxpy of incompatible vectors");
 #endif
 	copyonwrite();
-	cblas_zaxpy(nn, (void *)(&CONE), rhs.v, 1, v, 1);
+	cblas_zaxpy(nn, &CONE, rhs.v, 1, v, 1);
 	return *this;
 }
 
@@ -366,7 +371,7 @@ NRVec< complex<double> >::operator-=(const NRVec< complex<double> > &rhs)
 	if (nn != rhs.nn) laerror("daxpy of incompatible vectors");
 #endif
 	copyonwrite();
-	cblas_zaxpy(nn, (void *)(&CMONE), (void *)rhs.v, 1, (void *)v, 1);
+	cblas_zaxpy(nn, &CMONE, rhs.v, 1, v, 1);
 	return *this;
 }
 
@@ -398,7 +403,7 @@ inline NRVec< complex<double> > &
 NRVec< complex<double> >::operator*=(const complex<double> &a)
 {
 	copyonwrite();
-	cblas_zscal(nn, (void *)(&a), (void *)v, 1);
+	cblas_zscal(nn, &a, v, 1);
 	return *this;
 }
 
@@ -432,7 +437,7 @@ NRVec< complex<double> >::operator*(const NRVec< complex<double> > &rhs) const
 	if (nn != rhs.nn) laerror("dot of incompatible vectors");
 #endif
 	complex<double> dot;
-	cblas_zdotc_sub(nn, (void *)v, 1, (void *)rhs.v, 1, (void *)(&dot));
+	cblas_zdotc_sub(nn, v, 1, rhs.v, 1, &dot);
 	return dot;
 }
 
@@ -468,7 +473,7 @@ inline const complex<double>
 NRVec< complex<double> >::dot(const complex<double> *y, const int stride) const
 {
 	complex<double> dot;
-	cblas_zdotc_sub(nn, y, stride, v, 1, (void *)(&dot));
+	cblas_zdotc_sub(nn, y, stride, v, 1, &dot);
 	return dot;
 }
 
@@ -527,7 +532,7 @@ inline const double  NRVec<double>::norm() const
 template<>
 inline const double NRVec< complex<double> >::norm() const
 {
-	return cblas_dznrm2(nn, (void *)v, 1);
+	return cblas_dznrm2(nn, v, 1);
 }
 
 // Max element of the array
@@ -539,7 +544,7 @@ inline const double  NRVec<double>::amax() const
 template<>
 inline const complex<double> NRVec< complex<double> >::amax() const
 {
-	return v[cblas_izamax(nn, (void *)v, 1)];
+	return v[cblas_izamax(nn, v, 1)];
 }
 
 
@@ -687,6 +692,6 @@ for(int i=0; i<rhs.size(); ++i)  r[i]=rhs[i];
 return r;
 }
 
-
+}//namespace
 
 #endif /* _LA_VEC_H_ */