diff --git a/fourindex.h b/fourindex.h
index e3dd4ff..ff87971 100644
--- a/fourindex.h
+++ b/fourindex.h
@@ -619,6 +619,11 @@ fourindex_dense<twoelectronrealmullikan,T,I>::fourindex_dense<twoelectronrealmul
 {
 if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
 typename fourindex<I,T>::iterator p;
+#ifdef DEBUG
+unsigned long I = SMat_index_1(p->index.indiv.i,p->index.indiv.j);
+unsigned long J = SMat_index_1(p->index.indiv.k,p->index.indiv.l);
+if (I<0 || I>=(unsigned long)NRSMat<T>::nn || J<0 || J>=(unsigned long)NRSMat<T>::nn) laerror("fourindex_dense index out of range in constructor");
+#endif
 for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j,p->index.indiv.k,p->index.indiv.l) = p->elem;
 }
 
@@ -627,8 +632,12 @@ fourindex_dense<twoelectronrealmullikan,T,I>::fourindex_dense<twoelectronrealmul
 {
 if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
 typename fourindex_ext<I,T>::iterator p;
-for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j
-,p->index.indiv.k,p->index.indiv.l) = p->elem;
+#ifdef DEBUG
+unsigned long I = SMat_index_1(p->index.indiv.i,p->index.indiv.j);
+unsigned long J = SMat_index_1(p->index.indiv.k,p->index.indiv.l);
+if (I<0 || I>=(unsigned long)NRSMat<T>::nn || J<0 || J>=(unsigned long)NRSMat<T>::nn) laerror("fourindex_dense index out of range in constructor");
+#endif
+for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j ,p->index.indiv.k,p->index.indiv.l) = p->elem;
 }
 
 
diff --git a/la_traits.h b/la_traits.h
index 2ebd4ca..24651ba 100644
--- a/la_traits.h
+++ b/la_traits.h
@@ -21,9 +21,6 @@ extern "C" {
 }
 #endif
 
-#ifdef _GLIBCPP_NO_TEMPLATE_EXPORT
-# define export 
-#endif
 
 //forward declarations
 template<typename C> class NRVec;
diff --git a/laerror.cc b/laerror.cc
index 977097e..da5dfe8 100644
--- a/laerror.cc
+++ b/laerror.cc
@@ -3,6 +3,8 @@
 #include "laerror.h"
 #include <stdio.h>
 #include <errno.h>
+#include <stdarg.h>
+
 
 #ifdef USE_TRACEBACK
 #include "traceback.h"
@@ -30,7 +32,7 @@ void laerror(const char *s1)
 //stub for f77 blas called from strassen routine
 extern "C" void xerbla_(const char name[6], int *n)
 {
-char msg[128];
+char msg[1024];
 strcpy(msg,"LAPACK or BLAS error in routine ");
 strncat(msg,name,6);
 sprintf(msg+strlen(msg),": illegal value of parameter #%d",*n);
@@ -38,3 +40,30 @@ laerror(msg);
 }
 
 
+//with atlas-cblas another error routine is necessary
+
+extern "C" void ATL_xerbla(int p, char *rout, char *form, ...)
+{
+  char msg0[1024], *msg;
+   va_list argptr;
+   va_start(argptr, form);
+   strcpy(msg0,"ATLAS error\n");
+   msg=msg0+strlen(msg0);
+   if (p) {sprintf(msg, "Parameter %d to routine %s was incorrect\n", p, rout); msg+=strlen(msg);}
+   vsprintf(msg, form, argptr);
+   va_end(argptr);
+   laerror(msg0);
+}
+
+int cblas_errprn(int ierr, int info, char *form, ...)
+{
+  char msg0[1024], *msg;
+   va_list argptr;
+   va_start(argptr, form);
+   sprintf(msg0,"CBLAS error %d %d\n",ierr,info);
+   msg=msg0+strlen(msg0);
+   vsprintf(msg, form, argptr);
+   va_end(argptr);
+   laerror(msg0);
+return 0;
+}
diff --git a/mat.cc b/mat.cc
index 04e7f52..1951229 100644
--- a/mat.cc
+++ b/mat.cc
@@ -779,8 +779,8 @@ void NRMat< complex<double> >::diagmultr(const NRVec< complex<double> > &rhs)
 
 
 
-#ifdef oldversion
 // Mat * Smat, decomposed to nn x Vec * Smat
+//NOTE: dsymm is not appropriate as it works on UNPACKED symmetric matrix
 template<>
 const NRMat<double> 
 NRMat<double>::operator*(const NRSMat<double> &rhs) const
@@ -811,36 +811,6 @@ NRMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
 	return result;
 }
 
