fsb_linalg.h

#pragma once
#include <math.h>
#include <stddef.h>
#include <stdbool.h>
 
#ifdef __cplusplus
#include <array>
#endif
 
#include "openblas/lapack.h"
 
#ifdef __cplusplus
extern "C"
{
#endif
 
#define FSB_MAX(a, b) (((a) > (b)) ? (a) : (b))
#define FSB_MIN(a, b) (((a) < (b)) ? (a) : (b))
 
typedef double double_t;
 
typedef enum FsbLapackErrorType
{
    EFSB_LAPACK_ERROR_NONE = 0,
    EFSB_LAPACK_ERROR_INPUT = 1,
    EFSB_LAPACK_ERROR_MEMORY = 2,
    EFSB_LAPACK_ERROR_CONVERGE = 3,
    EFSB_LAPACK_ERROR_QUERY = 4,
    EFSB_LAPACK_NOT_POSITIVE_DEFINITE = 5,
    EFSB_LAPACK_SINGULAR = 6,
    EFSB_LAPACK_ERROR_NOT_FULL_RANK = 7, 
} FsbLinalgErrorType;
 
FsbLinalgErrorType fsb_linalg_svd(
    const double_t mat[], size_t rows, size_t cols, bool u_full, bool v_full, size_t work_len, double_t work[], double_t unitary_u[],
    double_t sing_val[], double_t unitary_vt[]);
 
FsbLinalgErrorType fsb_linalg_matrix_eig(
    const double_t mat[], size_t dim, size_t work_len, double_t work[],
    double_t val_real[], double_t val_imag[],
    double_t vec_real[], double_t vec_imag[]);
 
FsbLinalgErrorType fsb_linalg_sym_lt_eig(
    const double_t mat[], size_t dim, size_t work_len, double_t work[],
    double_t val[], double_t vec[]);
 
FsbLinalgErrorType fsb_linalg_cholesky_decomposition(
    const double_t mat[], size_t dim, double_t mat_chol[]);
 
FsbLinalgErrorType fsb_linalg_cholesky_solve(const double_t mat[], const double_t b_vec[],
                                             size_t nrhs, size_t dim,
                                             size_t work_len, double_t work[],
                                             double_t x_vec[]);
 
bool fsb_linalg_is_posdef(const double_t mat[], size_t dim, size_t work_len, double_t work[]);
 
FsbLinalgErrorType fsb_linalg_matrix_sqr_solve(
    const double_t mat[], const double_t y_vec[], size_t nrhs, size_t dim, size_t work_len,
    size_t iwork_len, double_t work[], lapack_int iwork[], double_t x_vec[]);
 
FsbLinalgErrorType fsb_linalg_pseudoinverse(const double_t mat[], size_t rows, size_t columns, size_t work_len, double_t work[], double_t inv_mat[]);
 
FsbLinalgErrorType fsb_linalg_leastsquares_solve(
    const double_t mat[], size_t rows, size_t columns, const double_t b_vec[], size_t nrhs,
    size_t work_len, double_t work[], double_t x_vec[]);
 
void sample_dgels_example(void);
 
#ifdef __cplusplus
}
 
template <size_t Rows, size_t Cols, size_t WorkLen>
inline FsbLinalgErrorType fsb_linalg_svd_array(
    const std::array<double_t, Rows * Cols>& mat,
    bool u_full,
    bool v_full,
    std::array<double_t, WorkLen>& work,
    std::array<double_t, Rows * Rows>& unitary_u,
    std::array<double_t, FSB_MIN(Rows, Cols)>& sing_val,
    std::array<double_t, Cols * Cols>& unitary_vt)
{
    return fsb_linalg_svd(
        mat.data(), Rows, Cols,
        u_full, v_full,
        WorkLen, work.data(),
        unitary_u.data(), sing_val.data(), unitary_vt.data());
}
 
