From 5266767a4c71ece7245a6fc80b2c47e45d5c81d2 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Thu, 5 Nov 2026 15:42:29 +0100
Subject: [PATCH] float support for smat

---
 smat.cc | 451 +++++++++++++++++++++++++++++++++++++++++++++++++++++++-
 smat.h  | 231 ++++++++++++++++++++++++++---
 2 files changed, 657 insertions(+), 25 deletions(-)

diff --git a/smat.cc b/smat.cc
index a948a50..af17b97 100644
--- a/smat.cc
+++ b/smat.cc
@@ -167,7 +167,7 @@ const T* NRSMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const
 	if(divide){
 		for(register int i=0; i<nn; i++){
 			const T a = v[(size_t)i*(i+1)/2+i];
-			if(a != 0.) r[i] /= a;
+			if(a != (T)0) r[i] /= a;
 		}
 	}else{
 	        for(register int i=0; i<nn; i++) r[i] = v[(size_t)i*(i+1)/2+i];
@@ -215,6 +215,26 @@ const NRSMat<double> NRSMat<double>::operator-() const {
 	return result;
 }
 
+template <>
+const NRSMat<float> NRSMat<float>::operator-() const {
+	NRSMat<float> result(nn, getlocation());
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		memcpy(result.v, v, NN2*sizeof(float));
+		cblas_sscal(NN2, -1., result.v, 1);
+#ifdef CUDALA
+	}else{
+		cublasScopy(NN2, v, 1, result.v, 1);
+		TEST_CUBLAS("cublasScopy");
+
+		cublasSscal(NN2, -1., result.v, 1);
+		TEST_CUBLAS("cublasSscal");
+	}
+#endif
+	return result;
+}
+
 /***************************************************************************//**
  * implements unary minus operator for this hermitian matrix
  * @return modified copy of this matrix
@@ -240,6 +260,27 @@ const NRSMat<std::complex<double> > NRSMat<std::complex<double> >::operator-() c
 	return result;
 }
 
+template <>
+const NRSMat<std::complex<float> > NRSMat<std::complex<float> >::operator-() const {
+	NRSMat<std::complex<float> > result(nn, getlocation());
+#ifdef CUDALA
+        if(location == cpu) {
+#endif
+		memcpy(result.v, v, NN2*sizeof(std::complex<float>));
+		cblas_cscal(NN2, &CMONEF, result.v, 1);
+
+#ifdef CUDALA
+        }else{
+		cublasCcopy(NN2, (cuComplex*)v, 1, (cuComplex*)result.v, 1);
+		TEST_CUBLAS("cublasCcopy");
+
+		cublasCscal(NN2, CUMONE, (cuComplex*)result.v, 1);
+		TEST_CUBLAS("cublasCscal");
+	}
+#endif
+	return result;
+}
+
 /***************************************************************************//**
  * @return the sum of the diagonal elements
  ******************************************************************************/
@@ -266,6 +307,16 @@ void NRSMat<double>::randomize(const double &x) {
 	}
 }
 
