float support for sparse matrices

la_traits.h

@@ -468,6 +468,8 @@ template<> struct LA_precision<double>
 	static double eigensolver_eps() {return 1e-6;} //default convergence threshold of davidson and lanczos
 	static double newton_threshold() {return 1e-14;} //convergence of Polynomial::newton
 	static double det_tiny() {return 1e-300;} //determinant underflow guard in linear_solve
+	static double sparse_epsilon() {return 1e-15;} //elements treated as zero in the sparse matrices (SPARSEEPSILON)
+	static double cholesky_epsilon() {return 1e-16;} //pivot threshold of the sparse Cholesky (CHOLESKYEPSILON)
 	};
 template<> struct LA_precision<float>
 	{
@@ -478,6 +480,8 @@ template<> struct LA_precision<float>
 	static float eigensolver_eps() {return 1e-4f;}
 	static float newton_threshold() {return 1e-6f;}
 	static float det_tiny() {return 1e-37f;}
+	static float sparse_epsilon() {return 1e-7f;}
+	static float cholesky_epsilon() {return 1e-8f;}
 	};
 
 }//namespace

nonclass.cc

@@ -2691,6 +2691,18 @@ if(a.issymmetric()) result.setsymmetric();
 return result;
 }
 
+template<>
+const SparseMat<float> realpart<SparseMat< std::complex<float> > >(const SparseMat< std::complex<float> > &a)
+{
+SparseMat<float> result(a.nrows(), a.ncols());
+for(SparseMat< std::complex<float> >::iterator p=a.begin(); p!=a.end(); ++p)
+	{
+	result.add(p->row,p->col,p->elem.real());
+	}
+if(a.issymmetric()) result.setsymmetric();
+return result;
+}
+
 template<>
 const SparseMat<double> imagpart<SparseMat< std::complex<double> > >(const SparseMat< std::complex<double> > &a)
 {
@@ -2703,6 +2715,18 @@ if(a.issymmetric()) result.setsymmetric();
 return result;
 }
 
