From 1a713455f840c15672ca28c2622bf208d0016bfc Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 14 Oct 2021 21:53:00 +0200
Subject: [PATCH] Fixed DGEMM

---
 org/qmckl_blas.org | 29 +++++++++++++++++++++--------
 1 file changed, 21 insertions(+), 8 deletions(-)

diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org
index 9cc1df1..33856fb 100644
--- a/org/qmckl_blas.org
+++ b/org/qmckl_blas.org
@@ -23,7 +23,8 @@ int main() {
 
 ** ~qmckl_dgemm~
 
-   Matrix multiply l$C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}$ using Fortran ~matmul~ function.
+   Matrix multiply: $C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}$ using Fortran ~matmul~ function.
+
    TODO: Add description about the external library dependence.
 
    #+NAME: qmckl_dgemm_args
@@ -41,7 +42,7 @@ int main() {
    | double        | beta     | in  | Array containing the $n \times m$ matrix $B$ |
    | double        | C[][ldc] | out | Array containing the $n \times m$ matrix $B$ |
    | int64_t       | ldc      | in  | Leading dimension of array ~B~               |
-   
+
 *** Requirements
 
     - ~context~ is not ~QMCKL_NULL_CONTEXT~
@@ -75,7 +76,7 @@ int main() {
 	  const int64_t ldb,
 	  const double beta,
 	  double* const C,
-	  const int64_t ldc ); 
+	  const int64_t ldc );
     #+END_src
 
 
@@ -161,11 +162,23 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA,
   endif
 
    if (TransA) then
-     C = matmul(AT,B)
+     if (alpha == 1.d0 && beta == 0.d0) then
+       C = matmul(AT,B)
+     else
+       C = beta*C + alpha*matmul(AT,B)
+     endif
    else if (TransB) then
-     C = matmul(A,BT)
+     if (alpha == 1.d0 && beta == 0.d0) then
+       C = matmul(A,BT)
+     else
+       C = beta*C + alpha*matmul(A,BT)
+     endif
    else
-     C = matmul(A,B)
+     if (alpha == 1.d0 && beta == 0.d0) then
+       C = matmul(A,B)
+     else
+       C = beta*C + alpha*matmul(A,B)
+     endif
    endif
 end function qmckl_dgemm_f
     #+end_src
@@ -251,7 +264,7 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
   double precision              :: x, alpha, beta
 
   TransA = .False.
-  TransB = .False. 
+  TransB = .False.
   m = 1_8
   k = 4_8
   n = 6_8
@@ -302,7 +315,7 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
   deallocate(A,B,C,D)
 end function test_qmckl_dgemm
     #+end_src
-    
+
     #+begin_src c :comments link :tangle (eval c_test)
 qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
 assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));