77
88#include < array>
99#include < cstddef>
10- #include < type_traits>
1110
1211#include < GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp>
1312#include < GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp>
@@ -22,56 +21,61 @@ namespace GridKit
2221 /* *
2322 * @brief Enzyme automatic differentiation Hessian evaluator: Lagrangian Hessian, d2L/dx2
2423 *
25- * Lower triangle of the Hessian of
26- * \f$L = \sigma f(x) + \lambda^T g(x)\f$ by one forward sweep per local
27- * variable over the reverse-mode gradient of \f$L\f$. Only structural
28- * entries are stored, in global indices.
24+ * Hessian of \f$L = \sigma f(x) + \lambda_c^T c(x) + \lambda_d^T d(x)\f$
25+ * by one forward sweep per local variable over the reverse-mode gradient
26+ * of \f$L\f$. Only structural entries are stored, in global indices.
2927 *
30- * @tparam ModelT - model type with `VARIABLE_SIZE` local variables and
31- * `CONSTRAINT_SIZE` local constraints
28+ * @tparam ModelT - model type with internal and external local variables
29+ * and equality constraints, and internal inequality constraints
3230 */
3331 template <typename ModelT>
3432 struct D2LDx2
3533 {
36- using ScalarT = typename ModelT::ScalarT;
37- using IdxT = typename ModelT::IdxT;
38- using RealT = typename ModelT::RealT;
39-
40- static_assert (std::is_same_v<ScalarT, double > && std::is_same_v<IdxT, size_t >,
41- " D2LDx2 supports double values and size_t indices" );
42-
4334 /* *
4435 * @param[in] model - Pointer to the model to be differentiated
4536 * @param[in] n_var - Number of local variables
4637 * @param[in] variable_indices - Global index of each local variable
4738 * @param[in] x - Local variables
4839 * @param[in] sigma - Objective factor
49- * @param[in] lambda - Local constraint multipliers
50- * @param[in,out] entries - Hessian entries
40+ * @param[in] lambda_c - Local equality constraint multipliers
41+ * @param[in] lambda_d - Local inequality constraint multipliers
42+ * @param[out] rows - Row indices
43+ * @param[out] cols - Column indices
44+ * @param[out] vals - Values
45+ * @param[out] nnz - Number of nonzeros
5146 */
47+ template <typename ScalarT, typename IdxT, typename RealT>
5248 static void eval (const ModelT* model,
5349 const size_t n_var,
5450 const IdxT* variable_indices,
5551 const ScalarT* x,
5652 const RealT sigma,
57- const RealT* lambda,
58- CooEntries* entries)
53+ const RealT* lambda_c,
54+ const RealT* lambda_d,
55+ IdxT* rows,
56+ IdxT* cols,
57+ ScalarT* vals,
58+ IdxT& nnz)
5959 {
60- std::array<ScalarT, ModelT:: VARIABLE_SIZE > gradient{};
60+ std::array<ScalarT, ObjectiveT:: NUM_INPUTS > gradient{};
6161 for (size_t var_i = 0 ; var_i < n_var; ++var_i)
6262 {
6363 // Sparse storage. @see LowerSparseStorage.hpp
6464 ScalarT* seed = __enzyme_todense<ScalarT*>((void *) ident_load<ScalarT, IdxT>,
6565 (void *) ident_store<ScalarT, IdxT>,
6666 var_i);
67- ScalarT* dgradient = __enzyme_todense<ScalarT*>((void *) mapped_load ,
68- (void *) lower_store ,
67+ ScalarT* dgradient = __enzyme_todense<ScalarT*>((void *) sparse_load<ScalarT, IdxT> ,
68+ (void *) sparse_store<ScalarT, IdxT> ,
6969 var_i,
70+ 1.0 , // value scaling
7071 variable_indices,
7172 variable_indices,
72- entries);
73+ rows,
74+ cols,
75+ vals,
76+ &nnz);
7377
74- __enzyme_fwddiff<void >((void *) lagrangianGradient,
78+ __enzyme_fwddiff<void >((void *) lagrangianGradient<ScalarT, RealT> ,
7579 enzyme_const,
7680 model,
7781 enzyme_dup,
@@ -83,23 +87,38 @@ namespace GridKit
8387 enzyme_const,
8488 sigma,
8589 enzyme_const,
86- lambda);
90+ lambda_c,
91+ enzyme_const,
92+ lambda_d);
8793 }
8894 }
8995
9096 private:
91- __attribute__ ((always_inline)) inline static ScalarT lagrangian (const ModelT* model,
92- const ScalarT* x,
93- const RealT sigma,
94- const RealT* lambda)
97+ using ObjectiveT = ModelWrapper<ModelT, MemberFunctions::Objective>;
98+ using EqualitiesT = ModelWrapper<ModelT, MemberFunctions::Equalities>;
99+ using InequalitiesT = ModelWrapper<ModelT, MemberFunctions::Inequalities>;
100+
101+ template <typename ScalarT, typename RealT>
102+ [[gnu::always_inline]] inline static ScalarT lagrangian (const ModelT* model,
103+ const ScalarT* x,
104+ const RealT sigma,
105+ const RealT* lambda_c,
106+ const RealT* lambda_d)
95107 {
96- std::array<ScalarT, ModelT:: CONSTRAINT_SIZE > g {};
97- ModelWrapper<ModelT, MemberFunctions::Constraints>:: eval (model, x, g .data ());
108+ std::array<ScalarT, EqualitiesT:: NUM_OUTPUTS > c {};
109+ EqualitiesT:: eval (model, x, c .data ());
98110
99- ScalarT value = sigma * ModelWrapper<ModelT, MemberFunctions::Objective>::eval (model, x);
100- for (size_t i = 0 ; i < ModelT::CONSTRAINT_SIZE ; ++i)
111+ std::array<ScalarT, InequalitiesT::NUM_OUTPUTS > d{};
112+ InequalitiesT::eval (model, x, d.data ());
113+
114+ ScalarT value = sigma * ObjectiveT::eval (model, x);
115+ for (size_t i = 0 ; i < EqualitiesT::NUM_OUTPUTS ; ++i)
116+ {
117+ value += lambda_c[i] * c[i];
118+ }
119+ for (size_t i = 0 ; i < InequalitiesT::NUM_OUTPUTS ; ++i)
101120 {
102- value += lambda [i] * g [i];
121+ value += lambda_d [i] * d [i];
103122 }
104123 return value;
105124 }
@@ -109,13 +128,15 @@ namespace GridKit
109128 *
110129 * Only the tangent of `gradient` is used, so it is not reset.
111130 */
112- __attribute__ ((always_inline)) inline static void lagrangianGradient (const ModelT* model,
113- const ScalarT* x,
114- ScalarT* gradient,
115- const RealT sigma,
116- const RealT* lambda)
131+ template <typename ScalarT, typename RealT>
132+ [[gnu::always_inline]] inline static void lagrangianGradient (const ModelT* model,
133+ const ScalarT* x,
134+ ScalarT* gradient,
135+ const RealT sigma,
136+ const RealT* lambda_c,
137+ const RealT* lambda_d)
117138 {
118- __enzyme_autodiff<void >((void *) lagrangian,
139+ __enzyme_autodiff<void >((void *) lagrangian<ScalarT, RealT> ,
119140 enzyme_const,
120141 model,
121142 enzyme_dup,
@@ -124,7 +145,9 @@ namespace GridKit
124145 enzyme_const,
125146 sigma,
126147 enzyme_const,
127- lambda);
148+ lambda_c,
149+ enzyme_const,
150+ lambda_d);
128151 }
129152 };
130153 } // namespace Sparse
0 commit comments