+template<>
+const SparseMat<float> imagpart<SparseMat< std::complex<float> > >(const SparseMat< std::complex<float> > &a)
+{
+SparseMat<float> result(a.nrows(), a.ncols());
+for(SparseMat< std::complex<float> >::iterator p=a.begin(); p!=a.end(); ++p)
+	{
+	result.add(p->row,p->col,p->elem.imag());
+	}
+if(a.issymmetric()) result.setsymmetric();
+return result;
+}
+
 template<>
 const NRMat<double> realpart<NRMat< std::complex<double> > >(const NRMat< std::complex<double> > &a)
 {
@@ -3044,6 +3068,27 @@ const SparseSMat< std::complex<double> > complexmatrix<SparseSMat<double> >(cons
   return result;
 }
 
+template<>
+const SparseSMat< std::complex<float> > complexmatrix<SparseSMat<float> >(const SparseSMat<float> &re, const SparseSMat<float> &im) {
+  if(re.nrows()!=im.nrows() || re.ncols() != im.ncols()) laerror("incompatible sizes of real and imaginary parts");
+  SparseSMat< std::complex<float> > result(re.nrows(),re.ncols());
+  std::complex<float> tmp;
+
+  SparseSMat<float>::iterator pre(re);
+  for(; pre.notend(); ++pre) {
+    tmp = pre->elem;
+    result.add(pre->row,pre->col,tmp,false);
+  }
+
+  SparseSMat<float>::iterator pim(im);
+  for(; pim.notend(); ++pim) {
+    tmp = std::complex<float>(0,1)*(pim->elem);
+    result.add(pim->row,pim->col,tmp,false);
+  }
+
+  return result;
+}
+
 template<>
 const SparseSMat< std::complex<double> > realmatrix<SparseSMat<double> >(const SparseSMat<double> &re) {
   SparseSMat< std::complex<double> > result(re.nrows(),re.ncols());
@@ -3058,6 +3103,20 @@ const SparseSMat< std::complex<double> > realmatrix<SparseSMat<double> >(const S
   return result;
 }
 
+template<>
+const SparseSMat< std::complex<float> > realmatrix<SparseSMat<float> >(const SparseSMat<float> &re) {
+  SparseSMat< std::complex<float> > result(re.nrows(),re.ncols());
+  std::complex<float> tmp;
+
+  SparseSMat<float>::iterator pre(re);
+  for(; pre.notend(); ++pre) {
+    tmp = pre->elem;
+    result.add(pre->row,pre->col,tmp,false);
+  }
+
+  return result;
+}
+
 template<>
 const SparseSMat< std::complex<double> > imagmatrix<SparseSMat<double> >(const SparseSMat<double> &im) {
   SparseSMat< std::complex<double> > result(im.nrows(),im.ncols());
@@ -3073,6 +3132,21 @@ const SparseSMat< std::complex<double> > imagmatrix<SparseSMat<double> >(const S
   return result;
 }
 
+template<>
+const SparseSMat< std::complex<float> > imagmatrix<SparseSMat<float> >(const SparseSMat<float> &im) {
+  SparseSMat< std::complex<float> > result(im.nrows(),im.ncols());
+  std::complex<float> tmp;
+
+
+  SparseSMat<float>::iterator pim(im);
+  for(; pim.notend(); ++pim) {
+    tmp = std::complex<float>(0,1)*(pim->elem);
+    result.add(pim->row,pim->col,tmp,false);
+  }
+
+  return result;
+}
+
 
 
 

sparsemat.cc

@@ -334,7 +334,7 @@ return nonzero; //number of (in principle) nonzero elements
 
 
 template <class T>
-void SparseMat<T>::simplify(const double sparseepsilon)
+void SparseMat<T>::simplify(const typename LA_traits<T>::normtype sparseepsilon)
 {
 unsigned int n;
 if(!list) return;
@@ -431,7 +431,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
-if(std::abs(elem)<=SPARSEEPSILON) return;
+if(std::abs(elem)<=LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) return;
 if(!count) {count=new int;  *count=1; list=NULL;} //blank new matrix
 else copyonwrite(); //makes also unsort
 add(n,m,elem);
@@ -482,7 +482,7 @@ register matel<T> *l=rhs.list;
 while(l)
 	{
 	if(rhs.symmetric || lower && l->row <=l->col || !lower && l->row >=l->col)
-	    if(std::abs(l->elem)>SPARSEEPSILON) 
+	    if(std::abs(l->elem)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) 
 		add( l->row,l->col,sign=='+'?l->elem:- l->elem);
 	l=l->next;
 	}
@@ -501,14 +501,14 @@ register matel<T> *l=rhs.list;
 if(symmetrize)
 while(l)
         {
-        if(std::abs(l->elem)>SPARSEEPSILON) 
+        if(std::abs(l->elem)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) 
 		{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)
 	{
-	if(std::abs(l->elem)>SPARSEEPSILON) 
+	if(std::abs(l->elem)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) 
 		add( l->row,l->col,l->elem);
 	l=l->next;
 	}
@@ -528,14 +528,14 @@ register matel<T> *l=rhs.list;
 if(symmetrize)
 while(l)
         {
-        if(std::abs(l->elem)>SPARSEEPSILON) 
+        if(std::abs(l->elem)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) 
 		{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)
         {
-        if(std::abs(l->elem)>SPARSEEPSILON) 
+        if(std::abs(l->elem)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) 
 		add( l->row,l->col,- l->elem);
         l=l->next;
         }
@@ -559,7 +559,7 @@ SPMatindex i,j;
 for(i=0;i<nn;++i)
 	for(j=0; j<mm;++j)
 		{register T t(rhs(i,j));
-		if( std::abs(t)>SPARSEEPSILON)
+		if( std::abs(t)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon())
 		 add(i,j,t);
 		}
 }
@@ -598,12 +598,12 @@ if(nn==mm) //square
 else //diagonal of A^TA, assuming it has been simplified (only one entry per position), will be used as preconditioner only anyway
 	while(l)
 		{
-		(*rr)[l->col] += l->elem*l->elem *(l->col!=l->row && symmetric?2.:1.);
+		(*rr)[l->col] += l->elem*l->elem *(T)(l->col!=l->row && symmetric?2:1);
 		l= l->next;
 		}
 if(divide)
 	{
-	for(unsigned int i=0; i<mm; ++i) if((*rr)[i]!=0.) r[i]/=(*rr)[i];
+	for(unsigned int i=0; i<mm; ++i) if((*rr)[i]!=(T)0) r[i]/=(*rr)[i];
 	delete rr;
 	}
 return divide?NULL:&r[0];
@@ -693,7 +693,7 @@ 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 
-if(std::abs(a)<=SPARSEEPSILON) return *this;
+if(std::abs(a)<=LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) return *this;
 SPMatindex i;
 for(i=0;i<nn;++i) add(i,i,a);
 return *this;
@@ -741,7 +741,7 @@ for(i=0;i<nn;++i)
         for(j=0; j<=i;++j)
 		{register T t;
                 if(
-			std::abs(t=rhs(i,j))>SPARSEEPSILON
+			std::abs(t=rhs(i,j))>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()
 			) add(i,j,t);
 		}
 }
@@ -774,7 +774,7 @@ matel<T> *l=list;
 while(l) //include the mirror picture of elements into the list
 	{
 	if(
-		std::abs(l->elem)>SPARSEEPSILON && 
+		std::abs(l->elem)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon() && 
 	   l->row!=l->col) add(l->col,l->row,l->elem); //not OK for complex-hermitean
 	l=l->next;
 	}
@@ -812,12 +812,32 @@ while(l)
 return r;
 }
 
+template<>
+const float SparseMat<float>::dot(const NRMat<float> &rhs) const
+{
+float r=0;
+matel<float> *l=list;
+while(l)
+	{
+	r+= l->elem*rhs[l->row][l->col];
+	if(symmetric&&l->row!=l->col) r+=l->elem*rhs[l->col][l->row];
+	l=l->next;
+	}
+return r;
+}
+
 template<>
 void NRMat<std::complex<double> >::gemm(const std::complex<double>  &beta, const SparseMat<std::complex<double> > &a, const char transa, const NRMat<std::complex<double> > &b, const char transb, const std::complex<double>  &alpha)
 {
 laerror("not implemented yet");
 }
 
+template<>
+void NRMat<std::complex<float> >::gemm(const std::complex<float>  &beta, const SparseMat<std::complex<float> > &a, const char transa, const NRMat<std::complex<float> > &b, const char transb, const std::complex<float>  &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)
@@ -867,6 +887,54 @@ else
 
 }
 
+template<>
+void NRMat<float>::gemm(const float &beta, const SparseMat<float> &a, const char transa, const NRMat<float> &b, const char transb, const float &alpha)
+{
+bool transpa = tolower(transa)!='n'; //not OK for complex
+bool transpb = tolower(transb)!='n'; //not OK for complex
+if(nn!=(int)(transpa?a.ncols():a.nrows()) | mm!= (transpb?b.nrows():b.ncols()) || (int)(transpa?a.nrows():a.ncols()) != (transpb?b.ncols():b.nrows())) laerror("incompatible sizes in gemm");
+copyonwrite();
+if(beta!=(float)0) (*this) *= beta; 
+else memset(v,0,nn*mm*sizeof(float));
+
+matel<float> *l=a.getlist();
+if(alpha==(float)0 || !l) return;
+
+//raw pointers to the full matrices
+const float *bp= b;
+float *p= *this;
+
+int ldb=(transpb?b.ncols():1);
+int bstep=(transpb?1:b.ncols());
+int len=(transpb?b.nrows():b.ncols());
+
+if(a.issymmetric())
+	{
+	while(l)
+	        {
+		cblas_saxpy(len, alpha*l->elem , bp+l->row*bstep, ldb, p+l->col*mm, 1);
+	        if(l->row!=l->col) cblas_saxpy(len, alpha*l->elem , bp+l->col*bstep, ldb, p+l->row*mm, 1);
+	        l=l->next;
+	        }
+	}
+else
+	{
+	if(transpa)
+	while(l)
+	        {
+		cblas_saxpy(len, alpha*l->elem , bp+l->row*bstep, ldb, p+l->col*mm, 1);
+	        l=l->next;
+	        }
+	else
+	while(l)
+	        {
+		cblas_saxpy(len, alpha*l->elem , bp+l->col*bstep, ldb, p+l->row*mm, 1);
+	        l=l->next;
+	        }
+	}
+
+}
+
 
 
 template<class T>
@@ -1036,7 +1104,7 @@ if(transp)
 while(l)
         {
 	register T t=alpha*l->elem;
-	if(std::abs(t)>SPARSEEPSILON) 
+	if(std::abs(t)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) 
 				add( l->col,l->row,t);
         l=l->next;
         }
@@ -1044,7 +1112,7 @@ else
 while(l)
 	{
 	register T t=alpha*l->elem;
-	if(std::abs(t)>SPARSEEPSILON) 
+	if(std::abs(t)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) 
 				add( l->row,l->col,t);
 	l=l->next;
 	}
@@ -1329,7 +1397,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)
-		if(std::abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
+		if(std::abs(tmp=rhs(i,j))>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) add(n0+i,m0+j,tmp);
 return *this;
 }
 
@@ -1345,7 +1413,7 @@ T tmp;
 incsize(rhs.nrows(), rhs.ncols());
 for(SPMatindex i=0; i<(SPMatindex)rhs.nrows(); ++i)
         for(SPMatindex j=0; j<=i; ++j) //symmetric storage: one triangle only
-                if(std::abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
+                if(std::abs(tmp=rhs(i,j))>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) add(n0+i,m0+j,tmp);
 return *this;
 }
 
@@ -1364,7 +1432,7 @@ incsize(rhs.nrows(), rhs.ncols());
 register matel<T> *l=rhs.list;
 while(l)
 	{
-	if(std::abs(l->elem)>SPARSEEPSILON) 
+	if(std::abs(l->elem)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) 
 		add(n0+l->row,m0+l->col,l->elem);
 	l=l->next;
 	}
@@ -1429,8 +1497,12 @@ INSTANTIZE(std::complex<double>) //some functions are not OK for hermitean matri
 //// forced instantization in the corresponding object file
 template class SparseMat<double>;
 template class SparseMat<std::complex<double> >;
+template class SparseMat<float>;
+template class SparseMat<std::complex<float> >;
 template class SemiSparseMat<double>;
 template class SemiSparseMat<std::complex<double> >;
+template class SemiSparseMat<float>;
+template class SemiSparseMat<std::complex<float> >;
 
 
 #define INSTANTIZE(T) \
@@ -1443,10 +1515,8 @@ template void NRVec<T>::gemv(const T beta, const SparseMat<T> &a, const char tra
 
 INSTANTIZE(double)
 INSTANTIZE(std::complex<double>)
-
-//needed by the instantiation of NRVec<float> and NRVec<std::complex<float> >
-template void NRVec<float>::gemv(const float beta, const SparseMat<float> &a, const char trans, const float alpha, const NRVec<float> &x, const bool treat_as_symmetric);
-template void NRVec<std::complex<float> >::gemv(const std::complex<float> beta, const SparseMat<std::complex<float> > &a, const char trans, const std::complex<float> alpha, const NRVec<std::complex<float> > &x, const bool treat_as_symmetric);
+INSTANTIZE(float)
+INSTANTIZE(std::complex<float>)
 
 
 }//namespace

sparsemat.h

@@ -141,7 +141,7 @@ public:
 	unsigned int length() const;
 	void copyonwrite(bool detachonly=false, bool deep=true);
 	void clear() {bool sym=symmetric; resize(nn,mm); symmetric=sym;} //zero matrix of the same dimensions and symmetry
-	void simplify(const double sparseepsilon=SPARSEEPSILON);
+	void simplify(const typename LA_traits<T>::normtype sparseepsilon=LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon());
 	const T trace() const;
 	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;

sparsesmat.cc

@@ -239,7 +239,7 @@ for(SPMatindex i=0; i<nn; ++i)
 		}
 if(divide)
 	{
-	for(unsigned int i=0; i<nn; ++i) if((*rr)[i]!=0.) r[i]/=(*rr)[i];
+	for(unsigned int i=0; i<nn; ++i) if((*rr)[i]!=(T)0) r[i]/=(*rr)[i];
 	delete rr;
 	}
 return divide?NULL:&r[0];
