diff --git a/mat.cc b/mat.cc
index 7508baa..fcaedd9 100644
--- a/mat.cc
+++ b/mat.cc
@@ -5414,6 +5414,44 @@ for(int j=0; j<nn; ++j) (*this)(j,i) *= f;
 
 
 
+/***************************************************************************//**
+ * 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<>
+NRMat<std::complex<double> > complexify<double>(const NRMat<double> &rhs) {
+	NRMat<std::complex<double> > r(rhs.nrows(), rhs.ncols(), rhs.getlocation());
+#ifdef CUDALA
+	if(rhs.getlocation() == cpu){
+#endif
+		cblas_dcopy((size_t)rhs.nrows()*rhs.ncols(), &rhs(0,0), 1, (double*)(&r(0,0)), 2);
+#ifdef CUDALA
+	}else{
+		cublasDcopy((size_t)rhs.nrows()*rhs.ncols(), &rhs(0,0), 1, (double*)(&r(0,0)), 2);
+		TEST_CUBLAS("cublasDcopy");//"NRMat<std::complex<double> > complexify(const NRMat<double> &)"
+	}
+#endif
+	return r;
+}
+
+template<>
+NRMat<std::complex<float> > complexify<float>(const NRMat<float> &rhs) {
+	NRMat<std::complex<float> > r(rhs.nrows(), rhs.ncols(), rhs.getlocation());
+#ifdef CUDALA
+	if(rhs.getlocation() == cpu){
+#endif
+		cblas_scopy((size_t)rhs.nrows()*rhs.ncols(), &rhs(0,0), 1, (float*)(&r(0,0)), 2);
+#ifdef CUDALA
+	}else{
+		cublasScopy((size_t)rhs.nrows()*rhs.ncols(), &rhs(0,0), 1, (float*)(&r(0,0)), 2);
+		TEST_CUBLAS("cublasScopy");//"NRMat<std::complex<float> > complexify(const NRMat<float> &)"
+	}
+#endif
+	return r;
+}
+
+
 /***************************************************************************//**
  * forced instantization in the corresponding object file
  ******************************************************************************/
diff --git a/mat.h b/mat.h
index 19cb524..3009fb8 100644
--- a/mat.h
+++ b/mat.h
@@ -1454,6 +1454,10 @@ NRMat<std::complex<T> > complexify(const NRMat<T> &rhs) {
 	return r;
 }
 
+//BLAS versions for double and float are in mat.cc
+template<> NRMat<std::complex<double> > complexify<double>(const NRMat<double> &rhs);
+template<> NRMat<std::complex<float> > complexify<float>(const NRMat<float> &rhs);
+
 //this is general for any type, complexmatrix() uses blas for doubles
 template<typename T>
 NRMat<std::complex<T> > complexify(const NRMat<T> &rhsr, const NRMat<T> &rhsi) {