block-kind.hpp

Go to the documentation of this file.

#pragma once
 
#include "proxsuite-nlp/config.hpp"
 
#include <Eigen/Core>
#include <numeric>
#include <vector>
 
namespace proxsuite {
namespace nlp {
namespace linalg {
 
using isize = Eigen::Index;
using usize = std::make_unsigned<isize>::type;

enum BlockKind {
  Zero,
  Diag,
  TriL,
  TriU,
  Dense,
};

PROXSUITE_NLP_DLLAPI BlockKind trans(BlockKind a) noexcept;

PROXSUITE_NLP_DLLAPI BlockKind add(BlockKind a, BlockKind b) noexcept;

PROXSUITE_NLP_DLLAPI BlockKind mul(BlockKind a, BlockKind b) noexcept;

struct PROXSUITE_NLP_DLLAPI SymbolicBlockMatrix {
  BlockKind *m_data;
  isize *segment_lens;
  isize segments_count;
  isize outer_stride;
  bool performed_llt = false;
 
private:
  void alloc() {
    m_data = new BlockKind[usize(size())];
    segment_lens = new isize[usize(nsegments())];
  }
 
public:
  SymbolicBlockMatrix() {}
  SymbolicBlockMatrix(isize nc, isize os)
      : segments_count(nc), outer_stride(os) {
    alloc();
  }

  SymbolicBlockMatrix(BlockKind *data, isize *lens, isize nc, isize os)
      : m_data(data), segment_lens(lens), segments_count(nc), outer_stride(os) {
  }

  SymbolicBlockMatrix(SymbolicBlockMatrix const &other) = default;
  SymbolicBlockMatrix &operator=(SymbolicBlockMatrix const &other) = default;
  SymbolicBlockMatrix(const std::vector<isize> &lens)
      : segments_count(static_cast<isize>(lens.size())) {
    outer_stride = segments_count;
 
    alloc();
 
    std::fill_n(m_data, size(), BlockKind::Dense);
    std::copy_n(lens.begin(), nsegments(), segment_lens);
  }

  SymbolicBlockMatrix copy() const;
 
  isize nsegments() const noexcept { return segments_count; }
  isize size() const noexcept { return segments_count * outer_stride; }
  BlockKind *ptr(isize i, isize j) noexcept {
    return m_data + (i + j * outer_stride);
  }
  BlockKind const *ptr(isize i, isize j) const noexcept {
    return m_data + (i + j * outer_stride);
  }

  SymbolicBlockMatrix submatrix(isize i, isize n) noexcept;
  BlockKind &operator()(isize i, isize j) noexcept { return *ptr(i, j); }
  const BlockKind &operator()(isize i, isize j) const noexcept {
    return *ptr(i, j);
  }

  Eigen::ComputationInfo
  brute_force_best_permutation(SymbolicBlockMatrix const &in, isize *best_perm,
                               isize *iwork);
  bool check_if_symmetric() const noexcept;
  isize count_nnz() const noexcept;
  bool llt_in_place() noexcept;
 
  SymbolicBlockMatrix transpose() const {
    const auto &self = *this;
    SymbolicBlockMatrix out(self.copy());
    for (isize i = 0; i < nsegments(); ++i) {
      for (isize j = 0; j < nsegments(); ++j) {
        out(i, j) = trans(self(j, i));
      }
    }
    return out;
  }
};
 
PROXSUITE_NLP_DLLAPI void
print_sparsity_pattern(const SymbolicBlockMatrix &smat) noexcept;

PROXSUITE_NLP_DLLAPI void
symbolic_deep_copy(const SymbolicBlockMatrix &in, SymbolicBlockMatrix &out,
                   isize const *perm = nullptr) noexcept;
 
} // namespace linalg
} // namespace nlp
} // namespace proxsuite