UP | HOME

BLAS functions

Table of Contents

1 Data types

1.1 Vector

Variable Type Description
size int64_t Dimension of the vector
data double* Elements 
typedef struct qmckl_vector {
  int64_t size;
  double* data;
} qmckl_vector;

qmckl_vector
qmckl_vector_alloc( qmckl_context context, 
                    const int64_t size);

Allocates a new vector. If the allocation failed the size is zero. 

qmckl_vector
qmckl_vector_alloc( qmckl_context context, 
                    const int64_t size)
{
  /* Should always be true by contruction */
  assert (size > (int64_t) 0);

  qmckl_vector result;
  result.size = size;

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = size * sizeof(double);
  result.data = (double*) qmckl_malloc (context, mem_info);

  if (result.data == NULL) {
    result.size = (int64_t) 0;
  }

  return result;
}

qmckl_exit_code
qmckl_vector_free( qmckl_context context, 
                   qmckl_vector* vector);

qmckl_exit_code
qmckl_vector_free( qmckl_context context, 
                   qmckl_vector* vector)
{
  /* Always true */
  assert (vector->data != NULL);

  qmckl_exit_code rc;

  rc = qmckl_free(context, vector->data);
  if (rc != QMCKL_SUCCESS) {
    return rc;
  }

  vector->size = (int64_t) 0;
  vector->data = NULL;
  return QMCKL_SUCCESS;
}

1.2 Matrix

Variable Type Description
size int64_t[2] Dimension of each component
data double* Elements

The dimensions use Fortran ordering: two elements differing by one in the first dimension are consecutive in memory. 

typedef struct qmckl_matrix {
  int64_t size[2];
  double* data;
} qmckl_matrix;

qmckl_matrix
qmckl_matrix_alloc( qmckl_context context, 
                    const int64_t size1,
                    const int64_t size2);

Allocates a new matrix. If the allocation failed the sizes are zero. 

qmckl_matrix
qmckl_matrix_alloc( qmckl_context context, 
                    const int64_t size1,
                    const int64_t size2)
{
  /* Should always be true by contruction */
  assert (size1 * size2 > (int64_t) 0);

  qmckl_matrix result;

  result.size[0] = size1;
  result.size[1] = size2;

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = size1 * size2 * sizeof(double);
  result.data = (double*) qmckl_malloc (context, mem_info);

  if (result.data == NULL) {
    result.size[0] = (int64_t) 0;
    result.size[1] = (int64_t) 0;
  }

  return result;
}

qmckl_exit_code
qmckl_matrix_free( qmckl_context context, 
                   qmckl_matrix* matrix);

qmckl_exit_code
qmckl_matrix_free( qmckl_context context, 
                   qmckl_matrix* matrix)
{
  /* Always true */
  assert (matrix->data != NULL);

  qmckl_exit_code rc;

  rc = qmckl_free(context, matrix->data);
  if (rc != QMCKL_SUCCESS) {
    return rc;
  }
  matrix->data = NULL;
  matrix->size[0] = (int64_t) 0;
  matrix->size[1] = (int64_t) 0;

  return QMCKL_SUCCESS;
}

1.3 Tensor

Variable Type Description
order int64_t Order of the tensor
size int64_t[QMCKL_TENSOR_ORDER_MAX] Dimension of each component
data double* Elements

The dimensions use Fortran ordering: two elements differing by one in the first dimension are consecutive in memory. 

#define QMCKL_TENSOR_ORDER_MAX 16

typedef struct qmckl_tensor {
  int64_t order;
  int64_t size[QMCKL_TENSOR_ORDER_MAX];
  double* data;
} qmckl_tensor;

qmckl_tensor
qmckl_tensor_alloc( qmckl_context context, 
                    const int64_t order,
                    const int64_t* size);

Allocates memory for a tensor. If the allocation failed, the size is zero. 

