diff --git a/nonclass.cc b/nonclass.cc
index 9f845cc..d6a3dae 100644
--- a/nonclass.cc
+++ b/nonclass.cc
@@ -37,6 +37,8 @@ template void lawritemat(FILE *file,const T *a,int r,int c,const char *form0, \
 		int nodim,int modulo, int issym);
 INSTANTIZE(double)
 INSTANTIZE(std::complex<double>)
+INSTANTIZE(float)
+INSTANTIZE(std::complex<float>)
 INSTANTIZE(int)
 INSTANTIZE(short)
 INSTANTIZE(char)
diff --git a/sparsemat.cc b/sparsemat.cc
index 342eb05..90162ed 100644
--- a/sparsemat.cc
+++ b/sparsemat.cc
@@ -1444,5 +1444,9 @@ 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);
+
 
 }//namespace
diff --git a/vec.cc b/vec.cc
index ede85d1..e5579ae 100644
--- a/vec.cc
+++ b/vec.cc
@@ -106,6 +106,23 @@ const NRVec<double> NRVec<double>::operator-() const {
 	return result;
 }
 
+template<>
+const NRVec<float> NRVec<float>::operator-() const {
+	NRVec<float> result(*this);
+	result.copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_sscal(nn, -1.0, result.v, 1);
+#ifdef CUDALA
+	}else{
+		cublasSscal(nn, -1.0, result.v, 1);
+		TEST_CUBLAS("cublasSscal");
+	}
+#endif
+	return result;
+}
+
 
 
 /***************************************************************************//**
@@ -129,6 +146,23 @@ const NRVec<std::complex<double> > NRVec<std::complex<double> >::operator-() con
 	return result;
 }
 
+template<>
+const NRVec<std::complex<float> > NRVec<std::complex<float> >::operator-() const {
+	NRVec<std::complex<float> > result(*this);
+	result.copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_csscal(nn, -1.0, result.v, 1);
+#ifdef CUDALA
+	}else{
+		cublasCsscal(nn, -1.0, (cuComplex*)result.v, 1);
+		TEST_CUBLAS("cublasCsscal");
+	}
+#endif
+	return result;
+}
+
 /***************************************************************************//**
  * unary minus operator for vector of general type
  * @return the modified vector
@@ -199,6 +233,16 @@ void NRVec<double>::randomize(const double &x){
 	}
 }
 
+template<>
+void NRVec<float>::randomize(const float &x){
+	NOT_GPU(*this);
+	copyonwrite();
+
+	for(register int i=0; i<nn; ++i){
+		v[i] = x*RANDDOUBLESIGNED();
+	}
+}
+
 /***************************************************************************//**
  * fill the complex vector with pseudorandom numbers generated using uniform distribution
  * the real and imaginary parts are generated independently
@@ -217,6 +261,16 @@ void NRVec<std::complex<double> >::randomize(const double &x) {
 	}
 }
 
+template<>
+void NRVec<std::complex<float> >::randomize(const float &x) {
+	NOT_GPU(*this);
+	copyonwrite();
+
+	for(register int i=0; i<nn; ++i){
+		v[i] = std::complex<float>(x*RANDDOUBLESIGNED(), x*RANDDOUBLESIGNED());
+	}
+}
+
 /***************************************************************************//**
  * constructor creating complex vector from a real one
  * @param[in] rhs the real vector being converted into the complex one
@@ -252,6 +306,31 @@ NRVec<std::complex<double> >::NRVec(const NRVec<double> &rhs, bool imagpart): nn
 #endif
 }
 
+template<>
+NRVec<std::complex<float> >::NRVec(const NRVec<float> &rhs, bool imagpart): nn(rhs.size()){
+
+	 count = new int;
+	*count = 1;
+#ifdef CUDALA
+	location = rhs.getlocation();
+	if(location == cpu){
+#endif
+		v = (std::complex<float>*)new std::complex<float>[nn];
+		memset(v, 0, nn*sizeof(std::complex<float>));
+		cblas_scopy(nn, &rhs[0], 1, ((float *)v) + (imagpart?1:0), 2);
+#ifdef CUDALA
+	}else{
+		v = (std::complex<float>*) gpualloc(nn*sizeof(std::complex<float>));
+
+		cublasCscal(nn, CUZERO, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCscal");
+
+		cublasScopy(nn, &rhs[0], 1, ((float *)v) + (imagpart?1:0), 2);
+		TEST_CUBLAS("cublasScopy");
+	}
+#endif
+}
+
 /***************************************************************************//**
  * perform the <b>axpy</b> operation on the current real vector \f$\vec{v}\f$, i.e.
  * \f[ \vec{v} \leftarrow \vec{v} + \alpha\vec{x} \f]
@@ -278,6 +357,26 @@ void NRVec<double>::axpy(const double alpha, const NRVec<double> &x) {
 #endif
 }
 
+template<>
+void NRVec<float>::axpy(const float alpha, const NRVec<float> &x) {
+#ifdef DEBUG
+	if (nn != x.nn) laerror("incompatible vectors");
+#endif
+	SAME_LOC(*this, x);
+	copyonwrite();
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_saxpy(nn, alpha, x.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasSaxpy(nn, alpha, x.v, 1, v, 1);
+		TEST_CUBLAS("cublasSaxpy");
+	}
+#endif
+}
+
 /***************************************************************************//**
  * perform the <b>axpy</b> operation on the current complex vector \f$\vec{v}\f$, i.e.
  * \f[ \vec{v} \leftarrow \vec{v} + \alpha\vec{x} \f]
@@ -306,6 +405,28 @@ void NRVec<std::complex<double> >::axpy(const std::complex<double> alpha, const
 
 }
 
+template<>
+void NRVec<std::complex<float> >::axpy(const std::complex<float> alpha, const NRVec<std::complex<float> > &x){
+#ifdef DEBUG
+	if (nn != x.nn) laerror("incompatible vectors");
+#endif
+	SAME_LOC(*this, x);
+	copyonwrite();
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_caxpy(nn, &alpha, x.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		const cuComplex _alpha = make_cuComplex(alpha.real(), alpha.imag());
+		cublasCaxpy(nn, _alpha, (cuComplex*)x.v, 1, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCaxpy");
+	}
+#endif
+
+}
+
 /***************************************************************************//**
  * perform the <b>axpy</b> operation on the current real vector \f$\vec{v}\f$, i.e.
  * \f[ \vec{v} \leftarrow \vec{v} + \alpha\vec{x} \f]
@@ -321,6 +442,14 @@ void NRVec<double>::axpy(const double alpha, const double *x, const int stride){
 	cblas_daxpy(nn, alpha, x, stride, v, 1);
 }
 
+template<>
+void NRVec<float>::axpy(const float alpha, const float *x, const int stride){
+	NOT_GPU(*this);
+
+	copyonwrite();
+	cblas_saxpy(nn, alpha, x, stride, v, 1);
+}
+
 /***************************************************************************//**
  * perform the <b>axpy</b> operation on the current complex vector \f$\vec{v}\f$, i.e.
  * \f[ \vec{v} \leftarrow \vec{v} + \alpha\vec{x} \f]
@@ -336,6 +465,14 @@ void NRVec<std::complex<double> >::axpy(const std::complex<double> alpha, const
 	cblas_zaxpy(nn, &alpha, x, stride, v, 1);
 }
 
+template<>
+void NRVec<std::complex<float> >::axpy(const std::complex<float> alpha, const std::complex<float> *x, const int stride){
+	NOT_GPU(*this);
+
+	copyonwrite();
+	cblas_caxpy(nn, &alpha, x, stride, v, 1);
+}
+
 /***************************************************************************//**
  * assign real scalar value to every element of the current vector
  * @param[in] a scalar value to be assigned
@@ -357,6 +494,22 @@ copyonwrite();
 	return *this;
 }
 
+template<>
+NRVec<float>& NRVec<float>::operator=(const float &a){
+copyonwrite();
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_scopy(nn, &a, 0, v, 1);
+#ifdef CUDALA
+	}else{
+		smart_gpu_set(nn, (float)0, v);
+	}
+#endif	
+	return *this;
+}
+
 /***************************************************************************//**
  * assign complex scalar value to every element of the current vector
  * @param[in] a scalar value to be assigned
@@ -378,6 +531,22 @@ copyonwrite();
 	return *this;
 }
 
+template<>
+NRVec<std::complex<float> >& NRVec<std::complex<float> >::operator=(const std::complex<float> &a){
+copyonwrite();
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_ccopy(nn, &a, 0, v, 1);
+#ifdef CUDALA
+	}else{
+		smart_gpu_set(nn, (std::complex<float>)0, v);
+	}
+#endif	
+	return *this;
+}
+
 /***************************************************************************//**
  * assign scalar value to every element of the current vector of general type <code>T</code>
  * @param[in] a scalar value to be assigned
@@ -433,6 +602,38 @@ NRVec<double>& NRVec<double>::normalize(double *norm){
 	return *this;
 }
 
+template<>
+NRVec<float>& NRVec<float>::normalize(float *norm){
+	float tmp(0.0);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		tmp = cblas_snrm2(nn, v, 1);
+		if(norm) *norm = tmp;
+		#ifdef DEBUG
+			if(!tmp) laerror("attempt to normalize zero vector");
+		#endif
+		copyonwrite();
+		tmp = 1.0 / tmp;
+		cblas_sscal(nn, tmp, v, 1);
+#ifdef CUDALA
+	}else{
+		tmp = cublasSnrm2(nn, v, 1);
+		TEST_CUBLAS("cublasSnrm2");
+
+		if(norm) *norm = tmp;
+		#ifdef DEBUG
+			if(!tmp) laerror("attempt to normalize zero vector");
+		#endif
+		copyonwrite();
+		tmp = 1.0 / tmp;
+		cublasSscal(nn, tmp, v, 1);
+		TEST_CUBLAS("cublasSscal");
+	}
+#endif
+	return *this;
+}
+
 /***************************************************************************//**
  * normalize current complex vector (in the Euclidean norm)
  * @param[in] norm if not NULL, the norm of this vector is stored into *norm
@@ -471,6 +672,39 @@ NRVec<std::complex<double> > & NRVec<std::complex<double> >::normalize(double *n
 	return *this;
 }
 
+template<>
+NRVec<std::complex<float> > & NRVec<std::complex<float> >::normalize(float *norm){
+	float tmp(0.0);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		tmp = cblas_scnrm2(nn, v, 1);
+		if(norm) *norm = tmp;
+		#ifdef DEBUG
+		if(tmp == 0.0) laerror("attempt to normalize zero vector");
+		#endif
+		copyonwrite();
+		tmp = 1.0 / tmp;
+		cblas_csscal(nn, tmp, v, 1);
+#ifdef CUDALA
+	}else{
+		tmp = cublasScnrm2(nn, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasScnrm2");
+
+		if(norm) *norm = tmp;
+		#ifdef DEBUG
+			if(tmp == 0.0) laerror("attempt to normalize zero vector");
+		#endif
+		copyonwrite();
+
+		tmp = 1.0 / tmp;
+		cublasCsscal(nn, tmp, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCsscal");
+	}
+#endif		
+	return *this;
+}
+
 /***************************************************************************//**
  * perform the \b gemv operation on this real vector \f$\vec{y}\f$, i.e.
  * \f[\vec{y}\leftarrow \alpha\operatorname{op}(A)\cdot\vec{x}+\beta\vec{y}\f]
@@ -503,6 +737,28 @@ void NRVec<double>::gemv(const double beta, const NRMat<double> &A,
 #endif
 }
 
+template<>
+void NRVec<float>::gemv(const float beta, const NRMat<float> &A, 
+		const char trans, const float alpha, const NRVec &x) {
+#ifdef DEBUG
+	if((tolower(trans) == 'n'?A.ncols():A.nrows()) != x.size()){ laerror("incompatible vectors"); }
+	if((tolower(trans) == 'n'?A.nrows():A.ncols()) != nn){ laerror("incompatible result vector in gemv"); }
+#endif
+	SAME_LOC3(*this, x, A);
+	copyonwrite();
+
+#ifdef CUDALA
+	if(location==cpu){
+#endif
+		cblas_sgemv(CblasRowMajor, (tolower(trans)=='n' ? CblasNoTrans:CblasTrans), A.nrows(), A.ncols(), alpha, A, A.ncols(), x.v, 1, beta, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasSgemv((tolower(trans)=='n'?'T':'N'), A.ncols(), A.nrows(), alpha, A, A.ncols(), x.v, 1, beta, v, 1);
+		TEST_CUBLAS("cublasSgemv");
+	}
+#endif
+}
+
 /***************************************************************************//**
  * perform the \b gemv operation on this complex vector \f$\vec{y}\f$, i.e.
  * \f[\vec{y}\leftarrow \alpha\operatorname{op}(A)\cdot\vec{x}+\beta\vec{y}\f]
@@ -541,6 +797,34 @@ void NRVec<std::complex<double> >::gemv(const double beta, const NRMat<double> &
 #endif
 }
 
+template<>
+void NRVec<std::complex<float> >::gemv(const float beta, const NRMat<float> &A,
+                const char trans, const float alpha, const NRVec<std::complex<float> > &x) {
+#ifdef DEBUG
+	if ((tolower(trans) == 'n'?A.ncols():A.nrows()) != x.size()){ laerror("incompatible vectors"); }
+	if((tolower(trans) == 'n'?A.nrows():A.ncols()) != nn){ laerror("incompatible result vector in gemv"); }
+#endif
+	SAME_LOC3(*this, x, A);
+	copyonwrite();
+
+#ifdef CUDALA
+	if(location==cpu){
+#endif
+		cblas_sgemv(CblasRowMajor, (tolower(trans)=='n'?CblasNoTrans:CblasTrans),
+				A.nrows(), A.ncols(), alpha, A, A.ncols(), (float *)x.v, 2, beta, (float *)v, 2);
+		cblas_sgemv(CblasRowMajor, (tolower(trans)=='n'?CblasNoTrans:CblasTrans),
+				A.nrows(), A.ncols(), alpha, A, A.ncols(), ((float *)x.v) + 1, 2, beta, ((float *)v)+1, 2);
+#ifdef CUDALA
+	}else{
+		cublasSgemv((tolower(trans)=='n'?'T':'N'), A.ncols(), A.nrows(), alpha, A, A.ncols(), (float*)(x.v), 2, beta, (float *)v, 2);
+		TEST_CUBLAS("cublasSgemv");
+
+		cublasSgemv((tolower(trans)=='n'?'T':'N'), A.ncols(), A.nrows(), alpha, A, A.ncols(), ((float *)x.v) + 1, 2, beta, ((float *)v)+1, 2);
+		TEST_CUBLAS("cublasSgemv");
+	}
+#endif
+}
+
 
 /***************************************************************************//**
  * perform the \b gemv operation on this complex vector \f$\vec{y}\f$, i.e.
@@ -581,6 +865,35 @@ void NRVec<std::complex<double> >::gemv(const std::complex<double> beta,
 #endif
 }
 
+template<>
+void NRVec<std::complex<float> >::gemv(const std::complex<float> beta,
+		const NRMat<std::complex<float> > &A, const char trans, 
+		const std::complex<float> alpha, const NRVec<std::complex<float> > &x) {
+#ifdef DEBUG
+	if ((tolower(trans) == 'n'?A.ncols():A.nrows()) != x.size()){ laerror("incompatible vectors"); }
+	if((tolower(trans) == 'n'?A.nrows():A.ncols()) != nn){ laerror("incompatible result vector in gemv"); }
+#endif
+	SAME_LOC3(*this, x, A);
+	copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_cgemv(CblasRowMajor, (tolower(trans)=='n' ? CblasNoTrans : (tolower(trans)=='c'?CblasConjTrans:CblasTrans)), 
+				A.nrows(), A.ncols(), &alpha, A, A.ncols(), x.v, 1, &beta, v, 1);
+#ifdef CUDALA
+	}else{
+		if(tolower(trans)=='c') laerror("gemv with conjugate transposition not implemented on GPU");
+
+		const cuComplex _alpha = make_cuComplex(alpha.real(), alpha.imag());
+		const cuComplex  _beta = make_cuComplex(beta.real(), beta.imag());
+
+		cublasCgemv((tolower(trans)=='n'?'T':'N'), A.ncols(), A.nrows(),
+				_alpha, (cuComplex*)(A[0]), A.ncols(), (cuComplex*)(x.v), 1, _beta, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCgemv");
+	}
+#endif
+}
+
 /***************************************************************************//**
  * perform the \b gemv operation on this real vector \f$\vec{y}\f$, i.e.
  * \f[\vec{y}\leftarrow \alpha\operatorname{op}(A)\cdot\vec{x}+\beta\vec{y}\f]
@@ -613,6 +926,28 @@ void NRVec<double>::gemv(const double beta, const NRSMat<double> &A,
 #endif
 }
 
+template<>
+void NRVec<float>::gemv(const float beta, const NRSMat<float> &A,
+		const char trans, const float alpha, const NRVec &x) {
+#ifdef DEBUG
+	if(A.ncols() != x.size()){ laerror("incompatible dimensions"); }
+	if(A.nrows() != nn){ laerror("incompatible result vector in gemv"); }
+#endif
+	SAME_LOC3(*this, A, x);
+	copyonwrite();
+
+#ifdef CUDALA
+	if(location==cpu){
+#endif
+		cblas_sspmv(CblasRowMajor, CblasLower, A.ncols(), alpha, A, x.v, 1, beta, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasSspmv('U', A.ncols(), alpha, A, x.v, 1, beta, v, 1);
+		TEST_CUBLAS("cublasSspmv");
+	}
+#endif
+}
+
 
 /***************************************************************************//**
  * perform the \c gemv operation on this complex vector \f$\vec{y}\f$, i.e.
@@ -650,6 +985,32 @@ void NRVec<std::complex<double> >::gemv(const double beta, const NRSMat<double>
 #endif
 }
 
+template<>
+void NRVec<std::complex<float> >::gemv(const float beta, const NRSMat<float> &A,
+		const char trans, const float alpha, const NRVec<std::complex<float> > &x) {
+#ifdef DEBUG
+	if(A.ncols() != x.size()){ laerror("incompatible dimensions"); }
+	if(A.nrows() != nn){ laerror("incompatible result vector in gemv"); }
+#endif
+	SAME_LOC3(*this, A, x);
+	copyonwrite();
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_sspmv(CblasRowMajor, CblasLower, A.ncols(), alpha, A, (float *)x.v, 2, beta, (float *)v, 2);
+		cblas_sspmv(CblasRowMajor, CblasLower, A.ncols(), alpha, A, ((float *)x.v) + 1, 2, beta, ((float *)v) + 1, 2);
+#ifdef CUDALA
+	}else{
+		cublasSspmv('U', A.ncols(), alpha, A, (float*)(x.v), 2, beta, (float*)v, 2);
+		TEST_CUBLAS("cublasSspmv");
+		
+		cublasSspmv('U', A.ncols(), alpha, A, ((float*)(x.v)) + 1, 2, beta, ((float*)v) + 1, 2);
+		TEST_CUBLAS("cublasSspmv");
+	}
+#endif
+}
+
 /***************************************************************************//**
  * perform the \b gemv operation on this complex vector \f$\vec{y}\f$, i.e.
  * \f[\vec{y}\leftarrow \alpha\operatorname{op}(A)\cdot\vec{x}+\beta\vec{y}\f]
@@ -689,6 +1050,35 @@ void NRVec<std::complex<double> >::gemv(const std::complex<double> beta,
 #endif
 }
 
+template<>
+void NRVec<std::complex<float> >::gemv(const std::complex<float> beta,
+                const NRSMat<std::complex<float> > &A, const char trans,
+                const std::complex<float> alpha, const NRVec<std::complex<float> > &x){
+#ifdef DEBUG
+	if(A.ncols() != x.size()) laerror("incompatible dimensions");
+	if(A.nrows() != nn){ laerror("incompatible result vector in gemv"); }
+#endif
+	SAME_LOC3(*this, A, x);
+	copyonwrite();
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		if(tolower(trans) == 't') //row-major lower packed H is column-major upper packed H^T
+			cblas_chpmv(CblasColMajor, CblasUpper, A.ncols(), &alpha, A, x.v, 1, &beta, v, 1);
+		else
+			cblas_chpmv(CblasRowMajor, CblasLower, A.ncols(), &alpha, A, x.v, 1, &beta, v, 1);
+#ifdef CUDALA
+	}else{
+		const cuComplex _alpha = make_cuComplex(alpha.real(), alpha.imag());
+		const cuComplex  _beta = make_cuComplex(beta.real(), beta.imag());
+
+		cublasChpmv('U', A.ncols(), _alpha, (cuComplex*)((const std::complex<float>*)A), (cuComplex*)(x.v), 1, _beta, (cuComplex*)(this->v), 1);
+		TEST_CUBLAS("cublasChpmv");
+	}
+#endif
+}
+
 /***************************************************************************//**
  * computes the outer product of this real vector \f$\vec{a}\f$ with given
  * real vector \f$\vec{b}\f$ and scales the resulting matrix with factor \f$\alpha\f$, i.e.
@@ -716,6 +1106,24 @@ const NRMat<double> NRVec<double>::otimes(const NRVec<double> &b,const bool conj
 	return result;
 }
 
+template<>
+const NRMat<float> NRVec<float>::otimes(const NRVec<float> &b,const bool conj, const float &scale) const {
+
+	SAME_LOC(*this, b);
+	NRMat<float> result(0.0, nn, b.nn, this->getlocation());
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_sger(CblasRowMajor, nn, b.nn, scale, v, 1, b.v, 1, result, b.nn);
+#ifdef CUDALA
+	}else{
+		cublasSger(b.nn, nn, scale, b.v, 1, v, 1, result[0], b.nn);
+		TEST_CUBLAS("cublasSger");
+	}
+#endif
+	return result;
+}
+
 template<>
 const NRVec<double> NRVec<double>::otimes2vec(const NRVec<double> &b,const bool conj, const double &scale) const {
 
@@ -734,6 +1142,24 @@ const NRVec<double> NRVec<double>::otimes2vec(const NRVec<double> &b,const bool
 	return result;
 }
 
+template<>
+const NRVec<float> NRVec<float>::otimes2vec(const NRVec<float> &b,const bool conj, const float &scale) const {
+
+	SAME_LOC(*this, b);
+	NRVec<float> result(0.0, nn*b.nn, this->getlocation());
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_sger(CblasRowMajor, nn, b.nn, scale, v, 1, b.v, 1, result.v, b.nn);
+#ifdef CUDALA
+	}else{
+		cublasSger(b.nn, nn, scale, b.v, 1, v, 1, result.v, b.nn);
+		TEST_CUBLAS("cublasSger");
+	}
+#endif
+	return result;
+}
+
 /***************************************************************************//**
  * computes the outer product of this complex vector \f$\vec{a}\f$ with given
  * complex vector \f$\vec{b}\f$ and scales the resulting matrix with factor \f$\alpha\f$, i.e.
@@ -780,6 +1206,41 @@ NRVec<std::complex<double> >::otimes(const NRVec<std::complex<double> > &b, cons
 	return result;
 }
 
+template<>
+const NRMat<std::complex<float> > 
+NRVec<std::complex<float> >::otimes(const NRVec<std::complex<float> > &b, const bool conj, const std::complex<float> &scale) const {
+	
+	SAME_LOC(*this, b);
+	NRMat<std::complex<float> > result(0., nn, b.nn, this->getlocation());
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		if(conj){
+			cblas_cgerc(CblasRowMajor, nn, b.nn, &scale, v, 1, b.v, 1, result, b.nn);
+		}else{
+			cblas_cgeru(CblasRowMajor, nn, b.nn, &scale, v, 1, b.v, 1, result, b.nn);
+		}
+#ifdef CUDALA
+	}else{
+		if(conj){
+			const cuComplex alpha = make_cuComplex(scale.real(), -scale.imag());
+
+			cublasCgerc(b.nn, nn, alpha, (cuComplex*)(b.v), 1, (cuComplex*)(v), 1, (cuComplex*)(result[0]), b.nn);
+			TEST_CUBLAS("cublasCgerc");
+
+			result.conjugateme();
+		}else{
+			const cuComplex alpha = make_cuComplex(scale.real(), +scale.imag());
+
+			cublasCgeru(b.nn, nn, alpha, (cuComplex*)(b.v), 1, (cuComplex*)(v), 1, (cuComplex*)(result[0]), b.nn);
+			TEST_CUBLAS("cublasCgeru");
+		}
+	}
+#endif
+	return result;
+}
+
 template<>
 const NRVec<std::complex<double> > 
 NRVec<std::complex<double> >::otimes2vec(const NRVec<std::complex<double> > &b, const bool conj, const std::complex<double> &scale) const {
@@ -815,6 +1276,41 @@ NRVec<std::complex<double> >::otimes2vec(const NRVec<std::complex<double> > &b,
 	return result;
 }
 
+template<>
+const NRVec<std::complex<float> > 
+NRVec<std::complex<float> >::otimes2vec(const NRVec<std::complex<float> > &b, const bool conj, const std::complex<float> &scale) const {
+	
+	SAME_LOC(*this, b);
+	NRVec<std::complex<float> > result(0., nn*b.nn, this->getlocation());
+
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		if(conj){
+			cblas_cgerc(CblasRowMajor, nn, b.nn, &scale, v, 1, b.v, 1, result.v, b.nn);
+		}else{
+			cblas_cgeru(CblasRowMajor, nn, b.nn, &scale, v, 1, b.v, 1, result.v, b.nn);
+		}
+#ifdef CUDALA
+	}else{
+		if(conj){
+			const cuComplex alpha = make_cuComplex(scale.real(), -scale.imag());
+
+			cublasCgerc(b.nn, nn, alpha, (cuComplex*)(b.v), 1, (cuComplex*)(v), 1, (cuComplex*)(result.v), b.nn);
+			TEST_CUBLAS("cublasCgerc");
+
+			result.conjugateme();
+		}else{
+			const cuComplex alpha = make_cuComplex(scale.real(), +scale.imag());
+
+			cublasCgeru(b.nn, nn, alpha, (cuComplex*)(b.v), 1, (cuComplex*)(v), 1, (cuComplex*)(result.v), b.nn);
+			TEST_CUBLAS("cublasCgeru");
+		}
+	}
+#endif
+	return result;
+}
+
 
 template<>
 NRVec<std::complex<double> > complexify(const NRVec<double> &rhs) {
@@ -834,6 +1330,24 @@ NRVec<std::complex<double> > complexify(const NRVec<double> &rhs) {
 	return r;
 }
 
+template<>
+NRVec<std::complex<float> > complexify(const NRVec<float> &rhs) {
+	NRVec<std::complex<float> > r(rhs.size(), rhs.getlocation());
+
+#ifdef CUDALA
+	if(rhs.getlocation() == cpu){
+#endif
+		cblas_scopy(rhs.size(), &rhs[0], 1, (float *)(&r[0]), 2);
+#ifdef CUDALA
+	}else{
+		r = 0;
+		cublasScopy(rhs.size(), rhs.v, 1, (float*)(r.v), 2);
+		TEST_CUBLAS("cublasScopy");
+	}
+#endif
+	return r;
+}
+
 
 
 
@@ -1073,6 +1587,22 @@ const double NRVec<double>::sum() const {
         return result;
 }
 
+template <>
+const float NRVec<float>::sum() const {
+        float result=0;
+#ifdef CUDALA
+        if(location == cpu){
+#endif
+                        cblas_saxpy(nn, 1.0, v, 1, &result, 0);
+#ifdef CUDALA
+        }else{
+        laerror("not implemented");
+                }
+        }
+#endif
+        return result;
+}
+
 /***************************************************************************//**
  * sum up the all of the current double-precision complex vector
  * @return sum
@@ -1093,6 +1623,22 @@ const std::complex<double>  NRVec<std::complex<double> >::sum() const {
         return result;
 }
 
+template <>
+const std::complex<float>  NRVec<std::complex<float> >::sum() const {
+        std::complex<float>  result=0;
+#ifdef CUDALA
+        if(location == cpu){
+#endif
+                        cblas_caxpy(nn, &CONEF, v, 1, &result, 0);
+#ifdef CUDALA
+        }else{
+        laerror("not implemented");
+                }
+        }
+#endif
+        return result;
+}
+
 
 
 
@@ -1232,8 +1778,8 @@ INSTANTIZE_CONCAT(std::complex<double>)
 
 
 
-//template class NRVec<float>;
-//template class NRVec<std::complex<float> >;
+template class NRVec<float>;
+template class NRVec<std::complex<float> >;
 template class NRVec<double>;
 template class NRVec<std::complex<double> >;
 template class NRVec<char>;
diff --git a/vec.h b/vec.h
index f783223..ae630b5 100644
--- a/vec.h
+++ b/vec.h
@@ -55,6 +55,7 @@ template <> inline unsigned long long MYABS(const  unsigned long long &x) {retur
  * static constants used in several cblas-routines
  ******************************************************************************/
 const static std::complex<double> CONE = 1.0, CMONE = -1.0, CZERO = 0.0;
