float support for nonclass functions

la_traits.h

@@ -456,6 +456,11 @@ generate_traits_smat(NRSMat_from1)
 template<typename C>
 struct LA_traits : LA_traits_aux<C, typename isscalar<C>::scalar_type> {};
 
+
+//precision-dependent default thresholds (the double values are those used so far)
+template<typename T> struct LA_precision {static T svd_threshold() {return (T)1e-12;}};
+template<> struct LA_precision<float> {static float svd_threshold() {return 1e-5f;}};
+
 }//namespace
 
 #endif

mat.cc

@@ -5348,24 +5348,48 @@ NRMat<double> NRMat<double>::inverse() const
 return calcinverse(*this);
 }
 
+template<>
+NRMat<float> NRMat<float>::inverse() const
+{
+return calcinverse(*this);
+}
+
 template<>
 NRMat<std::complex<double> > NRMat<std::complex<double> >::inverse() const
 {
 return calcinverse(*this);
 }       
 
+template<>
+NRMat<std::complex<float> > NRMat<std::complex<float> >::inverse() const
+{
+return calcinverse(*this);
+}       
+
 template<>
 NRMat<double> NRMat<double>::svdinverse(const double thr) const
 {
 return calcsvdinverse(*this,thr,false);
 }
 
+template<>
+NRMat<float> NRMat<float>::svdinverse(const float thr) const
+{
+return calcsvdinverse(*this,thr,false);
+}
+
 template<>
 NRMat<std::complex<double> > NRMat<std::complex<double> >::svdinverse(const double thr) const
 {
 return calcsvdinverse(*this,thr,false);
 }
 
+template<>
+NRMat<std::complex<float> > NRMat<std::complex<float> >::svdinverse(const float thr) const
+{
+return calcsvdinverse(*this,thr,false);
+}
+
 
 
 
@@ -5375,12 +5399,24 @@ NRMat<double> NRMat<double>::absadjoint(const double thr) const
 return calcsvdinverse(*this,thr,true);
 }
 
+template<>
+NRMat<float> NRMat<float>::absadjoint(const float thr) const
+{
+return calcsvdinverse(*this,thr,true);
+}
+
 template<>
 NRMat<std::complex<double> > NRMat<std::complex<double> >::absadjoint(const double thr) const
 {
 return calcsvdinverse(*this,thr,true);
 }
 
+template<>
+NRMat<std::complex<float> > NRMat<std::complex<float> >::absadjoint(const float thr) const
+{
+return calcsvdinverse(*this,thr,true);
+}
+
 
 
 

noncblas.h

@@ -689,5 +689,8 @@ int cblas_errprn(int ierr, int info, char *form, ...);
 int clapack_dgesv(const  CBLAS_ORDER Order, const int N, const int NRHS,
                   double *A, const int lda, int *ipiv,
                   double *B, const int ldb); 
+int clapack_sgesv(const  CBLAS_ORDER Order, const int N, const int NRHS,
+                  float *A, const int lda, int *ipiv,
+                  float *B, const int ldb); 
 #endif
 #endif

nonclass.h