qmckl_tensor
qmckl_tensor_alloc( qmckl_context context, 
                    const int64_t  order,
                    const int64_t* size)
{
  /* Should always be true by contruction */
  assert (order > 0);
  assert (order <= QMCKL_TENSOR_ORDER_MAX);
  assert (size  != NULL);

  qmckl_tensor result;
  result.order = order;

  int64_t prod_size = (int64_t) 1;
  for (int64_t i=0 ; i<order ; ++i) {
    assert (size[i] > (int64_t) 0);
    result.size[i] = size[i];
    prod_size *= size[i];
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = prod_size * sizeof(double);

  result.data = (double*) qmckl_malloc (context, mem_info);

  if (result.data == NULL) {
    memset(&result, 0, sizeof(qmckl_tensor));
  }

  return result;
}

qmckl_exit_code
qmckl_tensor_free( qmckl_context context, 
                   qmckl_tensor* tensor);

qmckl_exit_code
qmckl_tensor_free( qmckl_context context, 
                   qmckl_tensor* tensor)
{
  /* Always true */
  assert (tensor->data != NULL);

  qmckl_exit_code rc;

  rc = qmckl_free(context, tensor->data);
  if (rc != QMCKL_SUCCESS) {
    return rc;
  }

  memset(tensor, 0, sizeof(qmckl_tensor));

  return QMCKL_SUCCESS;
}

1.4 Reshaping

Reshaping occurs in-place and the pointer to the data is copied.

1.4.1 Vector -> Matrix

qmckl_matrix
qmckl_matrix_of_vector(const qmckl_vector vector,
                       const int64_t size1,
                       const int64_t size2);

Reshapes a vector into a matrix. 

qmckl_matrix
qmckl_matrix_of_vector(const qmckl_vector vector,
                       const int64_t size1,
                       const int64_t size2)
{ 
  /* Always true */
  assert (size1 * size2 == vector.size);

  qmckl_matrix result;

  result.size[0] = size1;
  result.size[1] = size2;
  result.data    = vector.data;

  return result;
}

1.4.2 Vector -> Tensor

qmckl_tensor
qmckl_tensor_of_vector(const qmckl_vector vector,
                       const int64_t order,
                       const int64_t* size);

Reshapes a vector into a tensor. 

qmckl_tensor
qmckl_tensor_of_vector(const qmckl_vector vector,
                       const int64_t order,
                       const int64_t* size)
{ 
  qmckl_tensor result;

  int64_t prod_size = 1;
  for (int64_t i=0 ; i<order ; ++i) {
    result.size[i] = size[i];
    prod_size *= size[i];
  }
  assert (prod_size == vector.size);

  result.data = vector.data;

  return result;
}

1.4.3 Matrix -> Vector

qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix);

Reshapes a matrix into a vector. 

qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix)
{ 
  qmckl_vector result;

  result.size = matrix.size[0] * matrix.size[1];
  result.data = matrix.data;

  return result;
}

1.4.4 Matrix -> Tensor

qmckl_tensor
qmckl_tensor_of_matrix(const qmckl_matrix matrix,
                       const int64_t order,
                       const int64_t* size);

Reshapes a matrix into a tensor. 

qmckl_tensor
qmckl_tensor_of_matrix(const qmckl_matrix matrix,
                       const int64_t order,
                       const int64_t* size)
{ 
  qmckl_tensor result;

  int64_t prod_size = 1;
  for (int64_t i=0 ; i<order ; ++i) {
    result.size[i] = size[i];
    prod_size *= size[i];
  }
  assert (prod_size == matrix.size[0] * matrix.size[1]);

  result.data = matrix.data;

  return result;
}

1.4.5 Tensor -> Vector

qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor);

Reshapes a tensor into a vector. 

qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor)
{ 
  int64_t prod_size = (int64_t) tensor.size[0];
  for (int64_t i=1 ; i<tensor.order ; i++) {
    prod_size *= tensor.size[i];
  }

  qmckl_vector result;

  result.size = prod_size;
  result.data = tensor.data;

  return result;
}

