arael 0.3.0

Nonlinear optimization framework with compile-time symbolic differentiation
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170

// Eigen sparse Cholesky bridge for arael.
//
// Provides extern "C" functions wrapping:
// - Eigen::SimplicialLLT (f32 and f64)
// - Eigen::CholmodSimplicialLLT (f64 only, guarded by ARAEL_CHOLMOD define)
//
// CSC format: upper triangle, col_ptr as int64_t (Rust usize), row_idx as int32_t (Rust u32).
// Eigen uses int (32-bit) StorageIndex by default; col_ptr is narrowed once and cached.

#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
#include <cstdint>
#include <vector>

#ifdef ARAEL_CHOLMOD
#include <Eigen/CholmodSupport>
#endif

// -- Templated SimplicialLLT handle --

template<typename Scalar>
struct LLTHandle {
    using SpMat = Eigen::SparseMatrix<Scalar, Eigen::ColMajor, int>;
    using Solver = Eigen::SimplicialLLT<SpMat, Eigen::Upper>;

    Solver solver;
    bool pattern_analyzed;
    std::vector<int> col_ptr_i32;

    LLTHandle() : pattern_analyzed(false) {}
};

template<typename Scalar>
static bool llt_solve(
    void* handle_ptr,
    int n, int nnz,
    const int64_t* col_ptr,
    const int32_t* row_idx,
    const Scalar* vals,
    const Scalar* rhs,
    Scalar* solution)
{
    try {
        auto* h = static_cast<LLTHandle<Scalar>*>(handle_ptr);
        using SpMat = typename LLTHandle<Scalar>::SpMat;

        // Convert col_ptr from int64 to int32 (cached allocation)
        h->col_ptr_i32.resize(n + 1);
        for (int i = 0; i <= n; i++) {
            h->col_ptr_i32[i] = static_cast<int>(col_ptr[i]);
        }

        // Map CSC data into Eigen sparse matrix (no copy of vals/row_idx)
        Eigen::Map<const SpMat> mat(
            n, n, nnz,
            h->col_ptr_i32.data(),
            reinterpret_cast<const int*>(row_idx),
            vals);

        if (!h->pattern_analyzed) {
            h->solver.analyzePattern(mat);
            if (h->solver.info() != Eigen::Success) return false;
            h->pattern_analyzed = true;
        }

        h->solver.factorize(mat);
        if (h->solver.info() != Eigen::Success) return false;

        Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>> b(rhs, n);
        Eigen::Matrix<Scalar, Eigen::Dynamic, 1> x = h->solver.solve(b);
        if (h->solver.info() != Eigen::Success) return false;
        Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, 1>>(solution, n) = x;

        return true;
    } catch (...) {
        return false;
    }
}

// -- CholmodSimplicialLLT handle (f64 only) --

#ifdef ARAEL_CHOLMOD
struct CholmodHandle {
    using SpMat = Eigen::SparseMatrix<double, Eigen::ColMajor, int>;
    using Solver = Eigen::CholmodSimplicialLLT<SpMat, Eigen::Upper>;

    Solver solver;
    bool pattern_analyzed;
    std::vector<int> col_ptr_i32;

    CholmodHandle() : pattern_analyzed(false) {}
};

static bool cholmod_solve_impl(
    void* handle_ptr,
    int n, int nnz,
    const int64_t* col_ptr,
    const int32_t* row_idx,
    const double* vals,
    const double* rhs,
    double* solution)
{
    try {
        auto* h = static_cast<CholmodHandle*>(handle_ptr);
        using SpMat = CholmodHandle::SpMat;

        h->col_ptr_i32.resize(n + 1);
        for (int i = 0; i <= n; i++) {
            h->col_ptr_i32[i] = static_cast<int>(col_ptr[i]);
        }

        Eigen::Map<const SpMat> mat(
            n, n, nnz,
            h->col_ptr_i32.data(),
            reinterpret_cast<const int*>(row_idx),
            vals);

        if (!h->pattern_analyzed) {
            h->solver.analyzePattern(mat);
            h->pattern_analyzed = true;
        }

        h->solver.factorize(mat);
        if (h->solver.info() != Eigen::Success) return false;

        Eigen::Map<const Eigen::VectorXd> b(rhs, n);
        Eigen::Map<Eigen::VectorXd> x(solution, n);
        x = h->solver.solve(b);

        return h->solver.info() == Eigen::Success;
    } catch (...) {
        return false;
    }
}
#endif

// -- extern "C" interface --

extern "C" {

// f64 SimplicialLLT
void* eigen_llt_f64_create() { return new LLTHandle<double>(); }
void eigen_llt_f64_destroy(void* h) { delete static_cast<LLTHandle<double>*>(h); }
bool eigen_llt_f64_solve(void* h, int n, int nnz,
    const int64_t* cp, const int32_t* ri, const double* v,
    const double* rhs, double* sol) {
    return llt_solve<double>(h, n, nnz, cp, ri, v, rhs, sol);
}

// f32 SimplicialLLT
void* eigen_llt_f32_create() { return new LLTHandle<float>(); }
void eigen_llt_f32_destroy(void* h) { delete static_cast<LLTHandle<float>*>(h); }
bool eigen_llt_f32_solve(void* h, int n, int nnz,
    const int64_t* cp, const int32_t* ri, const float* v,
    const float* rhs, float* sol) {
    return llt_solve<float>(h, n, nnz, cp, ri, v, rhs, sol);
}

// f64 CholmodSimplicialLLT
#ifdef ARAEL_CHOLMOD
void* eigen_cholmod_f64_create() { return new CholmodHandle(); }
void eigen_cholmod_f64_destroy(void* h) { delete static_cast<CholmodHandle*>(h); }
bool eigen_cholmod_f64_solve(void* h, int n, int nnz,
    const int64_t* cp, const int32_t* ri, const double* v,
    const double* rhs, double* sol) {
    return cholmod_solve_impl(h, n, nnz, cp, ri, v, rhs, sol);
}
#endif

} // extern "C"