@@ -164,16 +164,18 @@ extern void singular_decomposition(NRMat<T> &a, NRMat<T> *u, NRVec<LA_traits<T>:
 
 declare_la(double)
 declare_la(std::complex<double>)
+declare_la(float)
+declare_la(std::complex<float>)
 
 
 
 //svd_solve without contamination from nullspace, b is n x nrhs, result returns in b, A is destroyed
 template<typename T>
-NRVec<T> svd_solve(NRMat<T> &a, const NRVec<T> &b, double thres=1e-12, int *nullity = NULL)
+NRVec<T> svd_solve(NRMat<T> &a, const NRVec<T> &b, typename LA_traits<T>::normtype thres=LA_precision<typename LA_traits<T>::normtype >::svd_threshold(), int *nullity = NULL)
 {
 if(b.size()!=a.nrows()) laerror("size mismatch in svd_solve");
 if(nullity) *nullity=0;
-NRVec<double> w(a.ncols());
+NRVec<typename LA_traits<T>::normtype > w(a.ncols());
 NRMat<T> u(a.nrows(),a.ncols());
 NRMat<T> vt(a.ncols(),a.ncols());
 singular_decomposition(a,&u,w,&vt,false);
@@ -191,11 +193,11 @@ return res;
 
 
 template<typename T>
-NRMat<T> svd_solve(NRMat<T> &a, const NRMat<T> &b, double thres=1e-12, int *nullity = NULL)
+NRMat<T> svd_solve(NRMat<T> &a, const NRMat<T> &b, typename LA_traits<T>::normtype thres=LA_precision<typename LA_traits<T>::normtype >::svd_threshold(), int *nullity = NULL)
 {
 if(b.nrows()!=a.nrows()) laerror("size mismatch in svd_solve");
 if(nullity) *nullity=0;
-NRVec<double> w(a.ncols());
+NRVec<typename LA_traits<T>::normtype > w(a.ncols());
 NRMat<T> u(a.nrows(),a.ncols());
 NRMat<T> vt(a.ncols(),a.ncols());
 singular_decomposition(a,&u,w,&vt,false);
@@ -221,10 +223,13 @@ return res;
 extern void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
 		NRMat<double> *vl=NULL, NRMat<double> *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
 		NRMat<double> *b=NULL, NRVec<double> *beta=NULL); //this used real storage of eigenvectors like dgeev: for a complex pair (i,i+1) the eigenvector of w_i is Re + i Im = vector i + i*vector i+1
+extern void gdiagonalize(NRMat<float> &a, NRVec<float> &wr, NRVec<float> &wi,
+		NRMat<float> *vl=NULL, NRMat<float> *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
+		NRMat<float> *b=NULL, NRVec<float> *beta=NULL);
 
 template<typename T>
-extern void gdiagonalize(NRMat<T> &a, NRVec< std::complex<double> > &w,
+extern void gdiagonalize(NRMat<T> &a, NRVec<typename LA_traits<T>::complextype > &w,
-		 NRMat< std::complex<double> >*vl=NULL, NRMat< std::complex<double> > *vr=NULL,
+		 NRMat<typename LA_traits<T>::complextype > *vl=NULL, NRMat<typename LA_traits<T>::complextype > *vr=NULL,
 		 const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
 		NRMat<T> *b=NULL, NRVec<T> *beta=NULL); //eigenvectors are stored in complex matrices for T both double and complex
 
@@ -232,9 +237,13 @@ extern void gdiagonalize(NRMat<T> &a, NRVec< std::complex<double> > &w,
 extern void gdiagonalize(NRMat<std::complex<double> > &a, NRVec<double> &wr, NRVec<double> &wi,
                 NRMat<std::complex<double> > *vl=NULL, NRMat<std::complex<double> > *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
                 NRMat<std::complex<double> > *b=NULL, NRVec<std::complex<double> > *beta=NULL);
+extern void gdiagonalize(NRMat<std::complex<float> > &a, NRVec<float> &wr, NRVec<float> &wi,
+                NRMat<std::complex<float> > *vl=NULL, NRMat<std::complex<float> > *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
+                NRMat<std::complex<float> > *b=NULL, NRVec<std::complex<float> > *beta=NULL);
 
 //diagonalization of symmetric real tridiagonal matrix
 extern void diagonalize3(NRVec<double> &d, NRVec<double> &d1, NRMat<double> *v, const bool corder=1, int n0=0);
+extern void diagonalize3(NRVec<float> &d, NRVec<float> &d1, NRMat<float> *v, const bool corder=1, int n0=0);
 
 //complex,real,imaginary parts of various entities
 template<typename T>
@@ -251,6 +260,8 @@ extern const typename LA_traits<T>::complextype complexmatrix (const T&, const T
 //Cholesky decomposition
 extern void cholesky(NRMat<double> &a, bool upper=1);
 extern void cholesky(NRMat<std::complex<double> > &a, bool upper=1);
+extern void cholesky(NRMat<float> &a, bool upper=1);
+extern void cholesky(NRMat<std::complex<float> > &a, bool upper=1);
 
 //inverse by means of linear solve, pass by value to preserve argument intact
 template<typename T>
@@ -268,13 +279,13 @@ const NRMat<T> calcinverse(NRMat<T> a, T *det=NULL)
 
 //inverse by means of SVD , pass by value to preserve argument intact
 template<typename T>
-const NRMat<T> calcsvdinverse(NRMat<T> a, double thr=0, bool do_absadjoint=false)
+const NRMat<T> calcsvdinverse(NRMat<T> a, typename LA_traits<T>::normtype thr=0, bool do_absadjoint=false)
 {
         if(a.nrows()!=a.ncols()) laerror("svdinverse() for non-square matrix");
 	int n=a.nrows();
         a.copyonwrite();
 	NRMat<T> u(n,n),v(n,n);
-	NRVec<double> w(n),s(n);
+	NRVec<typename LA_traits<T>::normtype > w(n),s(n);
 	singular_decomposition(a,&u,w,&v,true);
 	if(do_absadjoint) //absadjoint matrix - missing det(U)*det(V) to be full adjoint - a sign or complex phase
 		{
@@ -348,7 +359,7 @@ return determinant_destroy(a);
 //	the info parameter of gesvx is returned (see man dgesvx)
 //------------------------------------------------------------------------------
 template<class T>
-int linear_solve_x(NRMat<T> &A, T *B, const int rhsCount, const int eqCount, const bool eq, const bool saveA, double *rcond);
+int linear_solve_x(NRMat<T> &A, T *B, const int rhsCount, const int eqCount, const bool eq, const bool saveA, typename LA_traits<T>::normtype *rcond);
 
 
 //------------------------------------------------------------------------------
@@ -364,7 +375,7 @@ int linear_solve_x(NRMat<T> &A, T *B, const int rhsCount, const int eqCount, con
 //	the info parameter of dgesvx is returned (see man dgesvx)
 //------------------------------------------------------------------------------
 template<class T>
-int multiply_by_inverse(NRMat<T> &A, NRMat<T> &B, bool useEq=false, double *rcond=NULL);
+int multiply_by_inverse(NRMat<T> &A, NRMat<T> &B, bool useEq=false, typename LA_traits<T>::normtype *rcond=NULL);
 
 
 //general submatrix, INDEX will typically be NRVec<int> or even int*
@@ -414,6 +425,7 @@ return r;
 
 //auxiliary routine to adjust eigenvectors to guarantee real logarithm
 extern void adjustphases(NRMat<double> &v);
+extern void adjustphases(NRMat<float> &v);
 
 
 //declaration of template interface to cblas routines with full options available
@@ -443,6 +455,24 @@ inline double xdot<double>(int n, const double *x, int incx, const double *y, in
 return cblas_ddot(n,x,incx,y,incy);
 }
 
+template<>
+inline void xcopy<float> (int n, const float *x, int incx, float *y, int incy)
+{
+cblas_scopy(n, x, incx, y, incy);
+}
+
+template<>
+inline void xaxpy<float>(int n, const float &a, const float *x, int incx, float *y, int incy)
+{
+cblas_saxpy(n, a, x, incx, y, incy);
+}
+
+template<>
+inline float xdot<float>(int n, const float *x, int incx, const float *y, int incy)
+{
+return cblas_sdot(n,x,incx,y,incy);
+}
+
 
 
 //debugging aid: reconstruct an explicit matrix from the implicit version
@@ -468,12 +498,14 @@ return r;
 
 extern NRMat<double> realmatrixfunction(NRMat<double> a, double (*f) (double)); //a has to by in fact symmetric
 extern NRMat<std::complex<double> > complexmatrixfunction(NRMat<double> a, double (*fre) (double), double (*fim) (double)); //a has to by in fact symmetric
+extern NRMat<float> realmatrixfunction(NRMat<float> a, float (*f) (float)); //a has to by in fact symmetric
+extern NRMat<std::complex<float> > complexmatrixfunction(NRMat<float> a, float (*fre) (float), float (*fim) (float)); //a has to by in fact symmetric
 
 template<typename T> 
-NRMat<T> matrixfunction(NRSMat<T> a, double (*f) (double)) //of symmetric/hermitian matrix
+NRMat<T> matrixfunction(NRSMat<T> a, typename LA_traits<T>::normtype (*f) (typename LA_traits<T>::normtype)) //of symmetric/hermitian matrix
 {
 	int n = a.nrows();
-	NRVec<double> w(n);
+	NRVec<typename LA_traits<T>::normtype > w(n);
 	NRMat<T> v(n, n);
 	diagonalize(a, w, &v, 0);
 
@@ -488,18 +520,18 @@ NRMat<T> matrixfunction(NRSMat<T> a, double (*f) (double)) //of symmetric/hermit
 
 
 template<typename T>
-extern NRMat<T> matrixfunction(NRMat<T> a, std::complex<double> (*f)(const std::complex<double> &)) //of a general real/complex matrix
+extern NRMat<T> matrixfunction(NRMat<T> a, typename LA_traits<T>::complextype (*f)(const typename LA_traits<T>::complextype &)) //of a general real/complex matrix
 {
 	int n = a.nrows();
-	NRVec<std::complex<double> > w(n);
+	NRVec<typename LA_traits<T>::complextype > w(n);
-	NRMat<std::complex<double> > u(n,n),v(n,n);
+	NRMat<typename LA_traits<T>::complextype > u(n,n),v(n,n);
 
 #ifdef debugmf
-NRMat<std::complex<double> > a0=a;
+NRMat<typename LA_traits<T>::complextype > a0=a;
 #endif
 
         gdiagonalize<T>(a, w, &u, &v, false,n,0,false,NULL,NULL);//a gets destroyed, eigenvectors are rows
-        NRVec< std::complex<double> > z = diagofproduct(u, v, 1, 1);
+        NRVec< typename LA_traits<T>::complextype > z = diagofproduct(u, v, 1, 1);
 
 #ifdef debugmf
 std::cout <<"TEST matrixfunction\n"<<w<<u<<v<<z;
@@ -508,7 +540,7 @@ std::cout <<"TEST matrixfunction2 "<< a0*v.transpose(1) - v.transpose(1)*diagona
 std::cout <<"TEST matrixfunction3 "<< u*v.transpose(1)<<diagonalmatrix(z)<<std::endl;
 #endif
 
-NRVec< std::complex<double> > wz(n);
+NRVec< typename LA_traits<T>::complextype > wz(n);
 for (int i=0; i<a.nrows(); i++) wz[i] = w[i]/z[i];
 
 #ifdef debugmf
@@ -518,7 +550,7 @@ std::cout <<"TEST matrixfunction4 "<< a0<< v.transpose(true)*diagonalmatrix(wz)*
         for (int i=0; i<a.nrows(); i++) w[i] = (*f)(w[i])/z[i];
         u.diagmultl(w);
 
-        NRMat< std::complex<double> > r(n, n);
+        NRMat< typename LA_traits<T>::complextype > r(n, n);
         r.gemm(0.0, v, 'c', u, 'n', 1.0);
 	return (NRMat<T>) r; //convert back to real if applicable by the explicit decomplexifying constructor; it is NOT checked to which accuracy the imaginary part is actually zero
 }
@@ -528,6 +560,8 @@ std::cout <<"TEST matrixfunction4 "<< a0<< v.transpose(true)*diagonalmatrix(wz)*
 
 extern std::complex<double> sqrtinv(const std::complex<double> &);
 extern double sqrtinv(const double);
+extern std::complex<float> sqrtinv(const std::complex<float> &);
+extern float sqrtinv(const float);
 
 //functions on matrices
 inline NRMat<double>  sqrt(const NRSMat<double> &a) { return matrixfunction(a,&std::sqrt); }
@@ -541,6 +575,17 @@ extern NRMat<std::complex<double> > exp0(const NRMat<std::complex<double> > &a);
 extern NRMat<std::complex<double> > copytest(const NRMat<std::complex<double> > &a);
 extern NRMat<double> copytest(const NRMat<double> &a);
 extern NRMat<double> exp0(const NRMat<double> &a);
+inline NRMat<float>  sqrt(const NRSMat<float> &a) { return matrixfunction(a,&std::sqrt); }
+inline NRMat<float>  sqrtinv(const NRSMat<float> &a) { return matrixfunction(a,&sqrtinv); }
+inline NRMat<float>  realsqrt(const NRMat<float> &a) { return realmatrixfunction(a,&std::sqrt); }
+inline NRMat<float>  realsqrtinv(const NRMat<float> &a) { return realmatrixfunction(a,&sqrtinv); }
+inline NRMat<float>  log(const NRSMat<float> &a) { return matrixfunction(a,&std::log); }
+extern NRMat<float> log(const NRMat<float> &a);
+extern NRMat<std::complex<float> > log(const NRMat<std::complex<float> > &a);
+extern NRMat<std::complex<float> > exp0(const NRMat<std::complex<float> > &a);
+extern NRMat<std::complex<float> > copytest(const NRMat<std::complex<float> > &a);
+extern NRMat<float> copytest(const NRMat<float> &a);
+extern NRMat<float> exp0(const NRMat<float> &a);
 
 
 }//namespace

smat.cc

@@ -1324,9 +1324,15 @@ void NRSMat<T>::storesubmatrix(const NRVec<int> &selection, const NRSMat &rhs)
 template<>
 NRSMat<double>  NRSMat<double>::inverse() {return NRSMat<double>(NRMat<double>(*this).inverse());}
 
+template<>
+NRSMat<float>  NRSMat<float>::inverse() {return NRSMat<float>(NRMat<float>(*this).inverse());}
+
 template<>
 NRSMat<std::complex<double> >  NRSMat<std::complex<double> >::inverse() {return NRSMat<std::complex<double> >(NRMat<std::complex<double> >(*this).inverse());}
 
+template<>
+NRSMat<std::complex<float> >  NRSMat<std::complex<float> >::inverse() {return NRSMat<std::complex<float> >(NRMat<std::complex<float> >(*this).inverse());}
+
 
 
 /***************************************************************************//**