dmatrix.h

#ifndef _DMATRIX_H_
#define _DMATRIX_H_
#include <iostream>
#include <sstream>
#include <string>
#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);
  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;
  /* \brief Get a pointer to the location of a matrix component.
 
  Although a sharp knife it is useful at times to get a raw pointer to
  the data in a dmatrix.   A common one is using the BLAS to do vector
  operations for speed.   Users of this method must note that the data
  for a dmatrix is stored as a single rowXcolumn std::vector container.
  The matrix is stored in Fortran order (column1, column2, ...).
  The contiguous memory guarantee of std::vector allows vector operations
  with the BLAS to work by rows or columns.
 
  \param r is the row index of the desired address
  \param c is the column index of the desired memory address.
  \return pointer to component at row r and column c.
  \exception dmatrix_size_error will be throw if r or c are outside
    matrix dimensions.
    */
  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, 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