internally use std::valarray

Author: luk036

core

@@ -3,6 +3,7 @@
 
 #include <cstddef> // for size_t
 #include <utility> // for pair
+#include <ellalgo/ell_matrix.hpp>
 
 /**
  * @brief LDLT factorization for LMI
@@ -14,8 +15,8 @@
  *  - O(p^2) per iteration, independent of N
  */
 template <typename Arr036> class ldlt_ext {
-  using Vec = Arr036;
-  using Mat = Arr036;
+  using Vec = std::valarray<double>;
+  using Mat = Matrix;
   using Rng = std::pair<size_t, size_t>;
 
 public:
@@ -51,7 +52,7 @@ template <typename Arr036> class ldlt_ext {
    * such that $v = R^-1 e_p$ is a certificate vector
    * to make $v'*A[:p,:p]*v < 0$
    */
-  auto factorize(const Mat &A) -> bool {
+  auto factorize(const Arr036 &A) -> bool {
     return this->factor([&](size_t i, size_t j) { return A(i, j); });
   }
 
@@ -83,17 +84,17 @@ template <typename Arr036> class ldlt_ext {
       this->T(i, i) = d;
 
       if constexpr (Allow_semidefinite) {
-        if (d < 0.) {
+        if (d < 0.0) {
           // this->stop = i + 1;
           stop = i + 1;
           break;
         }
-        if (d == 0.) {
+        if (d == 0.0) {
           start = i + 1;
           // restart at i + 1, special as an LMI oracle
         }
       } else { // not Allow_semidefinite
-        if (d <= 0.) {
+        if (d <= 0.0) {
           stop = i + 1;
           break;
         }
@@ -125,7 +126,7 @@ template <typename Arr036> class ldlt_ext {
    * @param[in] A
    * @return double
    */
-  [[nodiscard]] auto sym_quad(const Vec &A) const -> double;
+  [[nodiscard]] auto sym_quad(const Arr036 &A) const -> double;
 
-  auto sqrt() -> Mat;
+  auto sqrt() -> Arr036;
 };

include/lmisolver/ldlt_ext.hpp

@@ -3,6 +3,11 @@
 #include <cmath>   // for sqrt
 #include <lmisolver/ldlt_ext.hpp>
 
+template <typename Arr036>
+ldlt_ext<Arr036>::ldlt_ext(size_t N)
+    : witness_vec(0.0, N), _n{N}, T{N} {
+}
+
 /**
  * @brief witness that certifies $A$ is not
  * symmetric positive definite (spd)
@@ -16,11 +21,11 @@ template <typename Arr036> auto ldlt_ext<Arr036>::witness() -> double {
   const auto &start = this->p.first;
   const auto &n = this->p.second;
   auto m = n - 1; // assume stop > 0
-  this->witness_vec(m) = 1.;
+  this->witness_vec[m] = 1.0;
   for (auto i = m; i > start; --i) {
-    this->witness_vec(i - 1) = 0.;
+    this->witness_vec[i - 1] = 0.0;
     for (auto k = i; k != n; ++k) {
-      this->witness_vec(i - 1) -= this->T(k, i - 1) * this->witness_vec(k);
+      this->witness_vec[i - 1] -= this->T(k, i - 1) * this->witness_vec[k];
     }
   }
   return -this->T(m, m);
@@ -33,7 +38,7 @@ template <typename Arr036> auto ldlt_ext<Arr036>::witness() -> double {
  * @return double
  */
 template <typename Arr036>
-auto ldlt_ext<Arr036>::sym_quad(const typename ldlt_ext<Arr036>::Vec &A) const
+auto ldlt_ext<Arr036>::sym_quad(const Arr036 &A) const
     -> double {
   auto res = double{};
   const auto &v = this->witness_vec;
@@ -43,9 +48,9 @@ auto ldlt_ext<Arr036>::sym_quad(const typename ldlt_ext<Arr036>::Vec &A) const
   for (auto i = start; i != stop; ++i) {
     auto s = double{};
     for (auto j = i + 1; j != stop; ++j) {
-      s += A(i, j) * v(j);
+      s += A(i, j) * v[j];
     }
-    res += v(i) * (A(i, i) * v(i) + 2 * s);
+    res += v[i] * (A(i, i) * v[i] + 2 * s);
   }
   return res;
 }
@@ -57,21 +62,16 @@ auto ldlt_ext<Arr036>::sym_quad(const typename ldlt_ext<Arr036>::Vec &A) const
 
 using Arr = xt::xarray<double, xt::layout_type::row_major>;
 
-template <typename Arr036>
-ldlt_ext<Arr036>::ldlt_ext(size_t N)
-    : witness_vec{xt::zeros<double>({N})}, _n{N}, T{xt::zeros<double>({N, N})} {
-}
-
 /**
  * @brief Return upper triangular matrix $R$ where $A = R^T R$
  *
  * @return typename ldlt_ext<Arr036>::Mat
  */
 template <typename Arr036>
-auto ldlt_ext<Arr036>::sqrt() -> typename ldlt_ext<Arr036>::Mat {
+auto ldlt_ext<Arr036>::sqrt() -> Arr036 {
   assert(this->is_spd());
 
-  ldlt_ext<Arr036>::Mat M = xt::zeros<double>({this->_n, this->_n});
+  Arr036 M = xt::zeros<double>({this->_n, this->_n});
   for (auto i = 0U; i != this->_n; ++i) {
     M(i, i) = std::sqrt(this->T(i, i));
     for (auto j = i + 1; j != this->_n; ++j) {
@@ -81,4 +81,4 @@ auto ldlt_ext<Arr036>::sqrt() -> typename ldlt_ext<Arr036>::Mat {
   return M;
 }
 
-template class ldlt_ext<Arr>;
+template class ldlt_ext<Arr>;

source/ldlt_ext.cpp

@@ -48,11 +48,12 @@ auto LmiOracle<Arr036>::assess_feas(const Arr036 &x)
   }
   auto ep = this->_Q.witness();
   auto g = Arr036{xt::zeros<double>({x.size()})};
+  // auto g = Arr036(x.size(), 0.0);
   for (auto i = 0U; i != n; ++i) {
     g(i) = this->_Q.sym_quad(this->_F[i]);
   }
   return {{std::move(g), ep}};
 }
 
 using Arr = xt::xarray<double, xt::layout_type::row_major>;
-template class LmiOracle<Arr>;
+template class LmiOracle<Arr>;

source/lmi_oracle.cpp

@@ -19,7 +19,7 @@ endif(xtensor_ADDED)
 
 CPMAddPackage(
   NAME EllAlgo
-  GIT_TAG 1.3.2
+  GIT_TAG 1.3.3
   GITHUB_REPOSITORY luk036/ellalgo-cpp
   OPTIONS "INSTALL_ONLY YES" # create an installable target
 )