@@ -352,12 +352,18 @@ INSTANTIZE(complex<double>)
 //// forced instantization of functions in the header in the corresponding object file
 template class SparseSMat<double>;
 template class SparseSMat<std::complex<double> >;
+template class SparseSMat<float>;
+template class SparseSMat<std::complex<float> >;
 template SparseSMat<double> NRMat<double>::operator*(const SparseSMat<double> &rhs) const;
+template SparseSMat<float> NRMat<float>::operator*(const SparseSMat<float> &rhs) const;
 template SparseSMat<std::complex<double> > NRMat<std::complex<double> >::operator*(const SparseSMat<std::complex<double> > &rhs) const;
+template SparseSMat<std::complex<float> > NRMat<std::complex<float> >::operator*(const SparseSMat<std::complex<float> > &rhs) const;
 
 /*activate this if needed
 template void SparseSMat<NRMat<double> >::put(int fd, bool dimen, bool transp, bool orcaformat) const;
+template void SparseSMat<NRMat<float> >::put(int fd, bool dimen, bool transp, bool orcaformat) const;
 template void SparseSMat<NRMat<double> >::get(int fd, bool dimen, bool transp, bool orcaformat);
+template void SparseSMat<NRMat<float> >::get(int fd, bool dimen, bool transp, bool orcaformat);
 */
 
 