1.4.6 Tensor -> Matrix

qmckl_matrix
qmckl_matrix_of_tensor(const qmckl_tensor tensor,
                       const int64_t size1,
                       const int64_t size2);

Reshapes a tensor into a vector. 

qmckl_matrix
qmckl_matrix_of_tensor(const qmckl_tensor tensor,
                       const int64_t size1,
                       const int64_t size2)
{ 
  /* Always true */
  int64_t prod_size = (int64_t) 1;
  for (int64_t i=0 ; i<tensor.order ; i++) {
    prod_size *= tensor.size[i];
  }
  assert (prod_size == size1 * size2);

  qmckl_matrix result;

  result.size[0] = size1;
  result.size[1] = size2;
  result.data = tensor.data;

  return result;
}

1.5 Access macros

#define qmckl_vec(v, i) v.data[i]
#define qmckl_mat(m, i, j) m.data[(i) + (j)*m.size[0]]

#define qmckl_ten3(t, i, j, k) t.data[(i) + m.size[0]*((j) + size[1]*(k))]
#define qmckl_ten4(t, i, j, k, l) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*(l)))]
#define qmckl_ten5(t, i, j, k, l, m) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*((l) + size[3]*(m))))]

1.6 Copy to/from to double*

qmckl_exit_code
qmckl_double_of_vector(const qmckl_context context,
                       const qmckl_vector vector,
                       double* const target,
                       const int64_t size_max);

Converts a vector to a double*. 

qmckl_exit_code
qmckl_double_of_vector(const qmckl_context context,
                       const qmckl_vector vector,
                       double* const target,
                       const int64_t size_max)
{
  /* Always true by construction */
  assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
  assert (vector.size > (int64_t) 0);
  assert (target != NULL);
  assert (size_max > (int64_t) 0);
  assert (size_max >= vector.size);
  for (int64_t i=0 ; i<vector.size ; ++i) {
    target[i] = vector.data[i];
  }
  return QMCKL_SUCCESS;

}

qmckl_exit_code
qmckl_double_of_matrix(const qmckl_context context,
                       const qmckl_matrix matrix,
                       double* const target,
                       const int64_t size_max);

Converts a matrix to a double*. 

qmckl_exit_code
qmckl_double_of_matrix(const qmckl_context context,
                       const qmckl_matrix matrix,
                       double* const target,
                       const int64_t size_max)
{
  qmckl_vector vector = qmckl_vector_of_matrix(matrix);
  return qmckl_double_of_vector(context, vector, target, size_max);
}

qmckl_exit_code
qmckl_double_of_tensor(const qmckl_context context,
                       const qmckl_tensor tensor,
                       double* const target,
                       const int64_t size_max);

Converts a tensor to a double*. 

qmckl_exit_code
qmckl_double_of_tensor(const qmckl_context context,
                       const qmckl_tensor tensor,
                       double* const target,
                       const int64_t size_max)
{
  qmckl_vector vector = qmckl_vector_of_tensor(tensor);
  return qmckl_double_of_vector(context, vector, target, size_max);
}

qmckl_exit_code
qmckl_vector_of_double(const qmckl_context context,
                       const double* target,
                       const int64_t size_max,
                       qmckl_vector* vector);

Converts a double* to a vector. 

qmckl_exit_code
qmckl_vector_of_double(const qmckl_context context,
                       const double* target,
                       const int64_t size_max,
                       qmckl_vector* vector_out)
{
  qmckl_vector vector = *vector_out;
  /* Always true by construction */
  assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);

  if (vector.size == 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_4,
                           "qmckl_double_of_vector",
                           "Vector not allocated");
  }

  if (vector.size != size_max) {
      return qmckl_failwith( context,
                             QMCKL_INVALID_ARG_4,
                             "qmckl_double_of_vector",
                             "Wrong vector size");
  }

  for (int64_t i=0 ; i<vector.size ; ++i) {
    vector.data[i] = target[i];
  }

  *vector_out = vector;
  return QMCKL_SUCCESS;

}

