mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-08-26 23:21:46 +02:00
2750 lines
93 KiB
Org Mode
2750 lines
93 KiB
Org Mode
#+TITLE: BLAS 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 :tangle (eval h_private_type)
|
|
#ifndef QMCKL_BLAS_HPT
|
|
#define QMCKL_BLAS_HPT
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_func)
|
|
#ifndef QMCKL_BLAS_HPF
|
|
#define QMCKL_BLAS_HPF
|
|
#include "qmckl_blas_private_type.h"
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c)
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#ifdef HAVE_STDINT_H
|
|
#include <stdint.h>
|
|
#elif HAVE_INTTYPES_H
|
|
#include <inttypes.h>
|
|
#endif
|
|
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include <stdbool.h>
|
|
#include <assert.h>
|
|
#include <math.h>
|
|
|
|
#include "qmckl.h"
|
|
#include "qmckl_context_private_type.h"
|
|
#include "qmckl_memory_private_type.h"
|
|
#include "qmckl_blas_private_type.h"
|
|
|
|
#include "qmckl_memory_private_func.h"
|
|
#include "qmckl_blas_private_func.h"
|
|
#+end_src
|
|
|
|
#+begin_src c :comments link :tangle (eval c_test) :noweb yes
|
|
#include "qmckl.h"
|
|
#include <stdio.h>
|
|
#include <assert.h>
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#include "qmckl_blas_private_type.h"
|
|
#include "qmckl_blas_private_func.h"
|
|
|
|
int main() {
|
|
qmckl_context context;
|
|
context = qmckl_context_create();
|
|
|
|
#+end_src
|
|
|
|
* -
|
|
:PROPERTIES:
|
|
:UNNUMBERED: t
|
|
:END:
|
|
|
|
Basic linear algebra data types and operations are described in this file.
|
|
The data types are private, so that HPC implementations can use
|
|
whatever data structures they prefer.
|
|
|
|
These data types are expected to be used internally in QMCkl. They
|
|
are not intended to be passed to external codes.
|
|
|
|
* Data types
|
|
|
|
** Vector
|
|
|
|
| Variable | Type | Description |
|
|
|----------+-----------+-------------------------|
|
|
| ~size~ | ~int64_t~ | Dimension of the vector |
|
|
| ~data~ | ~double*~ | Elements |
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_type) :exports none
|
|
typedef struct qmckl_vector {
|
|
double* restrict data;
|
|
int64_t size;
|
|
} qmckl_vector;
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_vector
|
|
qmckl_vector_alloc( qmckl_context context,
|
|
const int64_t size);
|
|
#+end_src
|
|
|
|
Allocates a new vector. If the allocation failed the size is zero.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_exit_code
|
|
qmckl_vector_free( qmckl_context context,
|
|
qmckl_vector* vector);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
qmckl_exit_code
|
|
qmckl_vector_free( qmckl_context context,
|
|
qmckl_vector* vector)
|
|
{
|
|
if (vector == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_vector_free",
|
|
"Null pointer");
|
|
}
|
|
|
|
/* 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;
|
|
}
|
|
#+end_src
|
|
|
|
** 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.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_type) :exports none
|
|
typedef struct qmckl_matrix {
|
|
double* restrict data;
|
|
int64_t size[2];
|
|
} qmckl_matrix;
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_matrix
|
|
qmckl_matrix_alloc( qmckl_context context,
|
|
const int64_t size1,
|
|
const int64_t size2);
|
|
#+end_src
|
|
|
|
Allocates a new matrix. If the allocation failed the sizes are zero.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_exit_code
|
|
qmckl_matrix_free( qmckl_context context,
|
|
qmckl_matrix* matrix);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
qmckl_exit_code
|
|
qmckl_matrix_free( qmckl_context context,
|
|
qmckl_matrix* matrix)
|
|
{
|
|
if (matrix == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_matrix_free",
|
|
"Null pointer");
|
|
}
|
|
|
|
/* 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;
|
|
}
|
|
#+end_src
|
|
|
|
** 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.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_type) :exports none
|
|
#define QMCKL_TENSOR_ORDER_MAX 16
|
|
|
|
typedef struct qmckl_tensor {
|
|
double* restrict data;
|
|
int64_t order;
|
|
int64_t size[QMCKL_TENSOR_ORDER_MAX];
|
|
} qmckl_tensor;
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_tensor
|
|
qmckl_tensor_alloc( qmckl_context context,
|
|
const int64_t order,
|
|
const int64_t* size);
|
|
#+end_src
|
|
|
|
Allocates memory for a tensor. If the allocation failed, the size
|
|
is zero.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_exit_code
|
|
qmckl_tensor_free (qmckl_context context,
|
|
qmckl_tensor* tensor);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
qmckl_exit_code
|
|
qmckl_tensor_free( qmckl_context context,
|
|
qmckl_tensor* tensor)
|
|
{
|
|
if (tensor == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_tensor_free",
|
|
"Null pointer");
|
|
}
|
|
|
|
/* 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;
|
|
}
|
|
#+end_src
|
|
|
|
** Reshaping
|
|
|
|
Reshaping occurs in-place and the pointer to the data is copied.
|
|
|
|
*** Vector -> Matrix
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_matrix
|
|
qmckl_matrix_of_vector(const qmckl_vector vector,
|
|
const int64_t size1,
|
|
const int64_t size2);
|
|
#+end_src
|
|
|
|
Reshapes a vector into a matrix.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
*** Vector -> Tensor
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_tensor
|
|
qmckl_tensor_of_vector(const qmckl_vector vector,
|
|
const int64_t order,
|
|
const int64_t* size);
|
|
#+end_src
|
|
|
|
Reshapes a vector into a tensor.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
*** Matrix -> Vector
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_vector
|
|
qmckl_vector_of_matrix(const qmckl_matrix matrix);
|
|
#+end_src
|
|
|
|
Reshapes a matrix into a vector.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
*** Matrix -> Tensor
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_tensor
|
|
qmckl_tensor_of_matrix(const qmckl_matrix matrix,
|
|
const int64_t order,
|
|
const int64_t* size);
|
|
#+end_src
|
|
|
|
Reshapes a matrix into a tensor.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
*** Tensor -> Vector
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_vector
|
|
qmckl_vector_of_tensor(const qmckl_tensor tensor);
|
|
#+end_src
|
|
|
|
Reshapes a tensor into a vector.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
*** Tensor -> Matrix
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_matrix
|
|
qmckl_matrix_of_tensor(const qmckl_tensor tensor,
|
|
const int64_t size1,
|
|
const int64_t size2);
|
|
#+end_src
|
|
|
|
Reshapes a tensor into a vector.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
** Access macros
|
|
|
|
Macros are provided to ease the access to vectors, matrices and
|
|
tensors. Matrices use column-major ordering, as in Fortran.
|
|
|
|
#+begin_src c :tangle no
|
|
double qmckl_vec (qmckl_vector v, int64_t i); // v[i]
|
|
double qmckl_mat (qmckl_matrix m, int64_t i, int64_t j) // m[j][i]
|
|
|
|
double qmckl_ten3(qmckl_tensor t, int64_t i, int64_t j, int64_t k); // t[k][j][i]
|
|
double qmckl_ten4(qmckl_tensor t, int64_t i, int64_t j, int64_t k, int64_t l);
|
|
double qmckl_ten5(qmckl_tensor t, int64_t i, int64_t j, int64_t k, int64_t l, int64_t m);
|
|
...
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :exports none
|
|
#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) + t.size[0]*((j) + t.size[1]*(k))]
|
|
#define qmckl_ten4(t, i, j, k, l) t.data[(i) + t.size[0]*((j) + t.size[1]*((k) + t.size[2]*(l)))]
|
|
#define qmckl_ten5(t, i, j, k, l, m) t.data[(i) + t.size[0]*((j) + t.size[1]*((k) + t.size[2]*((l) + t.size[3]*(m))))]
|
|
#+end_src
|
|
|
|
For example:
|
|
|
|
** Set all elements
|
|
|
|
*** Vector
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_vector
|
|
qmckl_vector_set(qmckl_vector vector, double value);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
qmckl_vector
|
|
qmckl_vector_set(qmckl_vector vector, double value)
|
|
{
|
|
for (int64_t i=0 ; i<vector.size ; ++i) {
|
|
qmckl_vec(vector, i) = value;
|
|
}
|
|
return vector;
|
|
}
|
|
#+end_src
|
|
|
|
*** Matrix
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_matrix
|
|
qmckl_matrix_set(qmckl_matrix matrix, double value);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
qmckl_matrix
|
|
qmckl_matrix_set(qmckl_matrix matrix, double value)
|
|
{
|
|
qmckl_vector vector = qmckl_vector_of_matrix(matrix);
|
|
for (int64_t i=0 ; i<vector.size ; ++i) {
|
|
qmckl_vec(vector, i) = value;
|
|
}
|
|
return qmckl_matrix_of_vector(vector, matrix.size[0], matrix.size[1]);
|
|
}
|
|
#+end_src
|
|
|
|
*** Tensor
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_tensor
|
|
qmckl_tensor_set(qmckl_tensor tensor, double value);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
qmckl_tensor
|
|
qmckl_tensor_set(qmckl_tensor tensor, double value)
|
|
{
|
|
qmckl_vector vector = qmckl_vector_of_tensor(tensor);
|
|
for (int64_t i=0 ; i<vector.size ; ++i) {
|
|
qmckl_vec(vector, i) = value;
|
|
}
|
|
return qmckl_tensor_of_vector(vector, tensor.order, tensor.size);
|
|
}
|
|
#+end_src
|
|
|
|
** Copy to/from to ~double*~
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_exit_code
|
|
qmckl_double_of_vector(const qmckl_context context,
|
|
const qmckl_vector vector,
|
|
double* const target,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
Converts a vector to a ~double*~.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_exit_code
|
|
qmckl_double_of_matrix(const qmckl_context context,
|
|
const qmckl_matrix matrix,
|
|
double* const target,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
Converts a matrix to a ~double*~.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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);
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_exit_code
|
|
qmckl_double_of_tensor(const qmckl_context context,
|
|
const qmckl_tensor tensor,
|
|
double* const target,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
Converts a tensor to a ~double*~.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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);
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_exit_code
|
|
qmckl_vector_of_double(const qmckl_context context,
|
|
const double* target,
|
|
const int64_t size_max,
|
|
qmckl_vector* vector);
|
|
#+end_src
|
|
|
|
Converts a ~double*~ to a vector.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_exit_code
|
|
qmckl_matrix_of_double(const qmckl_context context,
|
|
const double* target,
|
|
const int64_t size_max,
|
|
qmckl_matrix* matrix);
|
|
#+end_src
|
|
|
|
Converts a ~double*~ to a matrix.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_exit_code
|
|
qmckl_tensor_of_double(const qmckl_context context,
|
|
const double* target,
|
|
const int64_t size_max,
|
|
qmckl_tensor* tensor);
|
|
#+end_src
|
|
|
|
Converts a ~double*~ to a tensor.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
** Tests :noexport:
|
|
|
|
#+begin_src c :comments link :tangle (eval c_test) :exports none
|
|
{
|
|
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);
|
|
|
|
}
|
|
#+end_src
|
|
* Matrix operations
|
|
|
|
** ~qmckl_dgemm~
|
|
|
|
Matrix multiplication with a BLAS interface:
|
|
|
|
\[
|
|
C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
|
|
\]
|
|
|
|
# TODO: Add description about the external library dependence.
|
|
|
|
#+NAME: qmckl_dgemm_args
|
|
| 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 | \alpha |
|
|
| ~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 | \beta |
|
|
| ~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
|
|
|
|
#+CALL: generate_c_header(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
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 );
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval f) :exports none
|
|
integer function qmckl_dgemm_f(context, TransA, TransB, &
|
|
m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) &
|
|
result(info)
|
|
use qmckl
|
|
#ifdef HAVE_LIBQMCKLDGEMM
|
|
use qmckl_dgemm_tiled_module
|
|
#endif
|
|
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,*)
|
|
double precision,allocatable,dimension(:,:) :: A1
|
|
double precision,allocatable,dimension(:,:) :: B1
|
|
double precision,allocatable,dimension(:,:) :: C1
|
|
integer*8 :: i, j, LDA1, LDB1, LDC1
|
|
|
|
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
|
|
|
|
#ifdef HAVE_LIBQMCKLDGEMM
|
|
! Copy A to A1
|
|
allocate(A1(k,m))
|
|
do j=1,m
|
|
do i=1,k
|
|
A1(i,j) = A(j,i)
|
|
end do
|
|
end do
|
|
|
|
! Copy B to B1
|
|
allocate(B1(n,k))
|
|
do j=1,k
|
|
do i=1,n
|
|
B1(i,j) = B(j,i)
|
|
end do
|
|
end do
|
|
|
|
! Copy C to C1
|
|
allocate(C1(n,m))
|
|
do j=1,m
|
|
do i=1,n
|
|
C1(i,j) = C(j,i)
|
|
end do
|
|
end do
|
|
|
|
LDA1 = size(A1,1)
|
|
LDB1 = size(B1,1)
|
|
LDC1 = size(C1,1)
|
|
|
|
info = qmckl_dgemm_tiled_avx2(int(m,8), int(n,8), int(k,8), &
|
|
A1, int(LDA1,8), B1, int(LDB1,8), C1, int(LDC1,8))
|
|
|
|
do j=1,n
|
|
do i=1,m
|
|
C(i,j) = alpha*C1(j,i) + beta*C(i,j)
|
|
end do
|
|
end do
|
|
|
|
deallocate(A1,B1,C1)
|
|
#else
|
|
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))
|
|
#endif
|
|
|
|
|
|
end function qmckl_dgemm_f
|
|
#+end_src
|
|
|
|
*** C interface :noexport:
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_dgemm &
|
|
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
character , intent(in) , value :: TransA
|
|
character , intent(in) , value :: TransB
|
|
integer (c_int64_t) , intent(in) , value :: m
|
|
integer (c_int64_t) , intent(in) , value :: n
|
|
integer (c_int64_t) , intent(in) , value :: k
|
|
real (c_double ) , intent(in) , value :: alpha
|
|
real (c_double ) , intent(in) :: A(lda,*)
|
|
integer (c_int64_t) , intent(in) , value :: lda
|
|
real (c_double ) , intent(in) :: B(ldb,*)
|
|
integer (c_int64_t) , intent(in) , value :: ldb
|
|
real (c_double ) , intent(in) , value :: beta
|
|
real (c_double ) , intent(out) :: C(ldc,*)
|
|
integer (c_int64_t) , intent(in) , value :: ldc
|
|
|
|
integer(c_int32_t), external :: qmckl_dgemm_f
|
|
info = qmckl_dgemm_f &
|
|
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
|
|
|
|
end function qmckl_dgemm
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_dgemm &
|
|
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
character , intent(in) , value :: TransA
|
|
character , intent(in) , value :: TransB
|
|
integer (c_int64_t) , intent(in) , value :: m
|
|
integer (c_int64_t) , intent(in) , value :: n
|
|
integer (c_int64_t) , intent(in) , value :: k
|
|
real (c_double ) , intent(in) , value :: alpha
|
|
real (c_double ) , intent(in) :: A(lda,*)
|
|
integer (c_int64_t) , intent(in) , value :: lda
|
|
real (c_double ) , intent(in) :: B(ldb,*)
|
|
integer (c_int64_t) , intent(in) , value :: ldb
|
|
real (c_double ) , intent(in) , value :: beta
|
|
real (c_double ) , intent(out) :: C(ldc,*)
|
|
integer (c_int64_t) , intent(in) , value :: ldc
|
|
|
|
end function qmckl_dgemm
|
|
end interface
|
|
#+end_src
|
|
|
|
*** Test :noexport:
|
|
|
|
#+begin_src f90 :tangle (eval f_test)
|
|
integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in), value :: context
|
|
|
|
double precision, allocatable :: A(:,:), B(:,:), C(:,:), D(:,:)
|
|
integer*8 :: m, n, k, LDA, LDB, LDC
|
|
integer*8 :: i,j,l
|
|
character :: TransA, TransB
|
|
double precision :: x, alpha, beta
|
|
|
|
TransA = 'N'
|
|
TransB = 'N'
|
|
m = 1_8
|
|
k = 4_8
|
|
n = 6_8
|
|
LDA = m
|
|
LDB = k
|
|
LDC = m
|
|
|
|
allocate( A(LDA,k), B(LDB,n) , C(LDC,n), D(LDC,n))
|
|
|
|
A = 0.d0
|
|
B = 0.d0
|
|
C = 0.d0
|
|
D = 0.d0
|
|
alpha = 1.0d0
|
|
beta = 0.0d0
|
|
do j=1,k
|
|
do i=1,m
|
|
A(i,j) = -10.d0 + dble(i+j)
|
|
end do
|
|
end do
|
|
|
|
do j=1,n
|
|
do i=1,k
|
|
B(i,j) = -10.d0 + dble(i+j)
|
|
end do
|
|
end do
|
|
|
|
test_qmckl_dgemm = qmckl_dgemm(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC)
|
|
|
|
if (test_qmckl_dgemm /= QMCKL_SUCCESS) return
|
|
|
|
test_qmckl_dgemm = QMCKL_FAILURE
|
|
|
|
x = 0.d0
|
|
do j=1,n
|
|
do i=1,m
|
|
do l=1,k
|
|
D(i,j) = D(i,j) + A(i,l)*B(l,j)
|
|
end do
|
|
x = x + (D(i,j) - C(i,j))**2
|
|
end do
|
|
end do
|
|
|
|
if (dabs(x) <= 1.d-12) then
|
|
test_qmckl_dgemm = QMCKL_SUCCESS
|
|
endif
|
|
|
|
deallocate(A,B,C,D)
|
|
|
|
end function test_qmckl_dgemm
|
|
#+end_src
|
|
|
|
#+begin_src c :comments link :tangle (eval c_test) :exports none
|
|
qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
|
|
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
|
|
#+end_src
|
|
|
|
** ~qmckl_dgemm_safe~
|
|
|
|
"Size-safe" proxy function with the same functionality as ~qmckl_dgemm~
|
|
but with 3 additional arguments. These arguments ~size_max_M~ (where ~M~ is a matix)
|
|
are required primarily for the Python API, where compatibility with
|
|
NumPy arrays implies that sizes of the input and output arrays are provided.
|
|
|
|
#+NAME: qmckl_dgemm_safe_args
|
|
| 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 | \alpha |
|
|
| ~A~ | ~double[][lda]~ | in | Array containing the matrix $A$ |
|
|
| ~size_max_A~ | ~int64_t~ | in | Size of the matrix $A$ |
|
|
| ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ |
|
|
| ~B~ | ~double[][ldb]~ | in | Array containing the matrix $B$ |
|
|
| ~size_max_B~ | ~int64_t~ | in | Size of the matrix $B$ |
|
|
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
|
|
| ~beta~ | ~double~ | in | \beta |
|
|
| ~C~ | ~double[][ldc]~ | out | Array containing the matrix $C$ |
|
|
| ~size_max_C~ | ~int64_t~ | in | Size of 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
|
|
- ~size_max_A >= m * k~
|
|
- ~size_max_B >= k * n~
|
|
- ~size_max_C >= m * n~
|
|
|
|
#+CALL: generate_c_header(table=qmckl_dgemm_safe_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm_safe")
|
|
|
|
#+RESULTS:
|
|
#+BEGIN_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_dgemm_safe (
|
|
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 size_max_A,
|
|
const int64_t lda,
|
|
const double* B,
|
|
const int64_t size_max_B,
|
|
const int64_t ldb,
|
|
const double beta,
|
|
double* const C,
|
|
const int64_t size_max_C,
|
|
const int64_t ldc );
|
|
#+END_src
|
|
|
|
#+begin_src f90 :tangle (eval f) :exports none
|
|
integer function qmckl_dgemm_safe_f(context, TransA, TransB, &
|
|
m, n, k, alpha, A, size_A, LDA, B, size_B, LDB, beta, C, size_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
|
|
integer*8 , intent(in) :: size_A
|
|
double precision , intent(in) :: A(lda,*)
|
|
integer*8 , intent(in) :: ldb
|
|
integer*8 , intent(in) :: size_B
|
|
double precision , intent(in) :: B(ldb,*)
|
|
integer*8 , intent(in) :: ldc
|
|
integer*8 , intent(in) :: size_C
|
|
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_10
|
|
return
|
|
endif
|
|
|
|
if (LDB <= 0) then
|
|
info = QMCKL_INVALID_ARG_13
|
|
return
|
|
endif
|
|
|
|
if (LDC <= 0) then
|
|
info = QMCKL_INVALID_ARG_17
|
|
return
|
|
endif
|
|
|
|
if (size_A <= 0) then
|
|
info = QMCKL_INVALID_ARG_9
|
|
return
|
|
endif
|
|
|
|
if (size_B <= 0) then
|
|
info = QMCKL_INVALID_ARG_12
|
|
return
|
|
endif
|
|
|
|
if (size_C <= 0) then
|
|
info = QMCKL_INVALID_ARG_16
|
|
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_safe_f
|
|
#+end_src
|
|
|
|
*** C interface :noexport:
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_dgemm_safe_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm_safe")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_dgemm_safe &
|
|
(context, TransA, TransB, m, n, k, alpha, A, size_A, lda, B, size_B, ldb, beta, C, size_C, ldc) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
character , intent(in) , value :: TransA
|
|
character , intent(in) , value :: TransB
|
|
integer (c_int64_t) , intent(in) , value :: m
|
|
integer (c_int64_t) , intent(in) , value :: n
|
|
integer (c_int64_t) , intent(in) , value :: k
|
|
real (c_double ) , intent(in) , value :: alpha
|
|
real (c_double ) , intent(in) :: A(lda,*)
|
|
integer (c_int64_t) , intent(in) , value :: size_A
|
|
integer (c_int64_t) , intent(in) , value :: lda
|
|
real (c_double ) , intent(in) :: B(ldb,*)
|
|
integer (c_int64_t) , intent(in) , value :: size_B
|
|
integer (c_int64_t) , intent(in) , value :: ldb
|
|
real (c_double ) , intent(in) , value :: beta
|
|
real (c_double ) , intent(out) :: C(ldc,*)
|
|
integer (c_int64_t) , intent(in) , value :: size_C
|
|
integer (c_int64_t) , intent(in) , value :: ldc
|
|
|
|
integer(c_int32_t), external :: qmckl_dgemm_safe_f
|
|
info = qmckl_dgemm_safe_f &
|
|
(context, TransA, TransB, m, n, k, alpha, A, size_A, lda, B, size_B, ldb, beta, C, size_C, ldc)
|
|
|
|
end function qmckl_dgemm_safe
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_dgemm_safe_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm_safe")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_dgemm_safe &
|
|
(context, TransA, TransB, m, n, k, alpha, A, size_A, lda, B, size_B, ldb, beta, C, size_C, ldc) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
character , intent(in) , value :: TransA
|
|
character , intent(in) , value :: TransB
|
|
integer (c_int64_t) , intent(in) , value :: m
|
|
integer (c_int64_t) , intent(in) , value :: n
|
|
integer (c_int64_t) , intent(in) , value :: k
|
|
real (c_double ) , intent(in) , value :: alpha
|
|
real (c_double ) , intent(in) :: A(lda,*)
|
|
integer (c_int64_t) , intent(in) , value :: size_A
|
|
integer (c_int64_t) , intent(in) , value :: lda
|
|
real (c_double ) , intent(in) :: B(ldb,*)
|
|
integer (c_int64_t) , intent(in) , value :: size_B
|
|
integer (c_int64_t) , intent(in) , value :: ldb
|
|
real (c_double ) , intent(in) , value :: beta
|
|
real (c_double ) , intent(out) :: C(ldc,*)
|
|
integer (c_int64_t) , intent(in) , value :: size_C
|
|
integer (c_int64_t) , intent(in) , value :: ldc
|
|
|
|
end function qmckl_dgemm_safe
|
|
end interface
|
|
#+end_src
|
|
|
|
** ~qmckl_matmul~
|
|
|
|
Matrix multiplication using the =qmckl_matrix= data type:
|
|
|
|
\[
|
|
C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
|
|
\]
|
|
|
|
# TODO: Add description about the external library dependence.
|
|
|
|
#+NAME: qmckl_matmul_args
|
|
| 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 | \alpha |
|
|
| ~A~ | ~qmckl_matrix~ | in | Matrix $A$ |
|
|
| ~B~ | ~qmckl_matrix~ | in | Matrix $B$ |
|
|
| ~beta~ | ~double~ | in | \beta |
|
|
| ~C~ | ~qmckl_matrix~ | out | Matrix $C$ |
|
|
|
|
#+CALL: generate_c_header(table=qmckl_matmul_args,rettyp="qmckl_exit_code",fname="qmckl_matmul")
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_private_func) :comments org
|
|
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 );
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c) :comments org :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
*** Test :noexport:
|
|
|
|
#+begin_src python :exports none :results output
|
|
import numpy as np
|
|
|
|
A = np.array([[ 1., 2., 3., 4. ],
|
|
[ 5., 6., 7., 8. ],
|
|
[ 9., 10., 11., 12. ]])
|
|
|
|
B = np.array([[ 1., -2., 3., 4., 5. ],
|
|
[ 5., -6., 7., 8., 9. ],
|
|
[ 9., 10., 11., 12., 13. ],
|
|
[10., 11., 12., 15., 14. ]])
|
|
|
|
C = 0.5 * A @ B
|
|
print(A.T)
|
|
print(B.T)
|
|
print(C.T)
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
#+begin_example
|
|
[[ 1. 5. 9.]
|
|
[ 2. 6. 10.]
|
|
[ 3. 7. 11.]
|
|
[ 4. 8. 12.]]
|
|
[[ 1. 5. 9. 10.]
|
|
[-2. -6. 10. 11.]
|
|
[ 3. 7. 11. 12.]
|
|
[ 4. 8. 12. 15.]
|
|
[ 5. 9. 13. 14.]]
|
|
[[ 39. 89. 139.]
|
|
[ 30. 56. 82.]
|
|
[ 49. 115. 181.]
|
|
[ 58. 136. 214.]
|
|
[ 59. 141. 223.]]
|
|
#+end_example
|
|
|
|
#+begin_src c :comments link :tangle (eval c_test) :exports none
|
|
{
|
|
double a[12] = { 1., 5., 9.,
|
|
2., 6., 10.,
|
|
3., 7., 11.,
|
|
4., 8., 12. };
|
|
|
|
double b[20] = { 1., 5., 9., 10.,
|
|
-2., -6., 10., 11.,
|
|
3., 7., 11., 12.,
|
|
4., 8., 12., 15.,
|
|
5., 9., 13., 14. };
|
|
|
|
double c[15] = { 39., 89., 139.,
|
|
30., 56., 82.,
|
|
49., 115., 181.,
|
|
58., 136., 214.,
|
|
59., 141., 223. };
|
|
|
|
double cnew[15];
|
|
|
|
qmckl_exit_code rc;
|
|
qmckl_matrix A = qmckl_matrix_alloc(context, 3, 4);
|
|
rc = qmckl_matrix_of_double(context, a, 12, &A);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
qmckl_matrix B = qmckl_matrix_alloc(context, 4, 5);
|
|
rc = qmckl_matrix_of_double(context, b, 20, &B);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
qmckl_matrix C = qmckl_matrix_alloc(context, 3, 5);
|
|
rc = qmckl_matmul(context, 'N', 'N', 0.5, A, B, 0., &C);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_double_of_matrix(context, C, cnew, 15);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
for (int i=0 ; i<15 ; ++i) {
|
|
printf("%f %f\n", cnew[i], c[i]);
|
|
assert (c[i] == cnew[i]);
|
|
}
|
|
}
|
|
#+end_src
|
|
** ~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
|
|
|
|
#+NAME: qmckl_adjugate_args
|
|
| 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
|
|
|
|
#+CALL: generate_c_header(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
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 );
|
|
#+end_src
|
|
|
|
For small matrices (\le 5\times 5), we use specialized routines
|
|
for performance motivations. For larger sizes, we rely on the
|
|
LAPACK library.
|
|
|
|
#+begin_src f90 :tangle (eval f) :exports none
|
|
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
|
|
|
|
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
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval f) :exports none
|
|
subroutine adjugate2(A,LDA,B,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
|
|
double precision :: C(2,2)
|
|
|
|
call cofactor2(A,LDA,C,2_8,na,det_l)
|
|
|
|
B(1,1) = C(1,1)
|
|
B(2,1) = C(1,2)
|
|
B(1,2) = C(2,1)
|
|
B(2,2) = C(2,2)
|
|
|
|
end subroutine adjugate2
|
|
|
|
subroutine adjugate3(a,LDA,B,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
|
|
double precision :: C(4,3)
|
|
|
|
call cofactor3(A,LDA,C,4_8,na,det_l)
|
|
|
|
B(1,1) = C(1,1)
|
|
B(1,2) = C(2,1)
|
|
B(1,3) = C(3,1)
|
|
B(2,1) = C(1,2)
|
|
B(2,2) = C(2,2)
|
|
B(2,3) = C(3,2)
|
|
B(3,1) = C(1,3)
|
|
B(3,2) = C(2,3)
|
|
B(3,3) = C(3,3)
|
|
|
|
end subroutine adjugate3
|
|
|
|
subroutine adjugate4(a,LDA,B,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
|
|
double precision :: C(4,4)
|
|
|
|
call cofactor4(A,LDA,C,4_8,na,det_l)
|
|
B(1,1) = C(1,1)
|
|
B(1,2) = C(2,1)
|
|
B(1,3) = C(3,1)
|
|
B(1,4) = C(4,1)
|
|
B(2,1) = C(1,2)
|
|
B(2,2) = C(2,2)
|
|
B(2,3) = C(3,2)
|
|
B(2,4) = C(4,2)
|
|
B(3,1) = C(1,3)
|
|
B(3,2) = C(2,3)
|
|
B(3,3) = C(3,3)
|
|
B(3,4) = C(4,3)
|
|
B(4,1) = C(1,4)
|
|
B(4,2) = C(2,4)
|
|
B(4,3) = C(3,4)
|
|
B(4,4) = C(4,4)
|
|
|
|
end subroutine adjugate4
|
|
|
|
subroutine adjugate5(A,LDA,B,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
|
|
double precision :: C(8,5)
|
|
|
|
call cofactor5(A,LDA,C,8_8,na,det_l)
|
|
|
|
B(1,1) = C(1,1)
|
|
B(1,2) = C(2,1)
|
|
B(1,3) = C(3,1)
|
|
B(1,4) = C(4,1)
|
|
B(1,5) = C(5,1)
|
|
B(2,1) = C(1,2)
|
|
B(2,2) = C(2,2)
|
|
B(2,3) = C(3,2)
|
|
B(2,4) = C(4,2)
|
|
B(2,5) = C(5,2)
|
|
B(3,1) = C(1,3)
|
|
B(3,2) = C(2,3)
|
|
B(3,3) = C(3,3)
|
|
B(3,4) = C(4,3)
|
|
B(3,5) = C(5,3)
|
|
B(4,1) = C(1,4)
|
|
B(4,2) = C(2,4)
|
|
B(4,3) = C(3,4)
|
|
B(4,4) = C(4,4)
|
|
B(4,5) = C(5,4)
|
|
B(5,1) = C(1,5)
|
|
B(5,2) = C(2,5)
|
|
B(5,3) = C(3,5)
|
|
B(5,4) = C(4,5)
|
|
B(5,5) = C(5,5)
|
|
|
|
end subroutine adjugate5
|
|
|
|
subroutine cofactor2(a,LDA,b,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8 :: na
|
|
double precision :: det_l
|
|
|
|
det_l = a(1,1)*a(2,2) - a(1,2)*a(2,1)
|
|
b(1,1) = a(2,2)
|
|
b(2,1) = -a(2,1)
|
|
b(1,2) = -a(1,2)
|
|
b(2,2) = a(1,1)
|
|
end subroutine cofactor2
|
|
|
|
subroutine cofactor3(a,LDA,b,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
integer :: i
|
|
|
|
det_l = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) &
|
|
-a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1)) &
|
|
+a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
|
|
|
|
b(1,1) = a(2,2)*a(3,3) - a(2,3)*a(3,2)
|
|
b(2,1) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
|
|
b(3,1) = a(2,1)*a(3,2) - a(2,2)*a(3,1)
|
|
|
|
b(1,2) = a(1,3)*a(3,2) - a(1,2)*a(3,3)
|
|
b(2,2) = a(1,1)*a(3,3) - a(1,3)*a(3,1)
|
|
b(3,2) = a(1,2)*a(3,1) - a(1,1)*a(3,2)
|
|
|
|
b(1,3) = a(1,2)*a(2,3) - a(1,3)*a(2,2)
|
|
b(2,3) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
|
|
b(3,3) = a(1,1)*a(2,2) - a(1,2)*a(2,1)
|
|
|
|
end subroutine cofactor3
|
|
|
|
subroutine cofactor4(a,LDA,b,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
integer :: i,j
|
|
det_l = a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) &
|
|
-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) &
|
|
+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))) &
|
|
-a(1,2)*(a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) &
|
|
-a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)) &
|
|
+a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))) &
|
|
+a(1,3)*(a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) &
|
|
-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)) &
|
|
+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))) &
|
|
-a(1,4)*(a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)) &
|
|
-a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)) &
|
|
+a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))
|
|
|
|
b(1,1) = a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))
|
|
b(2,1) = -a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))+a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))-a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))
|
|
b(3,1) = a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))
|
|
b(4,1) = -a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))+a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))-a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))
|
|
|
|
b(1,2) = -a(1,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))+a(1,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(1,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))
|
|
b(2,2) = a(1,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(1,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(1,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))
|
|
b(3,2) = -a(1,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))+a(1,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))-a(1,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))
|
|
b(4,2) = a(1,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))-a(1,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))+a(1,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))
|
|
|
|
b(1,3) = a(1,2)*(a(2,3)*a(4,4)-a(2,4)*a(4,3))-a(1,3)*(a(2,2)*a(4,4)-a(2,4)*a(4,2))+a(1,4)*(a(2,2)*a(4,3)-a(2,3)*a(4,2))
|
|
b(2,3) = -a(1,1)*(a(2,3)*a(4,4)-a(2,4)*a(4,3))+a(1,3)*(a(2,1)*a(4,4)-a(2,4)*a(4,1))-a(1,4)*(a(2,1)*a(4,3)-a(2,3)*a(4,1))
|
|
b(3,3) = a(1,1)*(a(2,2)*a(4,4)-a(2,4)*a(4,2))-a(1,2)*(a(2,1)*a(4,4)-a(2,4)*a(4,1))+a(1,4)*(a(2,1)*a(4,2)-a(2,2)*a(4,1))
|
|
b(4,3) = -a(1,1)*(a(2,2)*a(4,3)-a(2,3)*a(4,2))+a(1,2)*(a(2,1)*a(4,3)-a(2,3)*a(4,1))-a(1,3)*(a(2,1)*a(4,2)-a(2,2)*a(4,1))
|
|
|
|
b(1,4) = -a(1,2)*(a(2,3)*a(3,4)-a(2,4)*a(3,3))+a(1,3)*(a(2,2)*a(3,4)-a(2,4)*a(3,2))-a(1,4)*(a(2,2)*a(3,3)-a(2,3)*a(3,2))
|
|
b(2,4) = a(1,1)*(a(2,3)*a(3,4)-a(2,4)*a(3,3))-a(1,3)*(a(2,1)*a(3,4)-a(2,4)*a(3,1))+a(1,4)*(a(2,1)*a(3,3)-a(2,3)*a(3,1))
|
|
b(3,4) = -a(1,1)*(a(2,2)*a(3,4)-a(2,4)*a(3,2))+a(1,2)*(a(2,1)*a(3,4)-a(2,4)*a(3,1))-a(1,4)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
|
|
b(4,4) = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2))-a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1))+a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
|
|
|
|
end subroutine cofactor4
|
|
|
|
subroutine cofactor5(A,LDA,B,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
integer :: i,j
|
|
|
|
det_l = a(1,1)*(a(2,2)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*( &
|
|
a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))- &
|
|
a(2,3)*(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)- &
|
|
a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))+a(2,4)*(a(3,2)*( &
|
|
a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+ &
|
|
a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,5)*(a(3,2)*(a(4,3)*a(5,4)- &
|
|
a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)* &
|
|
a(5,3)-a(4,3)*a(5,2))))-a(1,2)*(a(2,1)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)* &
|
|
a(5,4))-a(3,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)- &
|
|
a(4,4)*a(5,3)))-a(2,3)*(a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*( &
|
|
a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))+ &
|
|
a(2,4)*(a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)- &
|
|
a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(2,5)*(a(3,1)*( &
|
|
a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+ &
|
|
a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))))+a(1,3)*(a(2,1)*(a(3,2)*(a(4,4)* &
|
|
a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*( &
|
|
a(4,2)*a(5,4)-a(4,4)*a(5,2)))-a(2,2)*(a(3,1)*(a(4,4)*a(5,5)-a(4,5)* &
|
|
a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)- &
|
|
a(4,4)*a(5,1)))+a(2,4)*(a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*( &
|
|
a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))- &
|
|
a(2,5)*(a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)- &
|
|
a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))-a(1,4)*(a(2,1)*( &
|
|
a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)* &
|
|
a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,2)*(a(3,1)*(a(4,3)* &
|
|
a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*( &
|
|
a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,3)*(a(3,1)*(a(4,2)*a(5,5)-a(4,5)* &
|
|
a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)- &
|
|
a(4,2)*a(5,1)))-a(2,5)*(a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*( &
|
|
a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))+ &
|
|
a(1,5)*(a(2,1)*(a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)* &
|
|
a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,2)*( &
|
|
a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)* &
|
|
a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,3)*(a(3,1)*(a(4,2)* &
|
|
a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*( &
|
|
a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(2,4)*(a(3,1)*(a(4,2)*a(5,3)-a(4,3)* &
|
|
a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)- &
|
|
a(4,2)*a(5,1))))
|
|
|
|
b(1,1) = &
|
|
(a(2,2)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))-a(2,3)* &
|
|
(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))+a(2,4)* &
|
|
(a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,5)* &
|
|
(a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))))
|
|
b(2,1) = &
|
|
(-a(2,1)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))+a(2,3)* &
|
|
(a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))-a(2,4)* &
|
|
(a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,5)* &
|
|
(a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))))
|
|
b(3,1) = &
|
|
(a(2,1)*(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))-a(2,2)* &
|
|
(a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))+a(2,4)* &
|
|
(a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(2,5)* &
|
|
(a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
b(4,1) = &
|
|
(-a(2,1)*(a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))+a(2,2)* &
|
|
(a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(2,3)* &
|
|
(a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))+a(2,5)* &
|
|
(a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
b(5,1) = &
|
|
(a(2,1)*(a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,2)* &
|
|
(a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,3)* &
|
|
(a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(2,4)* &
|
|
(a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
|
|
b(1,2) = &
|
|
(-a(1,2)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))+a(1,3)* &
|
|
(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))-a(1,4)* &
|
|
(a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))+a(1,5)* &
|
|
(a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))))
|
|
b(2,2) = &
|
|
(a(1,1)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))-a(1,3)* &
|
|
(a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))+a(1,4)* &
|
|
(a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(1,5)* &
|
|
(a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))))
|
|
b(3,2) = &
|
|
(-a(1,1)*(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))+a(1,2)* &
|
|
(a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))-a(1,4)* &
|
|
(a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))+a(1,5)* &
|
|
(a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
b(4,2) = &
|
|
(a(1,1)*(a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(1,2)* &
|
|
(a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(1,3)* &
|
|
(a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(1,5)* &
|
|
(a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
b(5,2) = &
|
|
(-a(1,1)*(a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))+a(1,2)* &
|
|
(a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(1,3)* &
|
|
(a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))+a(1,4)* &
|
|
(a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
|
|
b(1,3) = &
|
|
(a(1,2)*(a(2,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(2,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))-a(1,3)* &
|
|
(a(2,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))+a(1,4)* &
|
|
(a(2,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(1,5)* &
|
|
(a(2,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(2,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))))
|
|
b(2,3) = &
|
|
(-a(1,1)*(a(2,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(2,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))+a(1,3)* &
|
|
(a(2,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))-a(1,4)* &
|
|
(a(2,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(1,5)* &
|
|
(a(2,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))))
|
|
b(3,3) = &
|
|
(a(1,1)*(a(2,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))-a(1,2)* &
|
|
(a(2,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))+a(1,4)* &
|
|
(a(2,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(2,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(1,5)* &
|
|
(a(2,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(2,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
b(4,3) = &
|
|
(-a(1,1)*(a(2,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))+a(1,2)* &
|
|
(a(2,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(1,3)* &
|
|
(a(2,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(2,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))+a(1,5)* &
|
|
(a(2,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(2,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(2,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
b(5,3) = &
|
|
(a(1,1)*(a(2,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(2,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(1,2)* &
|
|
(a(2,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(1,3)* &
|
|
(a(2,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(2,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(1,4)* &
|
|
(a(2,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(2,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(2,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
|
|
b(1,4) = &
|
|
(-a(1,2)*(a(2,3)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))+a(2,5)*(a(3,3)*a(5,4)-a(3,4)*a(5,3)))+a(1,3)* &
|
|
(a(2,2)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,4)-a(3,4)*a(5,2)))-a(1,4)* &
|
|
(a(2,2)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,3)-a(3,3)*a(5,2)))+a(1,5)* &
|
|
(a(2,2)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))+a(2,4)*(a(3,2)*a(5,3)-a(3,3)*a(5,2))))
|
|
b(2,4) = &
|
|
(a(1,1)*(a(2,3)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))+a(2,5)*(a(3,3)*a(5,4)-a(3,4)*a(5,3)))-a(1,3)* &
|
|
(a(2,1)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,4)-a(3,4)*a(5,1)))+a(1,4)* &
|
|
(a(2,1)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,3)-a(3,3)*a(5,1)))-a(1,5)* &
|
|
(a(2,1)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,3)-a(3,3)*a(5,1))))
|
|
b(3,4) = &
|
|
(-a(1,1)*(a(2,2)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,4)-a(3,4)*a(5,2)))+a(1,2)* &
|
|
(a(2,1)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,4)-a(3,4)*a(5,1)))-a(1,4)* &
|
|
(a(2,1)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))-a(2,2)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,2)-a(3,2)*a(5,1)))+a(1,5)* &
|
|
(a(2,1)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))-a(2,2)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,2)-a(3,2)*a(5,1))))
|
|
b(4,4) = &
|
|
(a(1,1)*(a(2,2)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,3)-a(3,3)*a(5,2)))-a(1,2)* &
|
|
(a(2,1)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,3)-a(3,3)*a(5,1)))+a(1,3)* &
|
|
(a(2,1)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))-a(2,2)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,2)-a(3,2)*a(5,1)))-a(1,5)* &
|
|
(a(2,1)*(a(3,2)*a(5,3)-a(3,3)*a(5,2))-a(2,2)*(a(3,1)*a(5,3)-a(3,3)*a(5,1))+a(2,3)*(a(3,1)*a(5,2)-a(3,2)*a(5,1))))
|
|
b(5,4) = &
|
|
(-a(1,1)*(a(2,2)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))+a(2,4)*(a(3,2)*a(5,3)-a(3,3)*a(5,2)))+a(1,2)* &
|
|
(a(2,1)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,3)-a(3,3)*a(5,1)))-a(1,3)* &
|
|
(a(2,1)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))-a(2,2)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,2)-a(3,2)*a(5,1)))+a(1,4)* &
|
|
(a(2,1)*(a(3,2)*a(5,3)-a(3,3)*a(5,2))-a(2,2)*(a(3,1)*a(5,3)-a(3,3)*a(5,1))+a(2,3)*(a(3,1)*a(5,2)-a(3,2)*a(5,1))))
|
|
|
|
b(1,5) = &
|
|
(a(1,2)*(a(2,3)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))+a(2,5)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)))-a(1,3)* &
|
|
(a(2,2)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)))+a(1,4)* &
|
|
(a(2,2)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)))-a(1,5)* &
|
|
(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))))
|
|
b(2,5) = &
|
|
(-a(1,1)*(a(2,3)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))+a(2,5)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)))+a(1,3)* &
|
|
(a(2,1)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)))-a(1,4)* &
|
|
(a(2,1)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)))+a(1,5)* &
|
|
(a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))))
|
|
b(3,5) = &
|
|
(a(1,1)*(a(2,2)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)))-a(1,2)* &
|
|
(a(2,1)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)))+a(1,4)* &
|
|
(a(2,1)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))-a(2,2)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))-a(1,5)* &
|
|
(a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))))
|
|
b(4,5) = &
|
|
(-a(1,1)*(a(2,2)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)))+a(1,2)* &
|
|
(a(2,1)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)))-a(1,3)* &
|
|
(a(2,1)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))-a(2,2)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))+a(1,5)* &
|
|
(a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))-a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))+a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))))
|
|
b(5,5) = &
|
|
(a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)))-a(1,2)* &
|
|
(a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)))+a(1,3)* &
|
|
(a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))-a(1,4)* &
|
|
(a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))-a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))+a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))))
|
|
|
|
end
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_adjugate &
|
|
(context, n, A, lda, B, ldb, det_l) &
|
|
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 :: 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
|
|
real (c_double ) , intent(inout) :: det_l
|
|
|
|
integer(c_int32_t), external :: qmckl_adjugate_f
|
|
info = qmckl_adjugate_f &
|
|
(context, n, A, lda, B, ldb, det_l)
|
|
|
|
end function qmckl_adjugate
|
|
#+end_src
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_adjugate &
|
|
(context, n, A, lda, B, ldb, det_l) &
|
|
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 :: 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
|
|
real (c_double ) , intent(inout) :: det_l
|
|
|
|
end function qmckl_adjugate
|
|
end interface
|
|
#+end_src
|
|
|
|
|
|
#+begin_src f90 :tangle (eval 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
|
|
|
|
#+end_src
|
|
|
|
We first copy the array ~A~ into array ~B~.
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
B(1:na,1:na) = A(1:na,1:na)
|
|
#+end_src
|
|
|
|
Then, we compute the LU factorization $LU=A$
|
|
using the ~dgetrf~ routine.
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
call dgetrf(na, na, B, LDB, ipiv, inf )
|
|
#+end_src
|
|
|
|
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.
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
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
|
|
#+end_src
|
|
|
|
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.
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
if (iand(j,1_8) /= 0_8) then
|
|
det_l = -det_l
|
|
endif
|
|
#+end_src
|
|
|
|
Then, the inverse of $A$ is computed using ~dgetri~:
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
lwork = SIZE(work)
|
|
call dgetri(na, B, LDB, ipiv, work, lwork, inf )
|
|
#+end_src
|
|
|
|
and the adjugate matrix is computed as the product of the
|
|
determinant with the inverse:
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
B(:,:) = B(:,:)*det_l
|
|
|
|
end subroutine adjugate_general
|
|
#+end_src
|
|
|
|
*** Test :noexport:
|
|
|
|
#+begin_src f90 :tangle (eval f_test)
|
|
integer(qmckl_exit_code) function test_qmckl_adjugate(context) bind(C)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in), value :: context
|
|
|
|
double precision, allocatable :: A(:,:), B(:,:)
|
|
integer*8 :: m, n, k, LDA, LDB
|
|
integer*8 :: i,j,l
|
|
double precision :: x, det_l, det_l_ref
|
|
|
|
LDA = 6_8
|
|
LDB = 6_8
|
|
|
|
allocate( A(LDA,6), B(LDB,6))
|
|
|
|
A = 0.1d0
|
|
A(1,1) = 1.0d0;
|
|
A(2,2) = 2.0d0;
|
|
A(3,3) = 3.0d0;
|
|
A(4,4) = 4.0d0;
|
|
A(5,5) = 5.0d0;
|
|
A(6,6) = 6.0d0;
|
|
|
|
test_qmckl_adjugate = QMCKL_SUCCESS
|
|
|
|
#+end_src
|
|
|
|
#+begin_src python :results output :output drawer
|
|
import numpy as np
|
|
import numpy.linalg as la
|
|
N=6
|
|
|
|
A = np.zeros( (N,N) )
|
|
A += 0.1
|
|
for i in range(N):
|
|
A[i][i] = i+1
|
|
|
|
def adj(A):
|
|
return la.det(A) * la.inv(A)
|
|
|
|
print ("#+begin_src f90 :tangle (eval f_test)")
|
|
|
|
for i in range(N):
|
|
print ("! N = ", i+1)
|
|
print (f" test_qmckl_adjugate = qmckl_adjugate(context, {i+1}_8, A, LDA, B, LDB, det_l)")
|
|
print (f" if (test_qmckl_adjugate /= QMCKL_SUCCESS) return")
|
|
print (f" if (dabs((det_l - ({la.det(A[0:i+1,0:i+1])}d0))/det_l) > 1.d-13) then")
|
|
print (f" print *, 'N = {i+1}: det = ', det_l, {la.det(A[0:i+1,0:i+1])}d0")
|
|
print (f" test_qmckl_adjugate = {i+1}")
|
|
print (f" return")
|
|
print (f" end if")
|
|
print (f"")
|
|
print (f" x = 0.d0")
|
|
for j in range(i+1):
|
|
C = adj(A[0:i+1,0:i+1])
|
|
for k in range(i+1):
|
|
print (f" x = x + (B({j+1},{k+1}) - ({C[j,k]}) )**2")
|
|
print (f" if (dabs(x / det_l) > 1.d-12) then")
|
|
print (f" print *, 'N = {i+1}: x = ', x")
|
|
print (f" test_qmckl_adjugate = {i+1}")
|
|
print (f" return")
|
|
print (f" end if")
|
|
print (f"")
|
|
|
|
|
|
print ("#+end_src")
|
|
# print(adj(A[0:i+1,0:i+1]))
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
#+begin_example
|
|
,#+begin_src f90 :tangle (eval f_test)
|
|
! N = 1
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 1_8, A, LDA, B, LDB, det_l)
|
|
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
|
|
if (dabs((det_l - (1.0d0))/det_l) > 1.d-13) then
|
|
print *, 'N = 1: det = ', det_l, 1.0d0
|
|
test_qmckl_adjugate = 1
|
|
return
|
|
end if
|
|
|
|
x = 0.d0
|
|
x = x + (B(1,1) - (1.0) )**2
|
|
if (dabs(x / det_l) > 1.d-12) then
|
|
print *, 'N = 1: x = ', x
|
|
test_qmckl_adjugate = 1
|
|
return
|
|
end if
|
|
|
|
! N = 2
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 2_8, A, LDA, B, LDB, det_l)
|
|
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
|
|
if (dabs((det_l - (1.99d0))/det_l) > 1.d-13) then
|
|
print *, 'N = 2: det = ', det_l, 1.99d0
|
|
test_qmckl_adjugate = 2
|
|
return
|
|
end if
|
|
|
|
x = 0.d0
|
|
x = x + (B(1,1) - (1.9999999999999998) )**2
|
|
x = x + (B(1,2) - (-0.09999999999999999) )**2
|
|
x = x + (B(2,1) - (-0.09999999999999999) )**2
|
|
x = x + (B(2,2) - (0.9999999999999999) )**2
|
|
if (dabs(x / det_l) > 1.d-12) then
|
|
print *, 'N = 2: x = ', x
|
|
test_qmckl_adjugate = 2
|
|
return
|
|
end if
|
|
|
|
! N = 3
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 3_8, A, LDA, B, LDB, det_l)
|
|
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
|
|
if (dabs((det_l - (5.942000000000001d0))/det_l) > 1.d-13) then
|
|
print *, 'N = 3: det = ', det_l, 5.942000000000001d0
|
|
test_qmckl_adjugate = 3
|
|
return
|
|
end if
|
|
|
|
x = 0.d0
|
|
x = x + (B(1,1) - (5.990000000000001) )**2
|
|
x = x + (B(1,2) - (-0.29000000000000004) )**2
|
|
x = x + (B(1,3) - (-0.19000000000000003) )**2
|
|
x = x + (B(2,1) - (-0.29000000000000004) )**2
|
|
x = x + (B(2,2) - (2.9900000000000007) )**2
|
|
x = x + (B(2,3) - (-0.09000000000000001) )**2
|
|
x = x + (B(3,1) - (-0.19000000000000003) )**2
|
|
x = x + (B(3,2) - (-0.09) )**2
|
|
x = x + (B(3,3) - (1.9900000000000002) )**2
|
|
if (dabs(x / det_l) > 1.d-12) then
|
|
print *, 'N = 3: x = ', x
|
|
test_qmckl_adjugate = 3
|
|
return
|
|
end if
|
|
|
|
! N = 4
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 4_8, A, LDA, B, LDB, det_l)
|
|
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
|
|
if (dabs((det_l - (23.669700000000006d0))/det_l) > 1.d-13) then
|
|
print *, 'N = 4: det = ', det_l, 23.669700000000006d0
|
|
test_qmckl_adjugate = 4
|
|
return
|
|
end if
|
|
|
|
x = 0.d0
|
|
x = x + (B(1,1) - (23.91200000000001) )**2
|
|
x = x + (B(1,2) - (-1.1310000000000004) )**2
|
|
x = x + (B(1,3) - (-0.7410000000000001) )**2
|
|
x = x + (B(1,4) - (-0.5510000000000002) )**2
|
|
x = x + (B(2,1) - (-1.1310000000000002) )**2
|
|
x = x + (B(2,2) - (11.922000000000002) )**2
|
|
x = x + (B(2,3) - (-0.351) )**2
|
|
x = x + (B(2,4) - (-0.261) )**2
|
|
x = x + (B(3,1) - (-0.7410000000000002) )**2
|
|
x = x + (B(3,2) - (-0.351) )**2
|
|
x = x + (B(3,3) - (7.932000000000001) )**2
|
|
x = x + (B(3,4) - (-0.17100000000000004) )**2
|
|
x = x + (B(4,1) - (-0.5510000000000002) )**2
|
|
x = x + (B(4,2) - (-0.261) )**2
|
|
x = x + (B(4,3) - (-0.17100000000000004) )**2
|
|
x = x + (B(4,4) - (5.942000000000001) )**2
|
|
if (dabs(x / det_l) > 1.d-12) then
|
|
print *, 'N = 4: x = ', x
|
|
test_qmckl_adjugate = 4
|
|
return
|
|
end if
|
|
|
|
! N = 5
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 5_8, A, LDA, B, LDB, det_l)
|
|
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
|
|
if (dabs((det_l - (117.91554000000008d0))/det_l) > 1.d-13) then
|
|
print *, 'N = 5: det = ', det_l, 117.91554000000008d0
|
|
test_qmckl_adjugate = 5
|
|
return
|
|
end if
|
|
|
|
x = 0.d0
|
|
x = x + (B(1,1) - (119.31770000000006) )**2
|
|
x = x + (B(1,2) - (-5.541900000000004) )**2
|
|
x = x + (B(1,3) - (-3.6309000000000022) )**2
|
|
x = x + (B(1,4) - (-2.6999000000000017) )**2
|
|
x = x + (B(1,5) - (-2.1489000000000016) )**2
|
|
x = x + (B(2,1) - (-5.541900000000004) )**2
|
|
x = x + (B(2,2) - (59.435700000000026) )**2
|
|
x = x + (B(2,3) - (-1.7199000000000007) )**2
|
|
x = x + (B(2,4) - (-1.2789000000000006) )**2
|
|
x = x + (B(2,5) - (-1.0179000000000005) )**2
|
|
x = x + (B(3,1) - (-3.6309000000000027) )**2
|
|
x = x + (B(3,2) - (-1.7199000000000007) )**2
|
|
x = x + (B(3,3) - (39.53370000000002) )**2
|
|
x = x + (B(3,4) - (-0.8379000000000005) )**2
|
|
x = x + (B(3,5) - (-0.6669000000000004) )**2
|
|
x = x + (B(4,1) - (-2.699900000000002) )**2
|
|
x = x + (B(4,2) - (-1.2789000000000006) )**2
|
|
x = x + (B(4,3) - (-0.8379000000000004) )**2
|
|
x = x + (B(4,4) - (29.611700000000017) )**2
|
|
x = x + (B(4,5) - (-0.4959000000000003) )**2
|
|
x = x + (B(5,1) - (-2.1489000000000016) )**2
|
|
x = x + (B(5,2) - (-1.0179000000000005) )**2
|
|
x = x + (B(5,3) - (-0.6669000000000004) )**2
|
|
x = x + (B(5,4) - (-0.4959000000000003) )**2
|
|
x = x + (B(5,5) - (23.669700000000013) )**2
|
|
if (dabs(x / det_l) > 1.d-12) then
|
|
print *, 'N = 5: x = ', x
|
|
test_qmckl_adjugate = 5
|
|
return
|
|
end if
|
|
|
|
! N = 6
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 6_8, A, LDA, B, LDB, det_l)
|
|
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
|
|
if (dabs((det_l - (705.1783350000001d0))/det_l) > 1.d-13) then
|
|
print *, 'N = 6: det = ', det_l, 705.1783350000001d0
|
|
test_qmckl_adjugate = 6
|
|
return
|
|
end if
|
|
|
|
x = 0.d0
|
|
x = x + (B(1,1) - (714.5040400000001) )**2
|
|
x = x + (B(1,2) - (-32.697210000000005) )**2
|
|
x = x + (B(1,3) - (-21.422310000000003) )**2
|
|
x = x + (B(1,4) - (-15.929410000000006) )**2
|
|
x = x + (B(1,5) - (-12.678510000000003) )**2
|
|
x = x + (B(1,6) - (-10.529610000000003) )**2
|
|
x = x + (B(2,1) - (-32.69721) )**2
|
|
x = x + (B(2,2) - (355.65834) )**2
|
|
x = x + (B(2,3) - (-10.147409999999997) )**2
|
|
x = x + (B(2,4) - (-7.54551) )**2
|
|
x = x + (B(2,5) - (-6.005610000000001) )**2
|
|
x = x + (B(2,6) - (-4.987709999999999) )**2
|
|
x = x + (B(3,1) - (-21.422310000000003) )**2
|
|
x = x + (B(3,2) - (-10.147409999999999) )**2
|
|
x = x + (B(3,3) - (236.51663999999997) )**2
|
|
x = x + (B(3,4) - (-4.943610000000001) )**2
|
|
x = x + (B(3,5) - (-3.93471) )**2
|
|
x = x + (B(3,6) - (-3.267810000000001) )**2
|
|
x = x + (B(4,1) - (-15.929410000000003) )**2
|
|
x = x + (B(4,2) - (-7.54551) )**2
|
|
x = x + (B(4,3) - (-4.9436100000000005) )**2
|
|
x = x + (B(4,4) - (177.13894000000002) )**2
|
|
x = x + (B(4,5) - (-2.92581) )**2
|
|
x = x + (B(4,6) - (-2.42991) )**2
|
|
x = x + (B(5,1) - (-12.678510000000001) )**2
|
|
x = x + (B(5,2) - (-6.005609999999999) )**2
|
|
x = x + (B(5,3) - (-3.9347100000000004) )**2
|
|
x = x + (B(5,4) - (-2.92581) )**2
|
|
x = x + (B(5,5) - (141.58524) )**2
|
|
x = x + (B(5,6) - (-1.93401) )**2
|
|
x = x + (B(6,1) - (-10.529610000000003) )**2
|
|
x = x + (B(6,2) - (-4.98771) )**2
|
|
x = x + (B(6,3) - (-3.2678100000000003) )**2
|
|
x = x + (B(6,4) - (-2.42991) )**2
|
|
x = x + (B(6,5) - (-1.9340100000000005) )**2
|
|
x = x + (B(6,6) - (117.91554000000001) )**2
|
|
if (dabs(x / det_l) > 1.d-12) then
|
|
print *, 'N = 6: x = ', x
|
|
test_qmckl_adjugate = 6
|
|
return
|
|
end if
|
|
|
|
,#+end_src
|
|
#+end_example
|
|
|
|
#+begin_src f90 :tangle (eval f_test)
|
|
|
|
deallocate(A,B)
|
|
|
|
end function test_qmckl_adjugate
|
|
#+end_src
|
|
|
|
#+begin_src c :comments link :tangle (eval c_test)
|
|
qmckl_exit_code test_qmckl_adjugate(qmckl_context context);
|
|
assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
|
|
#+end_src
|
|
|
|
** ~qmckl_adjugate_safe~
|
|
|
|
"Size-safe" proxy function with the same functionality as ~qmckl_adjugate~
|
|
but with 2 additional arguments. These arguments ~size_max_M~ (where ~M~ is a matix)
|
|
are required primarily for the Python API, where compatibility with
|
|
NumPy arrays implies that sizes of the input and output arrays are provided.
|
|
|
|
|
|
#+NAME: qmckl_adjugate_safe_args
|
|
| 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$ |
|
|
| ~size_max_A~ | ~int64_t~ | in | Size of the matrix $A$ |
|
|
| ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ |
|
|
| ~B~ | ~double[][ldb]~ | out | Adjugate of $A$ |
|
|
| ~size_max_B~ | ~int64_t~ | in | Size of the matrix $B$ |
|
|
| ~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
|
|
- ~size_max_A >= m * m~
|
|
- ~size_max_B >= m * m~
|
|
|
|
#+CALL: generate_c_header(table=qmckl_adjugate_safe_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate_safe")
|
|
|
|
#+RESULTS:
|
|
#+BEGIN_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_adjugate_safe (
|
|
const qmckl_context context,
|
|
const int64_t n,
|
|
const double* A,
|
|
const int64_t size_max_A,
|
|
const int64_t lda,
|
|
double* const B,
|
|
const int64_t size_max_B,
|
|
const int64_t ldb,
|
|
double* det_l );
|
|
#+END_src
|
|
|
|
For small matrices (\le 5\times 5), we use specialized routines
|
|
for performance motivations. For larger sizes, we rely on the
|
|
LAPACK library.
|
|
|
|
#+begin_src f90 :tangle (eval f) :exports none
|
|
integer function qmckl_adjugate_safe_f(context, &
|
|
na, A, size_A, LDA, B, size_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) :: size_A
|
|
integer*8, intent(in) :: LDA
|
|
double precision, intent(out) :: B (LDB,*)
|
|
integer*8, intent(in) :: size_B
|
|
integer*8, intent(in) :: LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
|
|
integer, external :: qmckl_adjugate_f
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (size_A < na) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
if (size_B <= 0_8) then
|
|
info = QMCKL_INVALID_ARG_7
|
|
return
|
|
endif
|
|
|
|
info = qmckl_adjugate_f(context, na, A, LDA, B, LDB, det_l)
|
|
|
|
if (info == QMCKL_INVALID_ARG_4) then
|
|
info = QMCKL_INVALID_ARG_5
|
|
return
|
|
endif
|
|
|
|
if (info == QMCKL_INVALID_ARG_6) then
|
|
info = QMCKL_INVALID_ARG_8
|
|
return
|
|
endif
|
|
|
|
end function qmckl_adjugate_safe_f
|
|
#+end_src
|
|
|
|
*** C interface
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_adjugate_safe_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate_safe")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_adjugate_safe &
|
|
(context, n, A, size_A, lda, B, size_B, ldb, det_l) &
|
|
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 :: n
|
|
real (c_double ) , intent(in) :: A(lda,*)
|
|
integer (c_int64_t) , intent(in) , value :: lda
|
|
integer (c_int64_t) , intent(in) , value :: size_A
|
|
real (c_double ) , intent(out) :: B(ldb,*)
|
|
integer (c_int64_t) , intent(in) , value :: ldb
|
|
integer (c_int64_t) , intent(in) , value :: size_B
|
|
real (c_double ) , intent(inout) :: det_l
|
|
|
|
integer(c_int32_t), external :: qmckl_adjugate_safe_f
|
|
info = qmckl_adjugate_safe_f &
|
|
(context, n, A, size_A, lda, B, size_B, ldb, det_l)
|
|
|
|
end function qmckl_adjugate_safe
|
|
#+end_src
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_adjugate_safe_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate_safe")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_adjugate_safe &
|
|
(context, n, A, size_A, lda, B, size_B, ldb, det_l) &
|
|
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 :: n
|
|
real (c_double ) , intent(in) :: A(lda,*)
|
|
integer (c_int64_t) , intent(in) , value :: lda
|
|
integer (c_int64_t) , intent(in) , value :: size_A
|
|
real (c_double ) , intent(out) :: B(ldb,*)
|
|
integer (c_int64_t) , intent(in) , value :: ldb
|
|
integer (c_int64_t) , intent(in) , value :: size_B
|
|
real (c_double ) , intent(inout) :: det_l
|
|
|
|
end function qmckl_adjugate_safe
|
|
end interface
|
|
#+end_src
|
|
|
|
** ~qmckl_transpose~
|
|
|
|
Transposes a matrix: $A^\dagger_{ji} = A_{ij}$.
|
|
|
|
| Variable | Type | In/Out | Description |
|
|
|-----------+-----------------+--------+-------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~A~ | ~qmckl_matrix~ | in | Input matrix |
|
|
| ~At~ | ~qmckl_matrix~ | out | Transposed matrix |
|
|
|
|
#+begin_src c :tangle (eval h_private_func) :comments org
|
|
qmckl_exit_code
|
|
qmckl_transpose (qmckl_context context,
|
|
const qmckl_matrix A,
|
|
qmckl_matrix At );
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :tangle (eval c) :comments org :exports none
|
|
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;
|
|
}
|
|
#+end_src
|
|
|
|
*** Test :noexport:
|
|
|
|
#+begin_src c :comments link :tangle (eval c_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);
|
|
}
|
|
|
|
#+end_src
|
|
|
|
* End of files :noexport:
|
|
|
|
|
|
#+begin_src c :tangle (eval h_private_type)
|
|
#endif
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_func)
|
|
#endif
|
|
#+end_src
|
|
|
|
#+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
|