dmatrix.h

#ifndef _DMATRIX_H_
#define _DMATRIX_H_
#include <string>
#include <iostream>
#include <sstream>
#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 <boost/serialization/vector.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/binary_iarchive.hpp>
#include "mspass/utility/MsPASSError.h"
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);
};
 
} // end utility namespace  
} // end namespace mspass
#endif