qmckl_exit_code
qmckl_matrix_of_double(const qmckl_context context,
                       const double* target,
                       const int64_t size_max,
                       qmckl_matrix* matrix);

Converts a matrix to a double*. 

qmckl_exit_code
qmckl_matrix_of_double(const qmckl_context context,
                       const double* target,
                       const int64_t size_max,
                       qmckl_matrix* matrix)
{
  qmckl_vector vector = qmckl_vector_of_matrix(*matrix);
  qmckl_exit_code rc =
    qmckl_vector_of_double(context, target, size_max, &vector);
  *matrix = qmckl_matrix_of_vector(vector, matrix->size[0], matrix->size[1]);
  return rc;
}

qmckl_exit_code
qmckl_tensor_of_double(const qmckl_context context,
                       const double* target,
                       const int64_t size_max,
                       qmckl_tensor* tensor);

Converts a matrix to a double*. 

qmckl_exit_code
qmckl_tensor_of_double(const qmckl_context context,
                       const double* target,
                       const int64_t size_max,
                       qmckl_tensor* tensor)
{
  qmckl_vector vector = qmckl_vector_of_tensor(*tensor);
  qmckl_exit_code rc =
    qmckl_vector_of_double(context, target, size_max, &vector);
  *tensor = qmckl_tensor_of_vector(vector, tensor->order, tensor->size);
  return rc;
}

1.7 Tests

{
  int64_t m = 3;
  int64_t n = 4;
  int64_t p = m*n;
  qmckl_vector vec = qmckl_vector_alloc(context, p);

  for (int64_t i=0 ; i<p ; ++i) 
    qmckl_vec(vec, i) = (double) i;

  for (int64_t i=0 ; i<p ; ++i) 
    assert( vec.data[i] == (double) i );

  qmckl_matrix mat = qmckl_matrix_of_vector(vec, m, n);
  assert (mat.size[0] == m);
  assert (mat.size[1] == n);
  assert (mat.data == vec.data);

  for (int64_t j=0 ; j<n ; ++j)
    for (int64_t i=0 ; i<m ; ++i)
      assert ( qmckl_mat(mat, i, j) == qmckl_vec(vec, i+j*m)) ;

  qmckl_vector vec2 = qmckl_vector_of_matrix(mat);
  assert (vec2.size == p);
  assert (vec2.data == vec.data);
  for (int64_t i=0 ; i<p ; ++i) 
      assert ( qmckl_vec(vec2, i) == qmckl_vec(vec, i) ) ;

  qmckl_vector_free(context, &vec);

}

2 Matrix operations

2.1 qmckl_dgemm

Matrix multiplication:

\[ C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj} \]

Variable Type In/Out Description
context qmckl_context in Global state
TransA char in 'T' is transposed
TransB char in 'T' is transposed
m int64_t in Number of rows of the input matrix
n int64_t in Number of columns of the input matrix
k int64_t in Number of columns of the input matrix
alpha double in α
A double[][lda] in Array containing the matrix \(A\)
lda int64_t in Leading dimension of array A
B double[][ldb] in Array containing the matrix \(B\)
ldb int64_t in Leading dimension of array B
beta double in β
C double[][ldc] out Array containing the matrix \(C\)
ldc int64_t in Leading dimension of array C

Requirements:

  • context is not QMCKL_NULL_CONTEXT
  • m > 0
  • n > 0
  • k > 0
  • lda >= m
  • ldb >= n
  • ldc >= n
  • A is allocated with at least \(m \times k \times 8\) bytes
  • B is allocated with at least \(k \times n \times 8\) bytes
  • C is allocated with at least \(m \times n \times 8\) bytes
qmckl_exit_code qmckl_dgemm (
      const qmckl_context context,
      const char TransA,
      const char TransB,
      const int64_t m,
      const int64_t n,
      const int64_t k,
      const double alpha,
      const double* A,
      const int64_t lda,
      const double* B,
      const int64_t ldb,
      const double beta,
      double* const C,
      const int64_t ldc ); 

