test_lapack.c

Go to the documentation of this file.

// Copyright (c) 2026 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
 
#include "tests.h"
#include "yac_lapack_interface.h"
 
static int utest_nearly_equal(double a, double b, double tol) {
  return fabs(a - b) <= tol;
}
 
static void utest_check_identity_2x2(
  double const a_col_major[4], double const inv_col_major[4],
  double tol, char const * test_name) {
 
  // C = A * A^{-1} in column-major layout
  double c[4] = {
    a_col_major[0] * inv_col_major[0] + a_col_major[2] * inv_col_major[1],
    a_col_major[1] * inv_col_major[0] + a_col_major[3] * inv_col_major[1],
    a_col_major[0] * inv_col_major[2] + a_col_major[2] * inv_col_major[3],
    a_col_major[1] * inv_col_major[2] + a_col_major[3] * inv_col_major[3]
  };
 
  if (!utest_nearly_equal(c[0], 1.0, tol) ||
      !utest_nearly_equal(c[1], 0.0, tol) ||
      !utest_nearly_equal(c[2], 0.0, tol) ||
      !utest_nearly_equal(c[3], 1.0, tol)) {
    PUT_ERR_F("ERROR in %s: A*Ainv != I", test_name);
  }
}
 
static void utest_symmetrise_from_upper_2x2(double a_col_major[4]) {
 
  a_col_major[1] = a_col_major[2];
}
 
int main(void) {
 
  puts("----- yac_lapack_interface basic tests -----");
 
  {
    // A = [[4,2],[1,3]], x_ref = [1,2], b = A*x_ref = [8,7]
    double a[4] = {4.0, 1.0, 2.0, 3.0};
    double b[2] = {8.0, 7.0};
    lapack_int ipiv[2];
 
    lapack_int info =
      LAPACKE_dgesv(LAPACK_COL_MAJOR, 2, 1, a, 2, ipiv, b, 2);
 
    if (info != 0) PUT_ERR_F("ERROR in LAPACKE_dgesv: info=%d", (int)info);
 
    if (!utest_nearly_equal(b[0], 1.0, 1e-12) ||
        !utest_nearly_equal(b[1], 2.0, 1e-12))
      PUT_ERR_F(
        "ERROR in LAPACKE_dgesv: wrong solution [%1.16e, %1.16e]", b[0], b[1]);
  }
 
  {
    double a_ref[4] = {4.0, 1.0, 2.0, 3.0};
    double a_inv[4] = {4.0, 1.0, 2.0, 3.0};
    lapack_int ipiv[2];
    double work[16];
 
    lapack_int info = LAPACKE_dgetrf(LAPACK_COL_MAJOR, 2, 2, a_inv, 2, ipiv);
 
    if (info != 0) PUT_ERR_F("ERROR in LAPACKE_dgetrf: info=%d", (int)info);
 
    if (info == 0) {
      info = LAPACKE_dgetri_work(
        LAPACK_COL_MAJOR, 2, a_inv, 2, ipiv, work, 16);
      if (info != 0)
        PUT_ERR_F("ERROR in LAPACKE_dgetri_work: info=%d", (int)info);
    }
 
    if (info == 0) utest_check_identity_2x2(a_ref, a_inv, 1e-12, "dgetrf+dgetri");
  }
 
#ifndef YAC_LAPACK_NO_DSYTR
  {
    // symmetric matrix A = [[4,1],[1,3]]
    double a_ref[4] = {4.0, 1.0, 1.0, 3.0};
    double a_inv[4] = {4.0, 1.0, 1.0, 3.0};
    lapack_int ipiv[2];
    double work[16];
 
    lapack_int info =
      LAPACKE_dsytrf_work(LAPACK_COL_MAJOR, 'U', 2, a_inv, 2, ipiv, work, 16);
 
    if (info != 0)
      PUT_ERR_F("ERROR in LAPACKE_dsytrf_work: info=%d", (int)info);
 
    if (info == 0) {
      info = LAPACKE_dsytri_work(
        LAPACK_COL_MAJOR, 'U', 2, a_inv, 2, ipiv, work);
      if (info != 0)
        PUT_ERR_F("ERROR in LAPACKE_dsytri_work: info=%d", (int)info);
    }
 
    if (info == 0) {
      utest_symmetrise_from_upper_2x2(a_inv);
      utest_check_identity_2x2(a_ref, a_inv, 1e-12, "dsytrf+dsytri");
    }
  }
#endif
 
  return TEST_EXIT_CODE;
}