sparsesmat.h

@@ -200,7 +200,7 @@ count(new int(1))
 v= new std::map<SPMatindex,T> * [nn];
 memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
 int i,j;
-for(i=0; i<nn; ++i) for(j=0; j<=i; ++j) if(std::abs(rhs(i,j))>SPARSEEPSILON) (*this).add(i,j,rhs(i,j),true);
+for(i=0; i<nn; ++i) for(j=0; j<=i; ++j) if(std::abs(rhs(i,j))>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) (*this).add(i,j,rhs(i,j),true);
 }
 
 //from SparseMat: duplicates summed; a symmetric SparseMat (one triangle stored, complex symmetric) is unfolded without conjugation
@@ -227,7 +227,7 @@ count(new int(1))
 v= new std::map<SPMatindex,T> * [nn];
 memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
 int i,j;
-for(i=0; i<nn; ++i) for(j=0; j<mm; ++j) if(std::abs(rhs(i,j))>SPARSEEPSILON) (*this).add(i,j,rhs(i,j),false);
+for(i=0; i<nn; ++i) for(j=0; j<mm; ++j) if(std::abs(rhs(i,j))>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) (*this).add(i,j,rhs(i,j),false);
 }
 
 
@@ -437,7 +437,7 @@ for(SPMatindex i=0; i<nn; ++i) if(v[i])
 	std::list<SPMatindex> l;
 	typename std::map<SPMatindex,T>::iterator p;
 	for(p=v[i]->begin(); p!=v[i]->end(); ++p) 