+template<>
+void NRSMat<float>::randomize(const float &x) {
+	NOT_GPU(*this);
+	copyonwrite();
+
+	for(size_t i=0; i<NN2; ++i){
+		v[i] = x*RANDDOUBLESIGNED();
+	}
+}
+
 /***************************************************************************//**
  * Fill this hermitian matrix with pseudorandom numbers generated from uniform 
  * distribution. The real and imaginary parts are generated independently.
@@ -282,6 +333,18 @@ void NRSMat<std::complex<double> >::randomize(const double &x) {
 	}
 }
 
+template<>
+void NRSMat<std::complex<float> >::randomize(const float &x) {
+	copyonwrite();
+	for(register size_t i=0; i<NN2; ++i) v[i].real(x*RANDDOUBLESIGNED());
+	for(register size_t i=0; i<NN2; ++i) v[i].imag(x*RANDDOUBLESIGNED());
+	for(register int i=0; i<nn; ++i){
+		for(register int j=0; j<=i; ++j){
+		       	if(i == j) v[i*(size_t)(i+1)/2+j].imag(0.); //hermitean
+		}
+	}
+}
+
 /***************************************************************************//**
  * routine for formatted output via lawritemat
  * @param[in] file pointer to <tt>FILE</tt> structure representing the output file
@@ -390,6 +453,30 @@ const NRMat<double> NRSMat<double>::operator*(const NRMat<double> &rhs) const {
 	return result;
 }
 
+template<>
+const NRMat<float> NRSMat<float>::operator*(const NRMat<float> &rhs) const {
+#ifdef DEBUG
+	if(nn != rhs.nrows()) laerror("incompatible dimensions in NRMat<float> NRSMat<float>::operator*(const NRMat<float> &)");
+#endif
+	SAME_LOC(*this, rhs);
+	NRMat<float> result(nn, rhs.ncols(), getlocation());
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		for(register int k = 0; k<rhs.ncols(); k++){
+			cblas_sspmv(CblasRowMajor, CblasLower, nn, 1.0, v, rhs[0] + k, rhs.ncols(), 0.0, result[0] + k, rhs.ncols());
+		}
+#ifdef CUDALA
+	}else{
+		for(register int k = 0; k<rhs.ncols(); k++){
+			cublasSspmv('U', nn, 1.0, v, rhs[0] + k, rhs.ncols(), 0.0, result[0] + k, rhs.ncols());
+			TEST_CUBLAS("cublasSspmv");
+		}
+	}
+#endif
+	return result;
+}
+
 
 /***************************************************************************//**
  * multiply this real double-precision symmetric matrix \f$S\f$ stored in packed form
@@ -425,6 +512,34 @@ NRSMat<std::complex<double> >::operator*(const NRMat<std::complex<double> > &rhs
 	return result;
 }
 
+template<>
+const NRMat<std::complex<float> >
+NRSMat<std::complex<float> >::operator*(const NRMat<std::complex<float> > &rhs) const {
+#ifdef DEBUG
+	if (nn != rhs.nrows()) laerror("incompatible dimensions in NRSMat<std::complex<float> >::operator*(const NRMat<std::complex<float> > &)");
+#endif
+	SAME_LOC(*this, rhs);
+	NRMat<std::complex<float> > result(nn, rhs.ncols(), getlocation());
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		for(register int k=0; k<rhs.ncols(); k++){
+			cblas_chpmv(CblasRowMajor, CblasLower, nn, &CONEF, v, rhs[0]+k, rhs.ncols(), &CZEROF, result[0]+k, rhs.ncols());
+		}
+#ifdef CUDALA
+	}else{
+		for(register int k = 0; k<rhs.ncols(); k++){
+			cublasChpmv('U', nn, 
+					CUONE, (cuComplex*)v, (cuComplex*)(rhs[0] + k), rhs.ncols(),
+					CUZERO, (cuComplex*)(result[0] + k), rhs.ncols());
+
+			TEST_CUBLAS("cublasSspmv");
+		}
+	}
+#endif
+	return result;
+}
+
 /***************************************************************************//**
  * multiply this real double-precision symmetric matrix \f$S\f$ stored in packed form
  * with real double-precision symmetric matrix \f$T\f$
@@ -478,6 +593,54 @@ const NRMat<double> NRSMat<double>::operator*(const NRSMat<double> &rhs) const {
 	return result;
 }
 
+template<>
+const NRMat<float> NRSMat<float>::operator*(const NRSMat<float> &rhs) const {
+#ifdef DEBUG
+	if (nn != rhs.nn) laerror("incompatible dimensions in NRMat<float> NRSMat<float>::operator*(const NRSMat<float> &)");
+#endif
+	NRMat<float> result(0.0, nn, nn);
+	float *p, *q;
+
+	p = v;
+	for (int i=0; i<nn;i++) {
+		q = rhs.v;
+		for (int k=0; k<=i; k++) {
+			cblas_saxpy(k+1, *p++, q, 1, result[i], 1);
+			q += k+1;
+		}
+	}
+
+	p = v;
+	for (int i=0; i<nn;i++) {
+		q = rhs.v+1;
+		for (int j=1; j<nn; j++) {
+			result[i][j] += cblas_sdot(i+1<j ? i+1 : j, p, 1, q, 1);
+			q += j+1;
+		}
+		p += i+1;
+	}
+
+	p = v; 
+	q = rhs.v;
+	for (int i=0; i<nn; i++) {
+		cblas_sger(CblasRowMajor, i, i+1, 1., p, 1, q, 1, result, nn);
+		p += i+1;
+		q += i+1;
+	}
+	
+	q = rhs.v+3;
+	for (int j=2; j<nn; j++) {
+		p = v+1;
+		for (int i=1; i<j; i++) {
+			cblas_saxpy(i, *++q, p, 1, result[0]+j, nn);
+			p += i+1;
+		}
+		q += 2;
+	}
+
+	return result;
+}
+
 
 /***************************************************************************//**
  * multiply this complex double-precision symmetric matrix \f$G\f$ stored in packed form
@@ -497,6 +660,19 @@ NRSMat<std::complex<double> >::operator*(const NRSMat<std::complex<double> > &rh
 	return result;
 }
 
+template<>
+const NRMat<std::complex<float> > 
+NRSMat<std::complex<float> >::operator*(const NRSMat<std::complex<float> > &rhs) const {
+#ifdef DEBUG
+	if (nn != rhs.nn) laerror("incompatible dimensions in NRSMat<std::complex<float> >::operator*(const NRSMat<std::complex<float> > &)");
+#endif
+	SAME_LOC(*this, rhs);
+	NRMat<std::complex<float> > result(nn, nn, getlocation());
+	NRMat<std::complex<float> > rhsmat(rhs);
+	result = *this * rhsmat;
+	return result;
+}
+
 
 /***************************************************************************//**
  * compute inner product of this real symmetric matrix \f$A\f$ with given real symmetric matrix \f$B\f$
@@ -524,6 +700,25 @@ const double NRSMat<double>::dot(const NRSMat<double> &rhs) const {
 	return ret;
 }
 
+template<>
+const float NRSMat<float>::dot(const NRSMat<float> &rhs) const {
+	float ret(0.);
+#ifdef DEBUG
+	if (nn != rhs.nn) laerror("incompatible dimensions in float NRSMat<float>::dot(const NRSMat<float> &)");
+#endif
+	SAME_LOC(*this, rhs);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		ret = cblas_sdot(NN2, v, 1, rhs.v, 1);
+#ifdef CUDALA
+	}else{
+		ret = cublasSdot(NN2, v, 1, rhs.v, 1);
+	}
+#endif
+	return ret;
+}
+
 
 /***************************************************************************//**
  * compute inner product of this complex symmetric matrix \f$A\f$ with given complex symmetric matrix \f$B\f$
@@ -554,6 +749,28 @@ const std::complex<double> NRSMat<std::complex<double> >::dot(const NRSMat<std::
 	return dot;
 }
 
+template<>
+const std::complex<float> NRSMat<std::complex<float> >::dot(const NRSMat<std::complex<float> > &rhs) const {
+#ifdef DEBUG
+	if (nn != rhs.nn) laerror("incompatible dimensions in std::complex<float> NRSMat<std::complex<float> >::dot(const NRSMat<std::complex<float> > &)");
+#endif
+	std::complex<float> dot(0., 0.);
+	SAME_LOC(*this, rhs);
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_cdotc_sub(NN2, v, 1, rhs.v, 1, &dot);
+#ifdef CUDALA
+	}else{
+		const cuComplex _dot = cublasCdotc(NN2, (cuComplex*)v, 1, (cuComplex*)(rhs.v), 1);
+		dot = std::complex<float>(cuCreal(_dot), cuCimag(_dot));
+		TEST_CUBLAS("cublasCdotc");
+	}
+#endif
+	return dot;
+}
+
 
 /***************************************************************************//**
  * compute inner product of this real double-precision symmetric matrix \f$S\f$ of order \f$n\f$
@@ -581,6 +798,26 @@ const double NRSMat<double>::dot(const NRVec<double> &rhs) const {
 return ret;
 }
 
+template<>
+const float NRSMat<float>::dot(const NRVec<float> &rhs) const {
+	float ret(0.0);
+#ifdef DEBUG
+	if(NN2 != rhs.nn) laerror("incompatible dimensions in float NRSMat<float>::dot(const NRVec<float> &)");
+#endif
+	SAME_LOC(*this, rhs);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		ret = cblas_sdot(NN2, v, 1, rhs.v, 1);
+#ifdef CUDALA
+	}else{
+		ret = cublasSdot(NN2, v, 1, rhs.v, 1);
+		TEST_CUBLAS("cublasSdot");
+	}
+#endif
+return ret;
+}
+
 
 /***************************************************************************//**
  * compute inner product of this complex double-precision hermitian matrix \f$H\f$ of order \f$n\f$
@@ -610,6 +847,28 @@ NRSMat<std::complex<double> >::dot(const NRVec<std::complex<double> > &rhs) cons
 	return dot;
 }
 
+template<>
+const std::complex<float> 
+NRSMat<std::complex<float> >::dot(const NRVec<std::complex<float> > &rhs) const {
+#ifdef DEBUG
+	if(NN2 != rhs.nn) laerror("incompatible dimensions in std::complex<float>  NRSMat<std::complex<float> >::dot(const NRVec<std::complex<float> > &)");
+#endif
+	std::complex<float> dot(0., 0.);
+	SAME_LOC(*this, rhs);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_cdotc_sub(NN2, v, 1, rhs.v, 1, &dot);
+#ifdef CUDALA
+	}else{
+		const cuComplex _dot = cublasCdotc(NN2, (cuComplex*)v, 1, (cuComplex*)rhs.v, 1);
+		TEST_CUBLAS("cublasCdotc");
+		dot = std::complex<float>(cuCreal(_dot), cuCimag(_dot));
+	}
+#endif
+	return dot;
+}
+
 /***************************************************************************//**
  * compute the Frobenius norm of this real double-precision symmetric  matrix
  * @param[in] scalar subtract this scalar value from the diagonal elements before the norm computation
@@ -645,6 +904,37 @@ const double NRSMat<double>::norm(const double scalar) const {
 	return std::sqrt(sum);
 }
 
+template<>
+const float NRSMat<float>::norm(const float scalar) const {
+	if(!scalar){
+		float ret(0.);
+#ifdef CUDALA
+		if(location == cpu){
+#endif
+			ret = cblas_snrm2(NN2, v, 1);
+#ifdef CUDALA
+		}else{
+			ret = cublasSnrm2(NN2, v, 1);
+			TEST_CUBLAS("cublasSnrm2");
+		}
+#endif
+		return ret;
+	}
+
+	NOT_GPU(*this);
+
+	float sum(0.);
+	int k(0);
+	for(register int i=0; i<nn; ++i){
+		for(register int j=0; j<=i; ++j) {
+			register float tmp = v[k++];
+			if(i == j) tmp -= scalar;
+			sum += tmp*tmp;
+		}
+	}
+	return std::sqrt(sum);
+}
+
 /***************************************************************************//**
  * compute the Frobenius norm of this complex double-precision hermitian  matrix
  * @param[in] scalar subtract this scalar value from the diagonal elements before the norm computation
@@ -680,6 +970,37 @@ const double NRSMat< std::complex<double> >::norm(const std::complex<double> sca
 	return std::sqrt(sum);
 }
 
+template<>
+const float NRSMat< std::complex<float> >::norm(const std::complex<float> scalar) const {
+	if(!(scalar.real()) && !(scalar.imag())){
+		float ret(0.);
+#ifdef CUDALA
+		if(location == cpu){
+#endif
+			ret = cblas_scnrm2(NN2, v, 1);
+#ifdef CUDALA
+		}else{
+			ret = cublasScnrm2(NN2, (cuComplex*)v, 1);
+			TEST_CUBLAS("cublasScnrm2");
+		}
+#endif
+		return ret;
+	}
+
+	int k(0);
+	float sum(0.);
+	std::complex<float> tmp;
+
+	for(register int i=0; i<nn; ++i){
+		for(register int j=0; j<=i; ++j){
+			tmp = v[k++];
+			if (i == j) tmp -= scalar;
+			sum += tmp.real()*tmp.real() + tmp.imag()*tmp.imag();
+		}
+	}
+	return std::sqrt(sum);
+}
+
 
 /***************************************************************************//**
  * symmetric (hermitian) rank-2 update of this matrix stored in packed form
@@ -709,6 +1030,16 @@ void NRSMat<double>::rank2update(const double alpha, const NRVec<double> &x, con
 	cblas_dspr2(CblasRowMajor, CblasLower, nn, alpha, x, 1, y, 1, v);
 }
 
+template<>
+void NRSMat<float>::rank2update(const float alpha, const NRVec<float> &x, const NRVec<float> &y) {
+#ifdef DEBUG
+	if(nn != x.size() || nn != y.size()) laerror("incompatible dimensions in NRSMat<float>::rank2update");
+#endif
+	NOT_GPU(*this); NOT_GPU(x); NOT_GPU(y);
+	copyonwrite();
+	cblas_sspr2(CblasRowMajor, CblasLower, nn, alpha, x, 1, y, 1, v);
+}
+
 template<>
 void NRSMat<std::complex<double> >::rank2update(const std::complex<double> alpha, const NRVec<std::complex<double> > &x, const NRVec<std::complex<double> > &y) {
 #ifdef DEBUG
@@ -719,6 +1050,16 @@ void NRSMat<std::complex<double> >::rank2update(const std::complex<double> alpha
 	cblas_zhpr2(CblasRowMajor, CblasLower, nn, &alpha, x, 1, y, 1, v);
 }
 
+template<>
+void NRSMat<std::complex<float> >::rank2update(const std::complex<float> alpha, const NRVec<std::complex<float> > &x, const NRVec<std::complex<float> > &y) {
+#ifdef DEBUG
+	if(nn != x.size() || nn != y.size()) laerror("incompatible dimensions in NRSMat<std::complex<float> >::rank2update");
+#endif
+	NOT_GPU(*this); NOT_GPU(x); NOT_GPU(y);
+	copyonwrite();
+	cblas_chpr2(CblasRowMajor, CblasLower, nn, &alpha, x, 1, y, 1, v);
+}
+
 
 /***************************************************************************//**
  * for this real double-precision symmetric matrix \f$S\f$ stored in packed form,
@@ -744,6 +1085,25 @@ void NRSMat<double>::axpy(const double alpha, const NRSMat<double> &x) {
 #endif
 }
 
+template<>
+void NRSMat<float>::axpy(const float alpha, const NRSMat<float> &x) {
+#ifdef DEBUG
+	if(nn != x.nn) laerror("incompatible dimensions in void NRSMat<float>::axpy(const float, const NRSMat<float>&)");
+#endif
+	SAME_LOC(*this, x);
+	copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_saxpy(NN2, alpha, x.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasSaxpy(NN2, alpha, x.v, 1, v, 1);
+		TEST_CUBLAS("cublasSaxpy");
+	}
+#endif
+}
+
 
 /***************************************************************************//**
  * for this complex double-precision hermitian matrix \f$H\f$ stored in packed form,
@@ -771,6 +1131,27 @@ void NRSMat<std::complex<double> >::axpy(const std::complex<double> alpha, const
 	
 }
 
+template<>
+void NRSMat<std::complex<float> >::axpy(const std::complex<float> alpha, const NRSMat<std::complex<float> > & x) {
+#ifdef DEBUG
+	if(nn != x.nn) laerror("incompatible dimensions in void NRSMat<std::complex<float> >::axpy(const std::complex<float> , const NRSMat<std::complex<float> >&)");
+#endif
+	SAME_LOC(*this, x);
+	copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_caxpy(NN2, &alpha, x.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		const cuComplex _alpha = make_cuComplex(alpha.real(), alpha.imag());
+		cublasCaxpy(NN2, _alpha, (cuComplex*)x.v, 1, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCaxpy");
+	}
+#endif
+	
+}
+
 /***************************************************************************//**
  * create hermitian matrix \f$H\f$ from given real double-precision symmetric
  * matrix \f$S\f$
@@ -803,6 +1184,32 @@ NRSMat<std::complex<double> >::NRSMat(const NRSMat<double> &rhs, bool imagpart):
 #endif
 }
 
+template<>
+NRSMat<std::complex<float> >::NRSMat(const NRSMat<float> &rhs, bool imagpart): nn(rhs.nrows()), count(new int(1)) {
+	//inconsistent in general case?
+	const int nnp1 = nn*(nn + 1)/2;
+#ifdef CUDALA
+	location = rhs.getlocation();
+	if(location == cpu){
+#endif
+		v = new std::complex<float>[nnp1];
+		memset(v, 0, nnp1*sizeof(std::complex<float>));
+		cblas_scopy(nnp1, &rhs(0, 0), 1, ((float *)v) + (imagpart?1:0), 2);
+#ifdef CUDALA
+	}else{
+		v = (std::complex<float>*) gpualloc(nnp1*sizeof(std::complex<float>));
+
+		std::complex<float> *_val = gpuputcomplex(CZEROF);
+		cublasCcopy(nnp1, (cuComplex*)_val, 0, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCcopy");
+		gpufree(_val);
+
+		cublasScopy(nnp1, (float*)(&rhs(0,0)), 1, ((float*)v) + (imagpart?1:0), 2);
+		TEST_CUBLAS("cublasScopy"); 
+	}
+#endif
+}
+
 
 
 /***************************************************************************//**
@@ -978,13 +1385,51 @@ NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::conjugateme() {
 
 
 
+/***************************************************************************//**
+ * complexify a given real matrix \f$A\f$ (double and float, by BLAS; the complex default constructor zeroes the imaginary parts)
+ * @param[in] rhs matrix \f$A\f$ intended for this operation
+ * @return matrix \f$B\f$ where \f$\Re B=A\f$ and \f$\Im B = 0\f$
+ ******************************************************************************/
+template<>
+NRSMat<std::complex<double> > complexify<double>(const NRSMat<double> &rhs) {
+	NRSMat<std::complex<double> > r(rhs.nrows(), rhs.getlocation());
+#ifdef CUDALA
+	if(rhs.getlocation() == cpu){
+#endif
+		cblas_dcopy(rhs.size(), &(rhs[0]), 1, (double*)(&(r[0])), 2);
+#ifdef CUDALA
+	}else{
+		cublasDcopy(rhs.size(), &(rhs[0]), 1, (double*)(&(r[0])), 2);
+		TEST_CUBLAS("cublasDcopy");//"NRSMat<std::complex<double> > complexify(const NRSMat<double> &)"
+	}
+#endif
+	return r;
+}
+
+template<>
+NRSMat<std::complex<float> > complexify<float>(const NRSMat<float> &rhs) {
+	NRSMat<std::complex<float> > r(rhs.nrows(), rhs.getlocation());
+#ifdef CUDALA
+	if(rhs.getlocation() == cpu){
+#endif
+		cblas_scopy(rhs.size(), &(rhs[0]), 1, (float*)(&(r[0])), 2);
+#ifdef CUDALA
+	}else{
+		cublasScopy(rhs.size(), &(rhs[0]), 1, (float*)(&(r[0])), 2);
+		TEST_CUBLAS("cublasScopy");//"NRSMat<std::complex<float> > complexify(const NRSMat<float> &)"
+	}
+#endif
+	return r;
+}
+
+
 /***************************************************************************//**
  * forced instantization in the corresponding object file
  ******************************************************************************/
 template class NRSMat<double>;
 template class NRSMat<std::complex<double> >;
