use new EllAlgo 1.3.6

Author: luk036

CMakeLists.txt

@@ -12,7 +12,7 @@ option(INSTALL_ONLY "Enable for installation only" OFF)
 # Note: update this to your new project's name and version
 project(
   LmiSolver
-  VERSION 1.3.4
+  VERSION 1.3.5
   LANGUAGES CXX
 )
 

bench/BM_lmi.cpp

@@ -6,9 +6,11 @@
 #include <tuple>                        // for tuple
 #include <type_traits>                  // for move
 #include <vector>                       // for vector
-#include <xtensor/xlayout.hpp>          // for layout_type, layout_type::ro...
+#include <xtensor/xarray.hpp>           // for xarray_container
+#include <xtensor/xlayout.hpp>          // for layout_type, layout_type::row...
 #include <xtensor/xmath.hpp>            // for sum
-#include <xtensor/xoperation.hpp>       // for operator*
+#include <xtensor/xoperation.hpp>       // for xfunction_type_t, operator*
+#include <xtensor/xreducer.hpp>         // for xreducer
 #include <xtensor/xtensor_forward.hpp>  // for xarray
 
 #include "benchmark/benchmark.h" // for BENCHMARK, State, BENCHMARK_...
@@ -20,6 +22,7 @@
  */
 template <typename Oracle> class my_oracle {
   using Arr = xt::xarray<double, xt::layout_type::row_major>;
+  using Vec = std::valarray<double>;
   using Cut = std::pair<Arr, double>;
 
 private:
@@ -138,36 +141,37 @@ static void LMI_old(benchmark::State &state) {
 }
 BENCHMARK(LMI_old);
 
-/**
- * @brief
- *
- * @param[in,out] state
- */
-static void LMI_No_Trick(benchmark::State &state) {
-  using Arr = xt::xarray<double, xt::layout_type::row_major>;
-
-  // const auto c = Arr {1.0, -1.0, 1.0};
-  const auto F1 = std::vector<Arr>{{{-7.0, -11.0}, {-11.0, 3.0}},
-                                   {{7.0, -18.0}, {-18.0, 8.0}},
-                                   {{-2.0, -8.0}, {-8.0, 1.0}}};
-  const auto B1 = Arr{{33.0, -9.0}, {-9.0, 26.0}};
-  const auto F2 = std::vector<Arr>{
-      {{-21.0, -11.0, 0.0}, {-11.0, 10.0, 8.0}, {0.0, 8.0, 5.0}},
-      {{0.0, 10.0, 16.0}, {10.0, -10.0, -10.0}, {16.0, -10.0, 3.0}},
-      {{-5.0, 2.0, -17.0}, {2.0, -6.0, 8.0}, {-17.0, 8.0, 6.0}}};
-  const auto B2 = Arr{{14.0, 9.0, 40.0}, {9.0, 91.0, 10.0}, {40.0, 10.0, 15.0}};
-
-  while (state.KeepRunning()) {
-    auto P = my_oracle<LmiOracle<Arr>>(F1, B1, F2, B2, Arr{1.0, -1.0, 1.0});
-    auto E = Ell<Arr>(10.0, Arr{0.0, 0.0, 0.0});
-    E.no_defer_trick = true;
-    auto t = 1e100; // std::numeric_limits<double>::max()
-    [[maybe_unused]] const auto rslt = cutting_plane_optim(P, E, t);
-  }
-}
-
-// Register the function as a benchmark
-BENCHMARK(LMI_No_Trick);
+// /**
+//  * @brief
+//  *
+//  * @param[in,out] state
+//  */
+// static void LMI_No_Trick(benchmark::State &state) {
+//   using Arr = xt::xarray<double, xt::layout_type::row_major>;
+//
+//   // const auto c = Arr {1.0, -1.0, 1.0};
+//   const auto F1 = std::vector<Arr>{{{-7.0, -11.0}, {-11.0, 3.0}},
+//                                    {{7.0, -18.0}, {-18.0, 8.0}},
+//                                    {{-2.0, -8.0}, {-8.0, 1.0}}};
+//   const auto B1 = Arr{{33.0, -9.0}, {-9.0, 26.0}};
+//   const auto F2 = std::vector<Arr>{
+//       {{-21.0, -11.0, 0.0}, {-11.0, 10.0, 8.0}, {0.0, 8.0, 5.0}},
+//       {{0.0, 10.0, 16.0}, {10.0, -10.0, -10.0}, {16.0, -10.0, 3.0}},
+//       {{-5.0, 2.0, -17.0}, {2.0, -6.0, 8.0}, {-17.0, 8.0, 6.0}}};
+//   const auto B2 = Arr{{14.0, 9.0, 40.0}, {9.0, 91.0, 10.0},
+//   {40.0, 10.0, 15.0}};
+//
+//   while (state.KeepRunning()) {
+//     auto P = my_oracle<LmiOracle<Arr>>(F1, B1, F2, B2, Arr{1.0, -1.0, 1.0});
+//     auto E = Ell<Arr>(10.0, Arr{0.0, 0.0, 0.0});
+//     E.no_defer_trick = true;
+//     auto t = 1e100; // std::numeric_limits<double>::max()
+//     [[maybe_unused]] const auto rslt = cutting_plane_optim(P, E, t);
+//   }
+// }
+//
+// // Register the function as a benchmark
+// BENCHMARK(LMI_No_Trick);
 
 BENCHMARK_MAIN();
 

include/lmisolver/ldlt_ext.hpp

@@ -1,6 +1,7 @@
 // -*- coding: utf-8 -*-
 #pragma once
 
+#include <cassert> // for assert
 #include <cstddef> // for size_t
 #include <ellalgo/ell_matrix.hpp>
 #include <utility> // for pair
@@ -16,16 +17,15 @@
  */
 template <typename Arr036> class ldlt_ext {
   using Vec = std::valarray<double>;
-  using Mat = Matrix;
   using Rng = std::pair<size_t, size_t>;
 
 public:
   Rng p{0U, 0U};   //!< the rows where the process starts and stops
   Vec witness_vec; //!< witness vector
+  const size_t _n; //!< dimension
 
 private:
-  const size_t _n; //!< dimension
-  Mat T;           //!< temporary storage
+  Matrix T; //!< temporary storage
 
   // static Vec zeros_vec(size_t n);
   // static Mat zeros_mat(size_t n);
@@ -36,7 +36,7 @@ template <typename Arr036> class ldlt_ext {
    *
    * @param[in] N dimension
    */
-  explicit ldlt_ext(size_t N);
+  explicit ldlt_ext(size_t N) : witness_vec(0.0, N), _n{N}, T{N} {}
 
   ldlt_ext(const ldlt_ext &) = delete;
   ldlt_ext &operator=(const ldlt_ext &) = delete;
@@ -52,7 +52,21 @@ 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 Arr036 &A) -> bool {
+  // auto factorize(const Arr036 &A) -> bool {
+  //   return this->factor([&](size_t i, size_t j) { return A(i, j); });
+  // }
+
+  /**
+   * @brief Perform LDLT Factorization
+   *
+   * @param[in] A Symmetric Matrix
+   *
+   * If $A$ is positive definite, then $p$ is zero.
+   * If it is not, then $p$ is a positive integer,
+   * such that $v = R^-1 e_p$ is a certificate vector
+   * to make $v'*A[:p,:p]*v < 0$
+   */
+  template <typename Mat> auto factorize(const Mat &A) -> bool {
     return this->factor([&](size_t i, size_t j) { return A(i, j); });
   }
 
@@ -132,7 +146,35 @@ template <typename Arr036> class ldlt_ext {
    * @param[in] A
    * @return double
    */
-  [[nodiscard]] auto sym_quad(const Arr036 &A) const -> double;
+  template <typename Mat> auto sym_quad(const Mat &A) const -> double {
+    auto res = double{};
+    const auto &v = this->witness_vec;
+    // const auto& [start, stop] = this->p;
+    const auto &start = this->p.first;
+    const auto &stop = this->p.second;
+    for (auto i = start; i != stop; ++i) {
+      auto s = double{};
+      for (auto j = i + 1; j != stop; ++j) {
+        s += A(i, j) * v[j];
+      }
+      res += v[i] * (A(i, i) * v[i] + 2 * s);
+    }
+    return res;
+  }
+
+  /**
+   * @brief Return upper triangular matrix $R$ where $A = R^T R$
+   *
+   * @return typename ldlt_ext<Arr036>::Mat
+   */
+  template <typename Mat> auto sqrt(Mat &M) -> void {
+    assert(this->is_spd());
 
-  auto sqrt() -> Arr036;
+    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) {
+        M(i, j) = this->T(j, i) * M(i, i);
+      }
+    }
+  }
 };

source/ldlt_ext.cpp

@@ -3,9 +3,6 @@
 #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)
@@ -29,28 +26,7 @@ template <typename Arr036> auto ldlt_ext<Arr036>::witness() -> double {
   return -this->T(m, m);
 }
 
-/**
- * @brief Calculate v'*{A}(p,p)*v
- *
- * @param[in] A
- * @return double
- */
-template <typename Arr036>
-auto ldlt_ext<Arr036>::sym_quad(const Arr036 &A) const -> double {
-  auto res = double{};
-  const auto &v = this->witness_vec;
-  // const auto& [start, stop] = this->p;
-  const auto &start = this->p.first;
-  const auto &stop = this->p.second;
-  for (auto i = start; i != stop; ++i) {
-    auto s = double{};
-    for (auto j = i + 1; j != stop; ++j) {
-      s += A(i, j) * v[j];
-    }
-    res += v[i] * (A(i, i) * v[i] + 2 * s);
-  }
-  return res;
-}
+template auto ldlt_ext<std::valarray<double>>::witness() -> double;
 
 #include <xtensor/xarray.hpp>          // for xarray
 #include <xtensor/xcontainer.hpp>      // for xcontainer
@@ -59,22 +35,4 @@ auto ldlt_ext<Arr036>::sym_quad(const Arr036 &A) const -> double {
 
 using Arr = xt::xarray<double, xt::layout_type::row_major>;
 
-/**
- * @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() -> Arr036 {
-  assert(this->is_spd());
-
-  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) {
-      M(i, j) = this->T(j, i) * M(i, i);
-    }
-  }
-  return M;
-}
-
-template class ldlt_ext<Arr>;
+template auto ldlt_ext<Arr>::witness() -> double;

specific.cmake

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

test/source/test_ldlt.cpp

@@ -1,38 +1,49 @@
 #include <doctest/doctest.h> // for ResultBuilder, TestCase, CHECK
 
-#include <lmisolver/ldlt_ext.hpp>      // for ldlt_ext
-#include <xtensor/xarray.hpp>          // for xarray_container
-#include <xtensor/xcontainer.hpp>      // for xcontainer<>::inner_shape_type
-#include <xtensor/xlayout.hpp>         // for layout_type, layout_type::row...
-#include <xtensor/xtensor_forward.hpp> // for xarray
+#include <lmisolver/ldlt_ext.hpp> // for ldlt_ext
+// #include <xtensor/xarray.hpp>          // for xarray_container
+// #include <xtensor/xcontainer.hpp>      // for xcontainer<>::inner_shape_type
+// #include <xtensor/xlayout.hpp>         // for layout_type,
+// layout_type::row... #include <xtensor/xtensor_forward.hpp> // for xarray
 // #include <xtensor/xarray.hpp>
+#include <ellalgo/ell_matrix.hpp>
 
-using Arr = xt::xarray<double, xt::layout_type::row_major>;
+// using Arr = xt::xarray<double, xt::layout_type::row_major>;
+using Vec = std::valarray<double>;
 
 TEST_CASE("Cholesky test 1") {
-  const auto m1 = Arr{{25.0, 15.0, -5.0}, {15.0, 18.0, 0.0}, {-5.0, 0.0, 11.0}};
-  auto Q1 = ldlt_ext<Arr>(m1.shape()[0]);
+  Matrix m1(3U);
+  m1.row(0) = Vec{25.0, 15.0, -5.0};
+  m1.row(1) = Vec{15.0, 18.0, 0.0};
+  m1.row(2) = Vec{-5.0, 0.0, 11.0};
+  auto Q1 = ldlt_ext<Vec>(3);
   CHECK(Q1.factorize(m1));
 }
 
 TEST_CASE("Cholesky test 2") {
-  const auto m2 = Arr{{18.0, 22.0, 54.0, 42.0},
-                      {22.0, -70.0, 86.0, 62.0},
-                      {54.0, 86.0, -174.0, 134.0},
-                      {42.0, 62.0, 134.0, -106.0}};
-  auto Q2 = ldlt_ext<Arr>(m2.shape()[0]);
+  Matrix m2(4U);
+  m2.row(0) = Vec{18.0, 22.0, 54.0, 42.0};
+  m2.row(1) = Vec{22.0, -70.0, 86.0, 62.0};
+  m2.row(2) = Vec{54.0, 86.0, -174.0, 134.0};
+  m2.row(3) = Vec{42.0, 62.0, 134.0, -106.0};
+
+  auto Q2 = ldlt_ext<Vec>(4);
   Q2.factorize(m2);
   CHECK(!Q2.is_spd());
   // CHECK(Q2.p.second == 2);
 }
 
 TEST_CASE("Cholesky test 3") {
-  const auto m3 = Arr{{0.0, 15.0, -5.0}, {15.0, 18.0, 0.0}, {-5.0, 0.0, 11.0}};
-  auto Q3 = ldlt_ext<Arr>(m3.shape()[0]);
+  Matrix m3(3U);
+  m3.row(0) = Vec{0.0, 15.0, -5.0};
+  m3.row(1) = Vec{15.0, 18.0, 0.0};
+  m3.row(2) = Vec{-5.0, 0.0, 11.0};
+
+  auto Q3 = ldlt_ext<Vec>(3);
   Q3.factorize(m3);
   CHECK(!Q3.is_spd());
   const auto ep3 = Q3.witness();
-  Arr v = xt::zeros<double>({m3.shape()[0]});
+  Vec v(0.0, 3);
   Q3.set_witness_vec(v);
 
   // CHECK(Q3.p.second == 1);