diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index 556215de1c5bac47746ac4da5141107173232213..6261e39e41c0fae46ccd2bf497b7bb66ab493f13 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -111,6 +111,13 @@ void fullMatrix<std::complex<double> >::gemm(const fullMatrix<std::complex<doubl
                  &beta, _data, &LDC);
 }
 
+template<> 
+void fullMatrix<double>::axpy(const fullMatrix<double> &x,double alpha)
+{
+  int M = _r * _c, INCX = 1, INCY = 1;
+  F77NAME(daxpy)(&M, &alpha, x._data,&INCX, _data, &INCY);
+}
+
 template<> 
 void fullMatrix<double>::mult(const fullVector<double> &x, fullVector<double> &y) const
 {
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 8513c8e22f729755e253ad9f8fac77ab46269463..ea0caf004031d9261df498a6829c20a288b0c82f 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -329,6 +329,14 @@ class fullMatrix
   {
     for (int i = 0; i < _r * _c; i++) _data[i] *= a._data[i];
   }
+  void axpy(const fullMatrix<scalar> &x, scalar alpha=1.)
+#if !defined(HAVE_BLAS)
+  {
+    int n = _r * _c;
+    for (int i = 0; i < n; i++) _data[i] += alpha * x._data[i];
+  }
+#endif
+  ;
   void gemm_naive(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b, 
                   scalar alpha=1., scalar beta=1.)
   {