-//template class NRSMat<float>;
-//template class NRSMat<std::complex<float> >;
+template class NRSMat<float>;
+template class NRSMat<std::complex<float> >;
 template class NRSMat<long long>;
 template class NRSMat<long>;
 template class NRSMat<int>;
diff --git a/smat.h b/smat.h
index 6be29fd..3de96e6 100644
--- a/smat.h
+++ b/smat.h
@@ -387,6 +387,22 @@ inline NRSMat<double> & NRSMat<double>::operator*=(const double &a) {
 	return *this;
 }
 
+template<>
+inline NRSMat<float> & NRSMat<float>::operator*=(const float &a) {
+	copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_sscal(NN2, a, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasSscal(NN2, a, v, 1);
+		TEST_CUBLAS("cublasSscal");//"NRSMat<float>& NRSMat<float>::operator*=(const float &)"
+	}
+#endif
+	return *this;
+}
+
 /***************************************************************************//**
  * multiply this complex symmetric matrix with complex scalar value
  * @param[in] a complex multiplicative factor
@@ -410,6 +426,24 @@ NRSMat<std::complex<double> >::operator*=(const std::complex<double> &a) {
 	return *this;
 }
 
+template<>
+inline NRSMat<std::complex<float> > &
+NRSMat<std::complex<float> >::operator*=(const std::complex<float> &a) {
+	copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_cscal(NN2, &a, v, 1);
+#ifdef CUDALA
+	}else{
+		const cuComplex _a = make_cuComplex(a.real(), a.imag());
+		cublasCscal(NN2, _a, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCscal");//"NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator*=(const std::complex<float> &)"
+	}
+#endif
+	return *this;
+}
+
 
 /***************************************************************************//**
  * multiply this symmetric matrix of general type <code>T</code> stored in packed form
@@ -508,6 +542,27 @@ inline NRSMat<double>& NRSMat<double>::operator+=(const NRSMat<double> & rhs) {
 	return *this;
 }
 
+template<>
+inline NRSMat<float>& NRSMat<float>::operator+=(const NRSMat<float> & rhs) {
+#ifdef DEBUG
+	if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat<float>& NRSMat<float>::operator+=(const NRSMat<float> &)");
+#endif
+	SAME_LOC(*this, rhs);
+	copyonwrite();
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_saxpy(NN2, 1.0, rhs.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasSaxpy(NN2, 1.0, rhs.v, 1, v, 1);
+		TEST_CUBLAS("cublasSaxpy");//" NRSMat<float>& NRSMat<float>::operator+=(const NRSMat<float> &)"
+	}
+#endif
+	return *this;
+}
+
 /***************************************************************************//**
  * add up this complex symmetric matrix with given symmetric matrix
  * @param[in] rhs complex symmetric matrix to be added
@@ -534,6 +589,27 @@ inline NRSMat<std::complex<double> >& NRSMat<std::complex<double> >::operator+=(
 	return *this;
 }
 
+template<>
+inline NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator+=(const NRSMat<std::complex<float> > & rhs) {
+#ifdef DEBUG
+	if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator+=(const NRSMat<std::complex<float> > &)");
+#endif
+	SAME_LOC(*this, rhs);
+	copyonwrite();
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_caxpy(NN2, &CONEF, rhs.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasCaxpy(NN2, CUONE, (cuComplex*)(rhs.v), 1, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCaxpy");//"NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator+=(const NRSMat<std::complex<float> > &)"
+	}
+#endif
+	return *this;
+}
+
 /***************************************************************************//**
  * add up this symmetric matrix of general type <code>T</code> with given symmetric matrix
  * @param[in] rhs complex matrix of general type <code>T</code> to be added
@@ -578,6 +654,27 @@ inline NRSMat<double>& NRSMat<double>::operator-=(const NRSMat<double>& rhs) {
 	return *this;
 }
 
+template<>
+inline NRSMat<float>& NRSMat<float>::operator-=(const NRSMat<float>& rhs) {
+#ifdef DEBUG
+	if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat<float>& NRSMat<float>::operator-=(const NRSMat<float> &)");
+#endif
+	SAME_LOC(*this, rhs);
+	copyonwrite();
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_saxpy(NN2, -1.0, rhs.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasSaxpy(NN2, -1.0, rhs.v, 1, v, 1);
+		TEST_CUBLAS("cublasSaxpy");//" NRSMat<float>& NRSMat<float>::operator-=(const NRSMat<float> &)"
+	}
+#endif
+	return *this;
+}
+
 /***************************************************************************//**
  * subtracts given complex symmetric matrix from this complex symmetric matrix
  * @param[in] rhs complex symmetric matrix to be subtracted
@@ -604,6 +701,27 @@ inline NRSMat<std::complex<double> >& NRSMat<std::complex<double> >::operator-=(
 	return *this;
 }
 
+template<>
+inline NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator-=(const NRSMat<std::complex<float> >& rhs) {
+#ifdef DEBUG
+	if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator-=(const NRSMat<std::complex<float> > &)");
+#endif
+	SAME_LOC(*this, rhs);
+	copyonwrite();
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_caxpy(NN2, &CMONEF, rhs.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasCaxpy(NN2, CUMONE, (cuComplex*)(rhs.v), 1, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCaxpy");//"NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator-=(const NRSMat<std::complex<float> > &)"
+	}
+#endif
+	return *this;
+}
+
 /***************************************************************************//**
  * subtracts given symmetric matrix of general type <code>T</code> from this symmetric matrix of type <code>T</code>
  * @param[in] rhs symmetric matrix of general type <code>T</code> to be subtracted
@@ -871,6 +989,24 @@ inline const double NRSMat<double>::amax() const {
 	return ret;
 }
 
+template<>
+inline const float NRSMat<float>::amax() const {
+	float ret(0.0);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		ret = v[cblas_isamax(NN2, v, 1)];
+#ifdef CUDALA
+	}else{
+		const int pozice = cublasIsamax(NN2, v, 1) ;
+		TEST_CUBLAS("cublasIsamax");//"float NRSMat<float>::amax()"
+
+		gpuget(1, sizeof(float), v + pozice, &ret);
+	}
+#endif
+	return ret;
+}
+
 /***************************************************************************//**
  * for this real symmetric matrix \f$A\f$, determine the 
  * first element with smallest absolute value
@@ -901,6 +1037,31 @@ inline const double NRSMat<double>::amin() const {
 	return ret;
 }
 
+template<>
+inline const float NRSMat<float>::amin() const {
+	float ret(0.0);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		// idamin seems not to be supported
+		float val(0.0);
+		int index(-1);
+		ret = std::numeric_limits<float>::max();
+		for(register size_t i = 0; i < NN2; i++){
+			val = std::abs(v[i]);
+			if(val < ret){ index = i; ret = val; }
+		}
+		ret = v[index];
+#ifdef CUDALA
+	}else{
+		const int pozice = cublasIsamin(nn, v, 1);
+		TEST_CUBLAS("cublasIsamin");//"float NRSMat<float>::amin()"
+		gpuget(1, sizeof(float), v + pozice, &ret);
+	}
+#endif
+	return ret;
+}
+
 /***************************************************************************//**
  * for this complex symmetric matrix \f$A\f$, determine the 
  * first element with largest "absolute value"
@@ -923,6 +1084,23 @@ inline const std::complex<double> NRSMat< std::complex<double> >::amax() const{
 	return ret;
 }
 
+template<>
+inline const std::complex<float> NRSMat< std::complex<float> >::amax() const{
+	std::complex<float> ret(0., 0.);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		ret = v[cblas_icamax(NN2, v, 1)];
+#ifdef CUDALA
+	}else{
+		const int pozice = cublasIcamax(NN2, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasIcamax");//"std::complex<float> NRSMat<std::complex<float> >::amax()"
+		gpuget(1, sizeof(std::complex<float>), v + pozice, &ret);
+	}
+#endif
+	return ret;
+}
+
 /***************************************************************************//**
  * for this complex symmetric matrix \f$A\f$, determine the 
  * first element with smallest "absolute value"
@@ -956,6 +1134,34 @@ inline const std::complex<double> NRSMat<std::complex<double> >::amin() const{
 	return ret;
 }
 
+template<>
+inline const std::complex<float> NRSMat<std::complex<float> >::amin() const{
+	std::complex<float> ret(0., 0.);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		// izamin seems not to be supported
+		int index(0);
+		float val(0.0), min_val(0.0);
+		std::complex<float> z_val(0.0, 0.0);
+
+		min_val = std::numeric_limits<float>::max();
+		for(register size_t i = 0; i < NN2; i++){
+			z_val = v[i];
+			val = std::abs(z_val.real()) + std::abs(z_val.imag());
+			if(val < min_val){ index = i; min_val = val; }	
+		}
+		ret = v[index];
+#ifdef CUDALA
+	}else{
+		const int pozice = cublasIcamin(nn, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasIcamin");//"std::complex<float> NRSMat<std::complex<float> >::amin()"
+		gpuget(1, sizeof(std::complex<float>), v + pozice, &ret);
+	}
+#endif
+	return ret;
+}
+
 /***************************************************************************//**
  * @return pointer of general type T to the underlying data structure
  ******************************************************************************/