integer function qmckl_dgemm_f(context, TransA, TransB, &
     m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  character             , intent(in)  :: TransA, TransB
  integer*8             , intent(in)  :: m, n, k
  double precision      , intent(in)  :: alpha, beta
  integer*8             , intent(in)  :: lda
  double precision      , intent(in)  :: A(lda,*)
  integer*8             , intent(in)  :: ldb
  double precision      , intent(in)  :: B(ldb,*)
  integer*8             , intent(in)  :: ldc
  double precision      , intent(out) :: C(ldc,*)

  info = QMCKL_SUCCESS

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

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

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

  if (k <= 0_8) then
     info = QMCKL_INVALID_ARG_6
     return
  endif

  if (LDA <= 0) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (LDB <= 0) then
     info = QMCKL_INVALID_ARG_11
     return
  endif

  if (LDC <= 0) then
     info = QMCKL_INVALID_ARG_14
     return
  endif

  call dgemm(transA, transB, int(m,4), int(n,4), int(k,4), &
       alpha, A, int(LDA,4), B, int(LDB,4), beta, C, int(LDC,4))

end function qmckl_dgemm_f

2.2 qmckl_matmul

Matrix multiplication:

\[ C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj} \]

Variable Type In/Out Description
context qmckl_context in Global state
TransA char in 'T' is transposed
TransB char in 'T' is transposed
alpha double in α
A qmckl_matrix in Matrix \(A\)
B qmckl_matrix in Matrix \(B\)
beta double in β
C qmckl_matrix out Matrix \(C\) 
qmckl_exit_code qmckl_matmul (
      const qmckl_context context,
      const char TransA,
      const char TransB,
      const double alpha,
      const qmckl_matrix A,
      const qmckl_matrix B,
      const double beta,
      qmckl_matrix* const C ); 

qmckl_exit_code
qmckl_matmul (const qmckl_context context,
              const char TransA,
              const char TransB,
              const double alpha,
              const qmckl_matrix A,
              const qmckl_matrix B,
              const double beta,
              qmckl_matrix* const C )
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_INVALID_CONTEXT;
  }

  qmckl_exit_code rc = QMCKL_SUCCESS;

  if (TransA != 'N' && TransA != 'T') {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_matmul",
                           "TransA should be 'N' or 'T'");
  }

  if (TransB != 'N' && TransB != 'T') {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_matmul",
                           "TransB should be 'N' or 'T'");
  }

  if (A.size[0] < 1) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_5,
                           "qmckl_matmul",
                           "Invalid size for A");
  }

  if (B.size[0] < 1) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_6,
                           "qmckl_matmul",
                           "Invalid size for B");
  }

  if (C == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_8,
                           "qmckl_matmul",
                           "Null pointer");
  }

  int t = 0;
  if (TransA == 'T') t +=1;
  if (TransB == 'T') t +=2;
  /*
    | t | TransA | TransB |
    +---+--------+--------+
    | 0 | N      | N      |
    | 1 | T      | N      |
    | 2 | N      | T      |
    | 3 | T      | T      |
  */

  switch (t) {
  case 0:
    if (A.size[1] != B.size[0]) {
      return qmckl_failwith( context,
                             QMCKL_INVALID_ARG_2,
                             "qmckl_matmul",
                             "A and B have incompatible dimensions");
    }
    C->size[0] = A.size[0];
    C->size[1] = B.size[1];
    rc = qmckl_dgemm (context, 'N', 'N',
                      C->size[0], C->size[1], A.size[1],
                      alpha,
                      A.data, A.size[0],
                      B.data, B.size[0],
                      beta,
                      C->data, C->size[0]);
    break;
  case 1:
    if (A.size[0] != B.size[0]) { 
      return qmckl_failwith( context,
                             QMCKL_INVALID_ARG_2,
                             "qmckl_matmul",
                             "A and B have incompatible dimensions");
    }
    C->size[0] = A.size[1];
    C->size[1] = B.size[1];
    rc = qmckl_dgemm (context, 'T', 'N',
                      C->size[0], C->size[1], A.size[0],
                      alpha,
                      A.data, A.size[0],
                      B.data, B.size[0],
                      beta,
                      C->data, C->size[0]);
    break;
  case 2:
    if (A.size[1] != B.size[1]) {
      return qmckl_failwith( context,
                             QMCKL_INVALID_ARG_2,
                             "qmckl_matmul",
                             "A and B have incompatible dimensions");
    }
    C->size[0] = A.size[0];
    C->size[1] = B.size[0];
    rc = qmckl_dgemm (context, 'N', 'T',
                      C->size[0], C->size[1], A.size[1],
                      alpha,
                      A.data, A.size[0],
                      B.data, B.size[0],
                      beta,
                      C->data, C->size[0]);
    break;
  case 3:
    if (A.size[0] != B.size[1]) {
      return qmckl_failwith( context,
                             QMCKL_INVALID_ARG_2,
                             "qmckl_matmul",
                             "A and B have incompatible dimensions");
    }
    C->size[0] = A.size[1];
    C->size[1] = B.size[0];
    rc = qmckl_dgemm (context, 'T', 'T',
                      C->size[0], C->size[1], A.size[0],
                      alpha,
                      A.data, A.size[0],
                      B.data, B.size[0],
                      beta,
                      C->data, C->size[0]);
    break;
  }
  return rc;
}

