dmatrix.h

#ifndef _DMATRIX_H_
#define _DMATRIX_H_
#include <iostream>
#include <sstream>
#include <string>
#include <utility>
#include <vector>
/* Either text or binary can be specified here, but we use binary
 * to emphasize this class is normally serialized binary for
 * speed*/
#include "mspass/utility/MsPASSError.h"
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/serialization/vector.hpp>
namespace mspass {
namespace utility {
//==================================================================
class dmatrix_index_error : public MsPASSError {
public:
  dmatrix_index_error(const size_t nrmax, const size_t ncmax, const size_t ir,
                      const size_t ic) {
    row = ir;
    column = ic;
    nrr = nrmax;
    ncc = ncmax;
    std::ostringstream oss;
    oss << "dmatrix object:  indexing error" << std::endl
        << "Matrix index (" << row << "," << column
        << ")is outside range = " << nrr << "," << ncc << std::endl;
    message = oss.str();
    badness = ErrorSeverity::Invalid;
  };
  /* necessary baggage for some compilers - empty destructor */
  ~dmatrix_index_error() throw() {};
 
private:
  size_t row, column;
  size_t nrr, ncc;
};
class dmatrix_size_error : public MsPASSError {
public:
  dmatrix_size_error(const size_t nr1, const size_t nc1, const size_t nr2,
                     const size_t nc2) {
    nrow1 = nr1;
    ncol1 = nc1;
    nrow2 = nr2;
    ncol2 = nc2;
    std::ostringstream oss;
    oss << "dmatrix class:   size mismatch error in binary operator"
        << std::endl
        << "matrix on left is " << nrow1 << "X" << ncol1
        << " while matrix on right is " << nrow2 << "X" << ncol2 << std::endl;
    message = oss.str();
    badness = ErrorSeverity::Invalid;
  };
  /* necessary baggage for some compilers - empty destructor */
  ~dmatrix_size_error() throw() {};
 
private:
  size_t nrow1, ncol1, nrow2, ncol2;
};
class dmatrix {
public:
  dmatrix();
  dmatrix(const size_t nr, const size_t nc);
  dmatrix(const dmatrix &other);
  ~dmatrix();
  double operator()(const size_t rowindex, const size_t colindex) const;
  double &operator()(size_t r, size_t c);
  dmatrix &operator=(const dmatrix &other);
  void swap(dmatrix &other) noexcept {
    ary.swap(other.ary);
    std::swap(length, other.length);
    std::swap(nrr, other.nrr);
    std::swap(ncc, other.ncc);
  }
  bool storage_is_consistent() const noexcept {
    if (ary.size() != length)
      return false;
    if (nrr == 0 || ncc == 0)
      return nrr == 0 && ncc == 0 && length == 0;
    return length / nrr == ncc && length % nrr == 0;
  }
  dmatrix &operator+=(const dmatrix &other);
  dmatrix &operator-=(const dmatrix &other);
  dmatrix operator+(const dmatrix &other) const;
  dmatrix operator-(const dmatrix &other) const;
  // friend class dvector;
  friend dmatrix operator*(const dmatrix &A, const dmatrix &B);
  // Scale a matrix by a constant.  X=c*A where c is a constant.
 
  friend dmatrix operator*(const double &s, const dmatrix &A) noexcept;
  dmatrix operator*(double s) const noexcept;
  friend dmatrix tr(const dmatrix &A) noexcept;
  double *get_address(size_t r, size_t c) const;
  friend std::ostream &operator<<(std::ostream &os, dmatrix &A);
  size_t rows() const;
  size_t columns() const;
  std::vector<size_t> size() const;
  void zero();
 
protected:
  std::vector<double> ary; // initial size of container 0
  size_t length;
  size_t nrr;
  size_t ncc;
 
private:
  friend class boost::serialization::access;
  template <class Archive>
  void serialize(Archive &ar, const unsigned int version) {
    ar & nrr & ncc & length;
    ar & ary;
  }
};
class dvector : public dmatrix {
public:
  dvector() : dmatrix() {};
  dvector(size_t nrv) : dmatrix(nrv, 1) {};
  dvector(const dvector &other);
  dvector &operator=(const dvector &other);
  double &operator()(size_t rowindex);
  friend dvector operator*(const dmatrix &A, const dvector &x);
};
 
} // namespace utility
} // end namespace mspass
#endif