@@ -1215,28 +1421,9 @@ NRSMat<std::complex<T> > complexify(const NRSMat<T> &rhs) {
 }
 
 
-/***************************************************************************//**
- * complexify a given double-precision real matrix \f$A\f$
- * @param[in] rhs matrix \f$A\f$ intended for this operation
- * @return matrix \f$B\f$ where \f$\Re B=A\f$ and \f$\Im B = 0\f$
- ******************************************************************************/
-/*
-template<>
-NRSMat<std::complex<double> > complexify(const NRSMat<double> &rhs) {
-	NRSMat<std::complex<double> > r(rhs.nrows(), rhs.getlocation());
-#ifdef CUDALA
-	if(rhs.getlocation() == cpu){
-#endif
-		cblas_dcopy(rhs.size(), &(rhs[0]), 1, (double*)(&(r[0])), 2);
-#ifdef CUDALA
-	}else{
-		cublasDcopy(rhs.size(), &(rhs[0]), 1, (double*)(&(r[0])), 2);
-		TEST_CUBLAS("cublasDcopy");//"NRSMat<std::complex<double> > complexify(const NRSMat<double> &)"
-	}
-#endif
-	return r;
-}
-*/
+//BLAS versions for double and float are in smat.cc (a full specialization defined in the header would be multiply defined)
+template<> NRSMat<std::complex<double> > complexify<double>(const NRSMat<double> &rhs);
+template<> NRSMat<std::complex<float> > complexify<float>(const NRSMat<float> &rhs);
 /***************************************************************************//**
  * output operator
  * @param[in,out] s output stream