+const static std::complex<float> CONEF = 1.0f, CMONEF = -1.0f, CZEROF = 0.0f;
 #ifdef CUDALA
 const static cuDoubleComplex CUONE = {1.,0.}, CUMONE = {-1.,0.}, CUZERO = {0.,0.};
 #endif
@@ -1517,6 +1518,7 @@ NRVec<std::complex<T> > complexify(const NRVec<T> &rhs) {
 	return r;
 }
 template<> NRVec<std::complex<double> > complexify<double>(const NRVec<double> &rhs);
+template<> NRVec<std::complex<float> > complexify<float>(const NRVec<float> &rhs);
 
 /***************************************************************************//**
  * routine for moving vector data between CPU and GPU memory
@@ -1580,6 +1582,24 @@ inline NRVec<double>& NRVec<double>::operator+=(const double &a) {
 	return *this;
 }
 
+template<>
+inline NRVec<float>& NRVec<float>::operator+=(const float &a) {
+	copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_saxpy(nn, 1.0, &a, 0, v, 1);
+#ifdef CUDALA
+	}else{
+		float *d = gpuputdouble(a);
+		cublasSaxpy(nn, 1.0, d, 0, v, 1);
+		TEST_CUBLAS("cublasSaxpy");
+		gpufree(d);
+	}
+#endif
+	return *this;
+}
+
 /***************************************************************************//**
  * adds a complex scalar value \f$\alpha\f$ to all elements of this complex vector \f$\vec{x}\f$
  * \f[\vec{x}_i\leftarrow\vec{x}_i+\alpha\f]
@@ -1604,6 +1624,24 @@ inline NRVec<std::complex<double> >& NRVec<std::complex<double> >::operator+=(co
 	return *this;
 }
 
+template<>
+inline NRVec<std::complex<float> >& NRVec<std::complex<float> >::operator+=(const std::complex<float> &a) {
+	copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_caxpy(nn, &CONEF, &a, 0, v, 1);
+#ifdef CUDALA
+	}else{
+		std::complex<float> *d = gpuputcomplex(a);
+                cublasCaxpy(nn, CUONE, (cuComplex *)d, 0, (cuComplex *)v, 1);
+		TEST_CUBLAS("cublasCaxpy");
+                gpufree(d);
+        }
+#endif
+	return *this;
+}
+
 /***************************************************************************//**
  * subtracts a real scalar value \f$\alpha\f$ from all elements of this real vector \f$\vec{x}\f$
  * \f[\vec{x}_i\leftarrow\vec{x}_i-\alpha\f]
@@ -1628,6 +1666,24 @@ inline NRVec<double>& NRVec<double>::operator-=(const double &a) {
 	return *this;
 }
 
+template<>
+inline NRVec<float>& NRVec<float>::operator-=(const float &a) {
+	copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_saxpy(nn, -1.0, &a, 0, v, 1);
+#ifdef CUDALA
+	}else{
+		float *d = gpuputdouble(a);
+		cublasSaxpy(nn, -1.0, d, 0, v, 1);
+		TEST_CUBLAS("cublasSaxpy");
+		gpufree(d);
+	}
+#endif
+	return *this;
+}
+
 /***************************************************************************//**
  * subtracts a complex scalar value \f$\alpha\f$ from all elements of this complex vector \f$\vec{x}\f$
  * \f[\vec{x}_i\leftarrow\vec{x}_i-\alpha\f]
@@ -1652,6 +1708,24 @@ inline NRVec<std::complex<double> >& NRVec<std::complex<double> >::operator-=(co
 	return *this;
 }
 
+template<>
+inline NRVec<std::complex<float> >& NRVec<std::complex<float> >::operator-=(const std::complex<float> &a) {
+	copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_caxpy(nn, &CMONEF, &a, 0, v, 1);
+#ifdef CUDALA
+	}else{
+		std::complex<float> *d = gpuputcomplex(a);
+		cublasCaxpy(nn, CUMONE, (cuComplex *)d, 0, (cuComplex *)v, 1);
+		TEST_CUBLAS("cublasCaxpy");
+		gpufree(d);
+	}
+#endif
+	return *this;
+}
+
 /***************************************************************************//**
  * adds a real vector \f$\vec{y}\f$ to this real vector \f$\vec{x}\f$
  * \f[\vec{x}\leftarrow\vec{x}+\vec{y}\f]
@@ -1678,6 +1752,26 @@ inline NRVec<double>& NRVec<double>::operator+=(const NRVec<double> &rhs) {
         return *this;
 }
 
+template<>
+inline NRVec<float>& NRVec<float>::operator+=(const NRVec<float> &rhs) {
+#ifdef DEBUG
+        if (nn != rhs.nn) laerror("incompatible dimensions");
+#endif
+	SAME_LOC(*this, rhs);
+        copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_saxpy(nn, 1.0, rhs.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasSaxpy(nn, 1.0, rhs.v, 1, v, 1);
+		TEST_CUBLAS("cubasDaxpy");
+	}
+#endif
+        return *this;
+}
+
 /***************************************************************************//**
  * adds a complex vector \f$\vec{y}\f$ to this complex vector \f$\vec{x}\f$
  * \f[\vec{x}\leftarrow\vec{x}+\vec{y}\f]
@@ -1704,6 +1798,26 @@ inline NRVec<std::complex<double> >& NRVec<std::complex<double> >::operator+=(co
 	return *this;
 }
 
+template<>
+inline NRVec<std::complex<float> >& NRVec<std::complex<float> >::operator+=(const NRVec<std::complex<float> > &rhs) {
+#ifdef DEBUG
+        if (nn != rhs.nn) laerror("incompatible dimensions");
+#endif
+	SAME_LOC(*this, rhs);
+        copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_caxpy(nn, &CONEF, rhs.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasCaxpy(nn, CUONE, (cuComplex*)rhs.v, 1, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCaxpy");
+	}
+#endif
+	return *this;
+}
+
 /***************************************************************************//**
  * subtracts a real vector \f$\vec{y}\f$ from this real vector \f$\vec{x}\f$
  * \f[\vec{x}\leftarrow\vec{x}-\vec{y}\f]
@@ -1730,6 +1844,26 @@ inline NRVec<double> & NRVec<double>::operator-=(const NRVec<double> &rhs) {
         return *this;
 }
 
+template<>
+inline NRVec<float> & NRVec<float>::operator-=(const NRVec<float> &rhs) {
+#ifdef DEBUG
+        if (nn != rhs.nn) laerror("incompatible dimensions");
+#endif
+	SAME_LOC(*this,rhs);
+	copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_saxpy(nn, -1.0, rhs.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasSaxpy(nn, -1.0, rhs.v, 1, v, 1);
+		TEST_CUBLAS("cubasDaxpy");
+	}
+#endif
+        return *this;
+}
+
 /***************************************************************************//**
  * subtracts a complex vector \f$\vec{y}\f$ from this complex vector \f$\vec{x}\f$
  * \f[\vec{x}\leftarrow\vec{x}-\vec{y}\f]
@@ -1756,6 +1890,26 @@ inline NRVec<std::complex<double> >& NRVec<std::complex<double> >::operator-=(co
         return *this;
 }
 
+template<>
+inline NRVec<std::complex<float> >& NRVec<std::complex<float> >::operator-=(const NRVec<std::complex<float> > &rhs) {
+#ifdef DEBUG
+        if (nn != rhs.nn) laerror("incompatible dimensions");
+#endif
+	SAME_LOC(*this, rhs);
+	copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_caxpy(nn, &CMONEF, rhs.v, 1, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasCaxpy(nn, CUMONE, (cuComplex*)rhs.v, 1, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCaxpy");
+	}
+#endif
+        return *this;
+}
+
 /***************************************************************************//**
  * multiplies this real vector \f$\vec{x}\f$ by a real scalar value \f$\alpha\f$
  * \f[\vec{x}_i\leftarrow\alpha\vec{x}_i\f]
@@ -1778,9 +1932,28 @@ inline NRVec<double>& NRVec<double>::operator*=(const double &a) {
         return *this;
 }
 
+template<>
+inline NRVec<float>& NRVec<float>::operator*=(const float &a) {
+        copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_sscal(nn, a, v, 1);
+#ifdef CUDALA
+	}else{
+		cublasSscal(nn, a, v, 1);
+		TEST_CUBLAS("cublasSscal");
+	}
+#endif
+        return *this;
+}
+
 template<>
 inline NRVec<double>& NRVec<double>::operator/=(const double &a) {return *this *= (1./a);}
 
+template<>
+inline NRVec<float>& NRVec<float>::operator/=(const float &a) {return *this *= (1.f/a);}
+
 /***************************************************************************//**
  * multiplies this complex vector \f$\vec{x}\f$ by a complex scalar value \f$\alpha\f$
  * \f[\vec{x}_i\leftarrow\alpha\vec{x}_i\f]
@@ -1804,9 +1977,29 @@ inline NRVec<std::complex<double> >& NRVec<std::complex<double> >::operator*=(co
         return *this;
 }
 
+template<>
+inline NRVec<std::complex<float> >& NRVec<std::complex<float> >::operator*=(const std::complex<float> &a) {
+        copyonwrite();
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_cscal(nn, &a, v, 1);
+#ifdef CUDALA
+	}else{
+		const cuComplex alpha = make_cuComplex(a.real(), a.imag());
+		cublasCscal(nn, alpha, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasCscal");
+	}
+#endif
+        return *this;
+}
+
 template<>
 inline NRVec<std::complex<double> >& NRVec<std::complex<double> >::operator/=(const std::complex<double> &a) {return *this *= (1./a);}
 
+template<>
+inline NRVec<std::complex<float> >& NRVec<std::complex<float> >::operator/=(const std::complex<float> &a) {return *this *= (1.f/a);}
+
 /***************************************************************************//**
  * computes the inner product of this real vector \f$\vec{x}\f$ with given real vector \f$\vec{y]\f$
  * @param[in] rhs real vector \f$\vec{y}\f$
@@ -1832,6 +2025,26 @@ inline const double NRVec<double>::operator*(const NRVec<double> &rhs) const {
 	return ret;
 }
 
+template<>
+inline const float NRVec<float>::operator*(const NRVec<float> &rhs) const {
+	float ret(0.0);
+#ifdef DEBUG
+        if(nn != rhs.nn) laerror("incompatible dimensions");
+#endif
+	SAME_LOC(*this, rhs);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		ret = cblas_sdot(nn, v, 1, rhs.v, 1);
+#ifdef CUDALA
+	}else{
+		ret = cublasSdot(nn, v, 1, rhs.v, 1);
+		TEST_CUBLAS("cublasSdot");
+	}
+#endif
+	return ret;
+}
+
 /***************************************************************************//**
  * computes the inner product of this complex vector \f$\vec{x}\f$ with given complex vector \f$\vec{y}\f$
  * taking conjugation of vector \f$\vec{x}\f$ into account
@@ -1859,6 +2072,27 @@ inline const std::complex<double> NRVec<std::complex<double> >::operator*(const
         return dot;
 }
 
+template<>
+inline const std::complex<float> NRVec<std::complex<float> >::operator*(const NRVec< std::complex<float> > &rhs) const {
+#ifdef DEBUG
+        if(nn != rhs.nn) laerror("incompatible dimensions");
+#endif
+        std::complex<float> dot;
+	SAME_LOC(*this, rhs);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		cblas_cdotc_sub(nn, v, 1, rhs.v, 1, &dot);
+#ifdef CUDALA
+	}else{
+		const cuComplex val = cublasCdotc(nn, (cuComplex*)v, 1, (cuComplex*)rhs.v, 1);
+		TEST_CUBLAS("cublasCdotc");
+		dot = std::complex<float>(cuCreal(val), cuCimag(val)); 
+	}
+#endif
+        return dot;
+}
+
 /***************************************************************************//**
  * computes the inner product of this real vector \f$\vec{x}\f$ with given real data
  * @param[in] y pointer to the double-precision real array (sufficient length assumed)
@@ -1871,6 +2105,12 @@ inline const double NRVec<double>::dot(const double *y, const int stride, bool c
 	return cblas_ddot(nn,v, 1,y,stride);
 }
 
+template<>
+inline const float NRVec<float>::dot(const float *y, const int stride, bool conjugate) const {
+	NOT_GPU(*this);
+	return cblas_sdot(nn,v, 1,y,stride);
+}
+
 /***************************************************************************//**
  * computes the inner product of this complex vector \f$\vec{x}\f$ with given complex data
  * @param[in] y pointer to the double-precision complex array (sufficient length assumed)
@@ -1886,6 +2126,15 @@ inline const std::complex<double> NRVec<std::complex<double> >::dot(const std::c
         return dot;
 }
 
+template<>
+inline const std::complex<float> NRVec<std::complex<float> >::dot(const std::complex<float> *y, const int stride, bool conjugate) const {
+        std::complex<float> dot;
+	NOT_GPU(*this);
+        if(conjugate) cblas_cdotc_sub(nn, v,1,y, stride, &dot);
+	else cblas_cdotu_sub(nn, v,1,y, stride, &dot);
+        return dot;
+}
+
 /***************************************************************************//**
  * computes the sum of the absolute values of the elements of this real vector \f$\vec{x}\f$
  * @return \f$\sum_{i=1}^N\left|\vec{x}_i\right|\f$
@@ -1906,6 +2155,22 @@ inline const double NRVec<double>::asum() const {
 	return ret;
 }
 
+template<>
+inline const float NRVec<float>::asum() const {
+	float ret(0.0);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		ret = cblas_sasum(nn, v, 1);
+#ifdef CUDALA
+	}else{
+		ret = cublasSasum(nn, v, 1);
+		TEST_CUBLAS("cublasSasum");
+	}
+#endif
+	return ret;
+}
+
 
 /***************************************************************************//**
  * for this complex vector \f$\vec{x}\f$ compute the expression
@@ -1928,6 +2193,22 @@ inline const double NRVec<std::complex<double> >::asum() const {
 	return ret;
 }
 
+template<>
+inline const float NRVec<std::complex<float> >::asum() const {
+	float ret(0.0);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		ret = cblas_scasum(nn, v, 1);
+#ifdef CUDALA
+	}else{
+		ret = cublasScasum(nn, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasScasum");
+	}
+#endif
+	return ret;
+}
+
 /***************************************************************************//**
  * for this real vector \f$\vec{x}\f$ (of \f$N\f$ elements) determine the Frobenius norm
  * @return \f$\sum_{i=1}^N\left|\vec{x}_i\right|^2\f$
@@ -1948,6 +2229,22 @@ inline const double  NRVec<double>::norm(int length, int offset, int stride) con
 	return ret;
 }
 
+template<>
+inline const float  NRVec<float>::norm(int length, int offset, int stride) const {
+	float ret(0.);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		ret = cblas_snrm2((length>=0?length:nn), v+offset, stride);
+#ifdef CUDALA
+	}else{
+		ret = cublasSnrm2((length>=0?length:nn), v+offset, stride);
+		TEST_CUBLAS("cublasSnrm2");
+	}
+#endif
+	return ret;
+}
+
 /***************************************************************************//**
  * for this complex vector \f$\vec{x}\f$ (of \f$N\f$ elements) determine the Frobenius norm
  * @return \f$\sum_{i=1}^N\left|\vec{x}_i\right|^2\f$
@@ -1968,6 +2265,22 @@ inline const double NRVec< std::complex<double> >::norm(int length, int offset,
 	return ret;
 }
 
+template<>
+inline const float NRVec< std::complex<float> >::norm(int length, int offset, int stride) const {
+	float ret(0.);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		ret = cblas_scnrm2((length>=0?length:nn), v+offset, stride);
+#ifdef CUDALA
+	}else{
+		ret = cublasScnrm2((length>=0?length:nn), ((cuComplex*)v)+offset, stride);
+		TEST_CUBLAS("cublasSzrm2");
+	}
+#endif
+	return ret;
+}
+
 /***************************************************************************//**
  * for this real vector \f$\vec{x}\f$ determine the element with largest absolute value
  * @return \f$\vec{x}_i\f$ where \f$\left|\vec{x]_i\right|=\mathrm{max}_{j}\left|\vec{x}_{j}\right|\f$
@@ -1991,6 +2304,25 @@ inline const double NRVec<double>::amax() const {
 	return ret;
 }
 
+template<>
+inline const float NRVec<float>::amax() const {
+	if(nn==0) laerror("amax of an empty vector");
+	float ret(0.0);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		ret = v[cblas_isamax(nn, v, 1)];
+#ifdef CUDALA
+	}else{
+		const int pozice = cublasIsamax(nn, v, 1) ;
+		TEST_CUBLAS("cublasIsamax");
+
+		gpuget(1, sizeof(float), v + pozice, &ret);
+	}
+#endif
+	return ret;
+}
+
 template<>
 inline int NRVec<double>::iamax() const {
 	if(nn==0) laerror("iamax of an empty vector");
@@ -2008,6 +2340,23 @@ inline int NRVec<double>::iamax() const {
         return ret; 
 }       
 
+template<>
+inline int NRVec<float>::iamax() const {
+	if(nn==0) laerror("iamax of an empty vector");
+        int ret;
+#ifdef CUDALA
+        if(location == cpu){
+#endif  
+                ret = cblas_isamax(nn, v, 1);
+#ifdef CUDALA
+        }else{
+                ret = cublasIsamax(nn, v, 1) ;
+                TEST_CUBLAS("cublasIsamax");
+        }
+#endif  
+        return ret; 
+}       
+
 
 /***************************************************************************//**
  * for this real vector \f$\vec{x}\f$ determine the element with smallest absolute value
@@ -2039,6 +2388,32 @@ inline const double NRVec<double>::amin() const {
 	return ret;
 }
 
+template<>
+inline const float NRVec<float>::amin() const {
+	if(nn==0) laerror("amin of an empty vector");
+	float ret(0.0);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		//BLAS routine idamin seems no to be supported
+		float val(0.0);
+		int index(0);
+		ret = std::abs(v[0]);
+		for(register int i = 1; i < nn; 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");
+		gpuget(1, sizeof(float), v + pozice, &ret);
+	}
+#endif
+	return ret;
+}
+
 /***************************************************************************//**
  * for a given complex vector \f$\vec{v}\f$, determine the smallest index of the maximum 
  * magnitude element, i.e. maximal element in the 1-norm
@@ -2062,6 +2437,24 @@ inline const std::complex<double> NRVec<std::complex<double> >::amax() const {
 	return ret;
 }
 
+template<>
+inline const std::complex<float> NRVec<std::complex<float> >::amax() const {
+	if(nn==0) laerror("amax of an empty vector");
+	std::complex<float> ret(0., 0.);
+#ifdef CUDALA
+	if(location == cpu){
+#endif
+		ret = v[cblas_icamax(nn, v, 1) ];
+#ifdef CUDALA
+	}else{
+		const int pozice = cublasIcamax(nn, (cuComplex*)v, 1);
+		TEST_CUBLAS("cublasIcamax");
+		gpuget(1, sizeof(std::complex<float>), v + pozice, &ret);
+	}
+#endif
+	return ret;
+}
+
 template<>
 inline int NRVec<std::complex<double> >::iamax() const {
 	if(nn==0) laerror("iamax of an empty vector");
@@ -2079,6 +2472,23 @@ inline int NRVec<std::complex<double> >::iamax() const {
         return ret;
 }
 
+template<>
+inline int NRVec<std::complex<float> >::iamax() const {
+	if(nn==0) laerror("iamax of an empty vector");
+	int ret;
+#ifdef CUDALA
+        if(location == cpu){
+#endif
+                ret = cblas_icamax(nn, v, 1);
+#ifdef CUDALA
+        }else{
+                ret = cublasIcamax(nn, (cuComplex*)v, 1);
+                TEST_CUBLAS("cublasIcamax");
+        }
+#endif
+        return ret;
+}
+
 
 /***************************************************************************//**
  * for a given complex vector \f$\vec{v}\f$, determine the smallest index of the minimum 
@@ -2113,6 +2523,34 @@ inline const std::complex<double> NRVec<std::complex<double> >::amin() const {
 	return ret;
 }
 
+template<>
+inline const std::complex<float> NRVec<std::complex<float> >::amin() const {
+	if(nn==0) laerror("amin of an empty vector");
+	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(std::numeric_limits<float>::max());
+		std::complex<float> z_val(0.0, 0.0);
+
+		for(register int i=0; i < nn; 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");
+		gpuget(1, sizeof(std::complex<float>), v + pozice, &ret);
+	}
+#endif
+	return ret;
+}
+
 
 /***************************************************************************//**
  * routine for raw output