diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index b47830f41f661929f4a732bdeb8df93e6fdbbf93..342183ce78ad2b83a3529fa8039a3dd55883a133 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -177,6 +177,28 @@ void fullMatrix<double>::multWithATranspose(const fullVector<double> &x, const i
 
 }
 
+template<>
+void fullMatrix<double>::gemmWithAtranspose(const fullMatrix<double> &a, const fullMatrix<double> &b,
+                                             double alpha, double beta)
+{
+  int M = size2(), N = size2(), K = a.size1();
+  int LDA = a.size1(), LDB = b.size1(), LDC = size1();
+  F77NAME(dgemm)("T", "N", &M, &N, &K, &alpha, a._data, &LDA, b._data, &LDB,
+                 &beta, _data, &LDC);
+}
+
+template<>
+void fullMatrix<std::complex<double> >::gemmWithAtranspose(const fullMatrix<std::complex<double> > &a,
+                                                             const fullMatrix<std::complex<double> > &b,
+                                                             std::complex<double> alpha,
+                                                             std::complex<double> beta)
+{
+  int M = size2(), N = size2(), K = a.size1();
+  int LDA = a.size1(), LDB = b.size1(), LDC = size1();
+  F77NAME(zgemm)("T", "N", &M, &N, &K, &alpha, a._data, &LDA, b._data, &LDB,
+                 &beta, _data, &LDC);
+}
+
 #endif
 
 
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index cc5fa061736e5789caf3b1b8017f6a92321e872d..85a250c70e1f5a81c317175694e25c742a493345 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -552,7 +552,15 @@ class fullMatrix
      _data[cind+i] = x(i);
   }
 
+  void gemmWithAtranspose(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b,
+            scalar alpha=1., scalar beta=1.)
+#if !defined(HAVE_BLAS)
+  {
+    Msg::Error("gemmWithAtranspose is only available with blas. If blas is not installed please transpose a before used gemm_naive");
+  }
+#endif
+  ;
+
   static void registerBindings(binding *b);
 };
-
 #endif