2.3 qmckl_adjugate

Given a matrix \(\mathbf{A}\), the adjugate matrix \(\text{adj}(\mathbf{A})\) is the transpose of the cofactors matrix of \(\mathbf{A}\).

\[ \mathbf{B} = \text{adj}(\mathbf{A}) = \text{det}(\mathbf{A}) \, \mathbf{A}^{-1} \]

See also: https://en.wikipedia.org/wiki/Adjugate_matrix

Variable Type In/Out Description
context qmckl_context in Global state
n int64_t in Number of rows and columns of the input matrix
A double[][lda] in Array containing the \(n \times n\) matrix \(A\)
lda int64_t in Leading dimension of array A
B double[][ldb] out Adjugate of \(A\)
ldb int64_t in Leading dimension of array B
det_l double inout determinant of \(A\)

Requirements:

  • context is not QMCKL_NULL_CONTEXT
  • n > 0
  • lda >= m
  • A is allocated with at least \(m \times m \times 8\) bytes
  • ldb >= m
  • B is allocated with at least \(m \times m \times 8\) bytes
qmckl_exit_code qmckl_adjugate (
      const qmckl_context context,
      const int64_t n,
      const double* A,
      const int64_t lda,
      double* const B,
      const int64_t ldb,
      double* det_l );

For small matrices (≤ 5× 5), we use specialized routines for performance motivations. For larger sizes, we rely on the LAPACK library. 

integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  double precision, intent(in)          :: A (LDA,*)
  integer*8, intent(in)                 :: LDA
  double precision, intent(out)         :: B (LDB,*)
  integer*8, intent(in)                 :: LDB
  integer*8, intent(in)                 :: na
  double precision, intent(inout)       :: det_l

  integer    :: i,j

  info = QMCKL_SUCCESS

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

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

  if (LDA <= 0_8) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (LDA < na) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  select case (na)
  case (5)
     call adjugate5(A,LDA,B,LDB,na,det_l)
  case (4)
     call adjugate4(A,LDA,B,LDB,na,det_l)
  case (3)
     call adjugate3(A,LDA,B,LDB,na,det_l)
  case (2)
     call adjugate2(A,LDA,B,LDB,na,det_l)
  case (1)
    det_l = a(1,1)
    b(1,1) = 1.d0
  case default
     call adjugate_general(context, na, A, LDA, B, LDB, det_l)
  end select