-#else
-
-// Mat * Smat
-template<>
-const NRMat<double>
-NRMat<double>::operator*(const NRSMat<double> &rhs) const
-{
-#ifdef DEBUG
-        if (mm != rhs.nrows()) laerror("incompatible dimension in Mat*SMat");
-#endif
-        NRMat<double> result(nn, rhs.ncols());
-	cblas_dsymm(CblasRowMajor, CblasRight, CblasLower, nn, mm, 1., &rhs[0],mm,(*this)[0],mm,0.,result[0],mm);
-        return result;
-}
-
-template<>
-const NRMat< complex<double> >
-NRMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
-{
-#ifdef DEBUG
-        if (mm != rhs.nrows()) laerror("incompatible dimension in Mat*SMat");
-#endif
-        NRMat< complex<double> > result(nn, rhs.ncols());
-	cblas_zhemm(CblasRowMajor, CblasRight, CblasLower, nn, mm, (void *)&CONE, &rhs[0],mm,(*this)[0],mm,(void *)&CZERO,result[0],mm);
-        return result;
-}
-
-
-
-#endif
 
 
 // sum of rows
diff --git a/noncblas.cc b/noncblas.cc
index de8e3a3..ccb6954 100644
--- a/noncblas.cc
+++ b/noncblas.cc
@@ -222,4 +222,5 @@ for(int i=1; i<N; ++i) for(int j=0; j<i; ++j)  {double t=A[j*lda+i]; A[j*lda+i]=
 return INFO;
 }
 
+
 #endif
diff --git a/noncblas.h b/noncblas.h
index 3a6c749..7e292ee 100644
--- a/noncblas.h
+++ b/noncblas.h
@@ -16,7 +16,7 @@
 #define CBLAS_INDEX int
 
 int cblas_errprn(int ierr, int info, char *form, ...);
-
+void ATL_xerbla(int p, char *rout, char *form, ...);
 /*
  * ===========================================================================
  * Prototypes for level 1 BLAS functions (complex are recast as routines)
diff --git a/smat.cc b/smat.cc
index de3b2fd..88395ca 100644
--- a/smat.cc
+++ b/smat.cc
@@ -61,6 +61,7 @@ LA_traits<T>::multiget(NN2,fd,v,dim);
 template <typename T>
 NRSMat<T>::NRSMat(const NRMat<T> &rhs)
 {
+nn=rhs.nrows();
 #ifdef DEBUG
 	if (nn != rhs.ncols()) laerror("attempt to convert non-square Mat to SMat");
 #endif
@@ -148,8 +149,10 @@ void NRSMat<T>::fscanf(FILE *f, const char *format)
  * BLAS specializations for double and complex<double>
  */
 
-#ifdef oldversion
+
+
 // SMat * Mat
+//NOTE: dsymm is not appropriate as it works on UNPACKED symmetric matrix
 template<>
 const NRMat<double> NRSMat<double>::operator*(const NRMat<double> &rhs) const
 {
@@ -178,35 +181,6 @@ NRSMat< complex<double> >::operator*(const NRMat< complex<double> > &rhs) const
 	return result;
 }
 
-#else
-
-// SMat * Mat
-template<>
-const NRMat<double> NRSMat<double>::operator*(const NRMat<double> &rhs) const
-{
-#ifdef DEBUG
-        if (nn != rhs.nrows()) laerror("incompatible dimensions in SMat*Mat");
-#endif
-        NRMat<double> result(nn, rhs.ncols());
-	cblas_dsymm(CblasRowMajor, CblasLeft, CblasLower, nn, rhs.ncols(), 1., (*this),nn, rhs[0],rhs.ncols(), 0.,result[0],rhs.ncols());
-        return result;
-}
-
-
-template<>
-const NRMat< complex<double> >
-NRSMat< complex<double> >::operator*(const NRMat< complex<double> > &rhs) const
-{
-#ifdef DEBUG
-        if (nn != rhs.nrows()) laerror("incompatible dimensions in SMat*Mat");
-#endif
-        NRMat< complex<double> > result(nn, rhs.ncols());
-        cblas_zhemm(CblasRowMajor, CblasLeft, CblasLower, nn, rhs.ncols(), &CONE, (*this),nn, rhs[0],rhs.ncols(), &CZERO,result[0],rhs.ncols());
-        return result;
-}
-
-
-#endif
 
 
 // SMat * SMat
diff --git a/sparsemat.cc b/sparsemat.cc
index ede2dad..8a12b79 100644
--- a/sparsemat.cc
+++ b/sparsemat.cc
@@ -1309,7 +1309,6 @@ return *this;
 }
 
 
-#ifdef _GLIBCPP_NO_TEMPLATE_EXPORT
 
 #define INSTANTIZE(T) \
 template SparseMat<T> & SparseMat<T>::oplusequal(const SparseMat<T> &rhs);\
@@ -1363,4 +1362,3 @@ INSTANTIZE(double)
 INSTANTIZE(complex<double>) //some functions are not OK for hermitean matrices, needs a revision!!!
 
 
-#endif
diff --git a/t.cc b/t.cc
index 66569c5..9cdce7e 100644
--- a/t.cc
+++ b/t.cc
@@ -1105,7 +1105,7 @@ for(int iter=1; iter<100 && norm>1e-8 ; ++iter)
 	}
 }
 
-if(1)
+if(0)
 {
 NRMat<double> a,b;
 cin >>a;
@@ -1115,5 +1115,14 @@ cout <<b*b;
 cout <<(b*b-a).norm();
 }
 
+if(1)
+{
+NRSMat<double> a;
+NRMat<double> b;
+cin >>a>>b;
+cout <<a*b;
+
+}
+
 
 }