-	if(std::abs(p->second) < SPARSEEPSILON) l.push_front(p->first); else ++count;
+	if(std::abs(p->second) < LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) l.push_front(p->first); else ++count;
 	typename std::list<SPMatindex>::iterator q;	
 	for(q=l.begin(); q!=l.end(); ++q) v[i]->erase(*q);	
 	if(v[i]->size() == 0) {delete v[i]; v[i]=NULL;}
@@ -548,7 +548,7 @@ for(SPMatindex i=0; i<nn; ++i)
 //std::cout <<"Permuted upper triangle matrix\n"<<r;
 //SparseSMat<T> r0(r);r.copyonwrite();
 
-//perform complex, positive semidefinite Cholesky with thresholding by SPARSEEPSILON
+//perform complex, positive semidefinite Cholesky with thresholding by LA_precision<>::sparse_epsilon()
 SPMatindex i,j,k;
 int rank=0;
 for(k=0; k<nn; ++k)
@@ -557,7 +557,7 @@ for(k=0; k<nn; ++k)
 	typename std::map<SPMatindex,T>::iterator p;
 	p= r.v[k]->find(k);
 	if(p==r.v[k]->end()) continue; //must not break due to the a priori  pivoting
-	if(LA_traits<T>::realpart(p->second) < CHOLESKYEPSILON) continue; //must not break due to the a priori  pivoting
+	if(LA_traits<T>::realpart(p->second) < LA_precision<typename LA_traits<T>::normtype>::cholesky_epsilon()) continue; //must not break due to the a priori  pivoting
 	++rank;
 	typename LA_traits<T>::normtype tmp = std::sqrt(LA_traits<T>::realpart(p->second));
 	p->second = tmp;
