#+TITLE: Utility functions
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org

* Headers                                                         :noexport:
  #+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src

  #+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
int main() {
  qmckl_context context;
  context = qmckl_context_create();

  #+end_src

* Matrix operations

** ~qmckl_transpose~

   Transposes a matrix: $B_{ji} = A_{ij}$

   #+NAME: qmckl_transpose_args
   | qmckl_context | context  | in  | Global state                                 |
   | int64_t       | m        | in  | Number of rows of the input matrix           |
   | int64_t       | n        | in  | Number of columns of the input matrix        |
   | double        | A[][lda] | in  | Array containing the $m \times n$ matrix $A$ |
   | int64_t       | lda      | in  | Leading dimension of array ~A~               |
   | double        | B[][ldb] | out | Array containing the $n \times m$ matrix $B$ |
   | int64_t       | ldb      | in  | Leading dimension of array ~B~               |

*** Requirements

    - ~context~ is not ~QMCKL_NULL_CONTEXT~
    - ~m > 0~
    - ~n > 0~
    - ~lda >= m~
    - ~ldb >= n~
    - ~A~ is allocated with at least $m \times n \times 8$ bytes
    - ~B~ is allocated with at least $n \times m \times 8$ bytes

*** C header

    #+CALL: generate_c_header(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_transpose (
          const qmckl_context context,
          const int64_t m,
          const int64_t n,
          const double* A,
          const int64_t lda,
          double* const B,
          const int64_t ldb ); 
    #+end_src

*** Source
    #+begin_src f90 :tangle (eval f)
integer function qmckl_transpose_f(context, m, n, A, LDA, B, LDB) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  integer*8  , intent(in)  :: m, n
  integer*8  , intent(in)  :: lda
  real*8     , intent(in)  :: A(lda,*)
  integer*8  , intent(in)  :: ldb
  real*8     , intent(out) :: B(ldb,*)

  integer*8 :: i,j

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (m <= 0_8) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  if (n <= 0_8) then
     info = QMCKL_INVALID_ARG_3
     return
  endif

  if (LDA < m) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  if (LDB < n) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  do j=1,m
     do i=1,n
       B(i,j) = A(j,i)
    end do
 end do

end function qmckl_transpose_f
    #+end_src

*** C interface                                                    :noexport:

    #+CALL: generate_c_interface(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_transpose &
        (context, m, n, A, lda, B, ldb) &
        bind(C) result(info)

      use, intrinsic :: iso_c_binding
      implicit none

      integer (c_int64_t) , intent(in)  , value :: context
      integer (c_int64_t) , intent(in)  , value :: m
      integer (c_int64_t) , intent(in)  , value :: n
      real    (c_double ) , intent(in)          :: A(lda,*)
      integer (c_int64_t) , intent(in)  , value :: lda
      real    (c_double ) , intent(out)         :: B(ldb,*)
      integer (c_int64_t) , intent(in)  , value :: ldb

      integer(c_int32_t), external :: qmckl_transpose_f
      info = qmckl_transpose_f &
             (context, m, n, A, lda, B, ldb)

    end function qmckl_transpose
    #+end_src

    #+CALL: generate_f_interface(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")

    #+RESULTS:
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    interface
      integer(c_int32_t) function qmckl_transpose &
          (context, m, n, A, lda, B, ldb) &
          bind(C)
        use, intrinsic :: iso_c_binding
        import
        implicit none

        integer (c_int64_t) , intent(in)  , value :: context
        integer (c_int64_t) , intent(in)  , value :: m
        integer (c_int64_t) , intent(in)  , value :: n
        real    (c_double ) , intent(in)          :: A(lda,*)
        integer (c_int64_t) , intent(in)  , value :: lda
        real    (c_double ) , intent(out)         :: B(ldb,*)
        integer (c_int64_t) , intent(in)  , value :: ldb

      end function qmckl_transpose
    end interface
    #+end_src

*** Test                                                           :noexport:
    #+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_transpose(context) bind(C)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context

  double precision, allocatable :: A(:,:), B(:,:)
  integer*8                     :: m, n, LDA, LDB
  integer*8                     :: i,j
  double precision              :: x

  m = 5
  n = 6
  LDA = m+3
  LDB = n+1

  allocate( A(LDA,n), B(LDB,m) )

  A = 0.d0
  B = 0.d0
  do j=1,n
     do i=1,m
        A(i,j) = -10.d0 + dble(i+j)
     end do
  end do

  test_qmckl_transpose = qmckl_transpose(context, m, n, A, LDA, B, LDB)

  if (test_qmckl_transpose /= QMCKL_SUCCESS) return

  test_qmckl_transpose = QMCKL_FAILURE

  x = 0.d0
  do j=1,n
     do i=1,m
        x = x + (A(i,j)-B(j,i))**2
     end do
  end do
  if (dabs(x) <= 1.d-15) then
     test_qmckl_transpose = QMCKL_SUCCESS
  endif

  deallocate(A,B)
end function test_qmckl_transpose
    #+end_src

    #+begin_src c :comments link :tangle (eval c_test)
qmckl_exit_code test_qmckl_transpose(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_transpose(context));
    #+end_src

* End of files                                                     :noexport:

   #+begin_src c :comments link :tangle (eval c_test)
  assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
  return 0;
}

   #+end_src


# -*- mode: org -*-
# vim: syntax=c