sparsematrix.h

Go to the documentation of this file.

#ifndef ADOLC_SPARSE_MATRIX_H
#define ADOLC_SPARSE_MATRIX_H
 
#include <adolc/drivers/absnormalformconcept.h>
#include <adolc/internal/common.h>
#include <array>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <utility>
#include <vector>
 
namespace ADOLC::Sparse {
 
class CoordinateFormatTripled {
public:
  using coordinate_type = unsigned int;
  struct Coordinates {
    coordinate_type rowIndex_{0};
    coordinate_type colIndex_{0};
  };
 
private:
  Coordinates coords_;
  double value_{0.0};
 
public:
  CoordinateFormatTripled() = default;
  ~CoordinateFormatTripled() = default;
 
  CoordinateFormatTripled(const CoordinateFormatTripled &) = default;
  CoordinateFormatTripled &operator=(const CoordinateFormatTripled &) = default;
 
  CoordinateFormatTripled(CoordinateFormatTripled &&) noexcept = default;
  CoordinateFormatTripled &
  operator=(CoordinateFormatTripled &&) noexcept = default;
 
  constexpr CoordinateFormatTripled(Coordinates coords, double value)
      : coords_{coords}, value_(value) {}
  constexpr CoordinateFormatTripled(coordinate_type rowIndex,
                                    coordinate_type colIndex, double value)
      : coords_{Coordinates{rowIndex, colIndex}}, value_(value) {}
 
  coordinate_type rowIndex() const noexcept { return coords_.rowIndex_; }
  coordinate_type colIndex() const noexcept { return coords_.colIndex_; }
  double value() const noexcept { return value_; }
};
 
class SparseMatrix {
public:
  using value_type = CoordinateFormatTripled;
  using storage_type = std::vector<value_type>;
 
private:
  storage_type entries_;
 
public:
  ~SparseMatrix() = default;
  SparseMatrix() = default;
  explicit SparseMatrix(storage_type &&entries) noexcept
      : entries_(std::move(entries)) {}
 
  template <size_t N>
  explicit SparseMatrix(
      const std::array<CoordinateFormatTripled, N> &entries) noexcept
      : entries_(entries.begin(), entries.end()) {}
  explicit SparseMatrix(size_t size) : entries_(size) {}
 
  SparseMatrix(const SparseMatrix &other) = default;
  SparseMatrix &operator=(const SparseMatrix &other) = default;
 
  SparseMatrix(SparseMatrix &&other) noexcept = default;
  SparseMatrix &operator=(SparseMatrix &&other) noexcept = default;
 
  storage_type &entries() { return entries_; }
  const storage_type &entries() const { return entries_; }
 
  size_t size() const { return entries_.size(); }
  bool empty() const { return entries_.empty(); }
 
  void resize(size_t size) { entries_.resize(size); }
  void clear() { entries_.clear(); }
 
  void push_back(CoordinateFormatTripled value) { entries_.push_back(value); }
 
  unsigned int rowIndex(size_t idx) const { return entries_[idx].rowIndex(); }
  unsigned int colIndex(size_t idx) const { return entries_[idx].colIndex(); }
  double value(size_t idx) const { return entries_[idx].value(); }
 
  value_type &operator[](size_t idx) { return entries_[idx]; }
  const value_type &operator[](size_t idx) const { return entries_[idx]; }
};
 
namespace detail {
 
enum class SparseANFBlock { Y, J, Z, L };
 
inline SparseANFBlock
classifySparseANFBlock(CoordinateFormatTripled::Coordinates coords, int depen,
                       int indep) {
  const bool isDependentRow =
      coords.rowIndex_ < static_cast<unsigned int>(depen);
  const bool isIndependentCol =
      coords.colIndex_ < static_cast<unsigned int>(indep);
 
  if (isDependentRow && isIndependentCol)
    return SparseANFBlock::Y;
  else if (isDependentRow)
    return SparseANFBlock::J;
  else if (isIndependentCol)
    return SparseANFBlock::Z;
  else
    return SparseANFBlock::L;
}
} // namespace detail
 
struct SparseShape : ANFShape {
  size_t y{0};
  size_t j{0};
  size_t z{0};
  size_t l{0};
 
  SparseShape() = default;
  SparseShape(size_t y, size_t j, size_t z, size_t l)
      : y(y), j(j), z(z), l(l) {};
};
 
struct SparseANF {
  using Shape = SparseShape;
 
  Shape shape{};
 
  SparseMatrix Y;
  SparseMatrix J;
  SparseMatrix Z;
  SparseMatrix L;
 
  std::vector<double> cz;
  std::vector<double> cy;
  std::vector<double> y;
  std::vector<double> z;
 
  ~SparseANF() = default;
  SparseANF() = default;
 
  explicit SparseANF(Shape dims) { resize(dims); }
 
  SparseANF(SparseANF &&) = default;
  SparseANF &operator=(SparseANF &&) = default;
 
  SparseANF(const SparseANF &) = default;
  SparseANF &operator=(const SparseANF &) = default;
 
  bool empty() const {
    return shape.j == 0 && shape.l == 0 && shape.y == 0 && shape.z == 0;
  }
 
  void clear() {
    shape = {};
    Y.clear();
    J.clear();
    Z.clear();
    L.clear();
    cy.clear();
    cz.clear();
  }
 
  Shape dims() const { return shape; }
 
  void updateCy() {
    cy = y;
    for (size_t i = 0; i < J.size(); ++i) {
      cy[J.rowIndex(i)] -= J.value(i) * std::fabs(z[J.colIndex(i)]);
    }
  }
 
  void updateCz() {
    cz = z;
    for (size_t i = 0; i < L.size(); ++i) {
      cz[L.rowIndex(i)] -= L.value(i) * std::fabs(z[L.colIndex(i)]);
    }
  }
 
  void resize(Shape dims) {
    shape = dims;
    Y.resize(shape.y);
    J.resize(shape.j);
    Z.resize(shape.z);
    L.resize(shape.l);
  }
};
 
} // namespace ADOLC::Sparse
 
#endif // ADOLC_SPARSE_MATRIX_H