@@ -572,7 +572,7 @@ for(k=0; k<nn; ++k)
 		{
 		T akj = LA_traits<T>::conjugate(linek[j]);
 		NRVec<int> linek_used(0,nn);
-		if(std::abs(akj)>SPARSEEPSILON) 
+		if(std::abs(akj)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) 
 			{
 			for(p=r.v[j]->begin(); p!=r.v[j]->end(); ++p)
 				{
@@ -581,7 +581,7 @@ for(k=0; k<nn; ++k)
 				}	
 			//subtract also elements nonzero in line k but non-existent in line j
 			for(i=j; i<nn; ++i) 
-			    if(!linek_used[i] && std::abs(linek[i]) > SPARSEEPSILON)
+			    if(!linek_used[i] && std::abs(linek[i]) > LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon())
 				{
 				(*r.v[j])[i] = -akj * linek[i];
 				}
@@ -611,7 +611,7 @@ for(k=0; k<nn; ++k)
 	typename std::map<SPMatindex,T> *vnew = new typename std::map<SPMatindex,T>;
 	for(p=r.v[k]->begin(); p!=r.v[k]->end(); ++p)
 		{
-        	if(std::abs(p->second) > SPARSEEPSILON) (*vnew)[pivot[p->first]] = p->second;
+        	if(std::abs(p->second) > LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) (*vnew)[pivot[p->first]] = p->second;
 		}
 	delete r.v[k];
 	r.v[k]=vnew;
@@ -636,7 +636,7 @@ for(SPMatindex i=0; i<lhs.size(); ++i)
 		{
 		T x=lhs[i]*(conjugate?LA_traits<T>::conjugate(rhs[j]):rhs[j])*scale;
 		if(i==j) x+=LA_traits<T>::conjugate(x); //diagonal of M + M^\dagger, add() mirrors only off-diagonal elements
-		if(std::abs(x)>SPARSEEPSILON) r.add(i,j,x);
+		if(std::abs(x)>LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) r.add(i,j,x);
 		}
 	}
 return r;