template <size_t Dim, size_t WorkLen>
inline FsbLinalgErrorType fsb_linalg_matrix_eig_array(
    const std::array<double_t, Dim * Dim>& mat,
    std::array<double_t, WorkLen>& work,
    std::array<double_t, Dim>& val_real,
    std::array<double_t, Dim>& val_imag,
    std::array<double_t, Dim * Dim>& vec_real,
    std::array<double_t, Dim * Dim>& vec_imag)
{
    return fsb_linalg_matrix_eig(
        mat.data(), Dim, WorkLen, work.data(),
        val_real.data(), val_imag.data(),
        vec_real.data(), vec_imag.data());
}
 
template <size_t Dim, size_t WorkLen>
inline FsbLinalgErrorType fsb_linalg_sym_lt_eig_array(
    const std::array<double_t, Dim * Dim>& mat,
    std::array<double_t, WorkLen>& work,
    std::array<double_t, Dim>& val,
    std::array<double_t, Dim * Dim>& vec)
{
    return fsb_linalg_sym_lt_eig(
        mat.data(), Dim, WorkLen, work.data(),
        val.data(), vec.data());
}
 
template <size_t Dim>
inline FsbLinalgErrorType fsb_linalg_cholesky_decomposition_array(
    const std::array<double_t, Dim * Dim>& mat,
    std::array<double_t, Dim * Dim>& mat_chol)
{
    return fsb_linalg_cholesky_decomposition(mat.data(), Dim, mat_chol.data());
}
 
template <size_t Dim, size_t WorkLen>
inline bool fsb_linalg_is_posdef_array(
    const std::array<double_t, Dim * Dim>& mat,
    std::array<double_t, WorkLen>& work)
{
    return fsb_linalg_is_posdef(mat.data(), Dim, WorkLen, work.data());
}
 
template <size_t Dim, size_t Nrhs, size_t WorkLen>
inline FsbLinalgErrorType fsb_linalg_cholesky_solve_array(
    const std::array<double_t, Dim * Dim>& mat,
    const std::array<double_t, Dim * Nrhs>& b_vec,
    std::array<double_t, WorkLen>& work,
    std::array<double_t, Dim * Nrhs>& x_vec)
{
    return fsb_linalg_cholesky_solve(
        mat.data(), b_vec.data(),
        Nrhs, Dim,
        WorkLen, work.data(),
        x_vec.data());
}
 
template <size_t Dim, size_t Nrhs, size_t WorkLen, size_t IWorkLen>
inline FsbLinalgErrorType fsb_linalg_matrix_sqr_solve_array(
    const std::array<double_t, Dim * Dim>& mat,
    const std::array<double_t, Dim * Nrhs>& y_vec,
    std::array<double_t, WorkLen>& work,
    std::array<lapack_int, IWorkLen>& iwork,
    std::array<double_t, Dim * Nrhs>& x_vec)
{
    return fsb_linalg_matrix_sqr_solve(
        mat.data(), y_vec.data(),
        Nrhs, Dim,
        WorkLen, IWorkLen,
        work.data(), iwork.data(),
        x_vec.data());
}
 
template <size_t Rows, size_t Cols, size_t WorkLen>
inline FsbLinalgErrorType fsb_linalg_pseudoinverse_array(
    const std::array<double_t, Rows * Cols>& mat,
    std::array<double_t, WorkLen>& work,
    std::array<double_t, Cols * Rows>& inv_mat)
{
    return fsb_linalg_pseudoinverse(
        mat.data(), Rows, Cols,
        WorkLen, work.data(),
        inv_mat.data());
}
 
template <size_t Rows, size_t Cols, size_t Nrhs, size_t WorkLen>
inline FsbLinalgErrorType fsb_linalg_leastsquares_solve_array(
    const std::array<double_t, Rows * Cols>& mat,
    const std::array<double_t, Rows * Nrhs>& b_vec,
    std::array<double_t, WorkLen>& work,
    std::array<double_t, Cols * Nrhs>& x_vec)
{
    return fsb_linalg_leastsquares_solve(
        mat.data(), Rows, Cols,
        b_vec.data(), Nrhs,
        WorkLen, work.data(),
        x_vec.data());
}
 
#endif