end function qmckl_adjugate_f

subroutine adjugate_general(context, na, A, LDA, B, LDB, det_l)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in) :: context
  double precision, intent(in)       :: A (LDA,na)
  integer*8, intent(in)              :: LDA
  double precision, intent(out)      :: B (LDB,na)
  integer*8, intent(in)              :: LDB
  integer*8, intent(in)              :: na
  double precision, intent(inout)    :: det_l

  double precision :: work(LDA*max(na,64))
  integer          :: inf
  integer          :: ipiv(LDA)
  integer          :: lwork
  integer(8)       :: i, j

We first copy the array A into array B. 

B(1:na,1:na) = A(1:na,1:na)

Then, we compute the LU factorization \(LU=A\) using the dgetrf routine. 

call dgetrf(na, na, B, LDB, ipiv, inf )

By convention, the determinant of \(L\) is equal to one, so the determinant of \(A\) is equal to the determinant of \(U\), which is simply computed as the product of its diagonal elements. 

det_l = 1.d0
j=0_8
do i=1,na
 j = j+min(abs(ipiv(i)-i),1)
 det_l = det_l*B(i,i)
enddo

As dgetrf returns \(PLU=A\) where \(P\) is a permutation matrix, the sign of the determinant is computed as \(-1^m\) where \(m\) is the number of permutations. 

if (iand(j,1_8) /= 0_8)  then
  det_l = -det_l
endif

Then, the inverse of \(A\) is computed using dgetri: 

lwork = SIZE(work)
call dgetri(na, B, LDB, ipiv, work, lwork, inf )

and the adjugate matrix is computed as the product of the determinant with the inverse: 

  B(:,:) = B(:,:)*det_l

end subroutine adjugate_general

2.4 qmckl_transpose

Transposes a matrix: \(A^\dagger_{ji} = A_{ij}\).

Variable Type In/Out Description
context qmcklcontext in Global state
A qmcklmatrix in Input matrix
At qmcklmatrix out Transposed matrix 
qmckl_exit_code
qmckl_transpose (qmckl_context context,
                 const qmckl_matrix A,
                 qmckl_matrix At ); 

qmckl_exit_code
qmckl_transpose (qmckl_context context,
                 const qmckl_matrix A,
                 qmckl_matrix At )
{
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_INVALID_CONTEXT;
  }

  if (A.size[0] < 1) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_transpose",
                           "Invalid size for A");
  }

  if (At.data == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_transpose",
                           "Output matrix not allocated");
  }

  if (At.size[0] != A.size[1] || At.size[1] != A.size[0]) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_transpose",
                           "Invalid size for At");
  }

  for (int64_t j=0 ; j<At.size[1] ; ++j) 
    for (int64_t i=0 ; i<At.size[0] ; ++i) 
      qmckl_mat(At, i, j) = qmckl_mat(A, j, i);

  return QMCKL_SUCCESS;
}

2.4.1 Test

{
  qmckl_matrix A;
  qmckl_matrix At;
  A  = qmckl_matrix_alloc(context, 2, 3);
  At = qmckl_matrix_alloc(context, 3, 2);
  for (int j=0 ; j<3 ; ++j)
    for (int i=0 ; i<2 ; ++i)
      qmckl_mat(A, i, j) = (double) 10*i+j;

  qmckl_exit_code rc = qmckl_transpose(context, A, At);
  assert(rc == QMCKL_SUCCESS);

  assert(A.size[0] == At.size[1]);
  assert(A.size[1] == At.size[0]);
  for (int j=0 ; j<3 ; ++j)
    for (int i=0 ; i<2 ; ++i)
      assert (qmckl_mat(A, i, j) == qmckl_mat(At, j, i));

  qmckl_matrix_free(context, &A);
  qmckl_matrix_free(context, &At);
}

Author: TREX CoE

Created: 2022-02-17 Thu 14:40

Validate