/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/GeneralMatrix.h,v 1.3 2004/11/24 20:41:59 kamenik Exp $ */

/* Tag $Name:  $ */

#ifndef GENERAL_MATRIX_H
#define GENERAL_MATRIX_H

#include "Vector.h"

#include <algorithm>

class GeneralMatrix;

class ConstGeneralMatrix
{
  friend class GeneralMatrix;
protected:
  ConstVector data;
  int rows;
  int cols;
  int ld;
public:
  ConstGeneralMatrix(const double *d, int m, int n)
    : data(d, m*n), rows(m), cols(n), ld(m)
  {
  }
  ConstGeneralMatrix(const GeneralMatrix &m);
  ConstGeneralMatrix(const GeneralMatrix &m, int i, int j, int nrows, int ncols);
  ConstGeneralMatrix(const ConstGeneralMatrix &m, int i, int j, int nrows, int ncols);
  virtual ~ConstGeneralMatrix()
  {
  }

  const double &
  get(int i, int j) const
  {
    return data[j*ld+i];
  }
  int
  numRows() const
  {
    return rows;
  }
  int
  numCols() const
  {
    return cols;
  }
  int
  getLD() const
  {
    return ld;
  }
  const double *
  base() const
  {
    return data.base();
  }
  const ConstVector &
  getData() const
  {
    return data;
  }

  double getNormInf() const;
  double getNorm1() const;
  /* x = scalar(a)*x + scalar(b)*this*d */
  void multVec(double a, Vector &x, double b, const ConstVector &d) const;
  /* x = scalar(a)*x + scalar(b)*this'*d */
  void multVecTrans(double a, Vector &x, double b, const ConstVector &d) const;
  /* x = x + this*d */
  void
  multaVec(Vector &x, const ConstVector &d) const
  {
    multVec(1.0, x, 1.0, d);
  }
  /* x = x + this'*d */
  void
  multaVecTrans(Vector &x, const ConstVector &d) const
  {
    multVecTrans(1.0, x, 1.0, d);
  }
  /* x = x - this*d */
  void
  multsVec(Vector &x, const ConstVector &d) const
  {
    multVec(1.0, x, -1.0, d);
  }
  /* x = x - this'*d */
  void
  multsVecTrans(Vector &x, const ConstVector &d) const
  {
    multVecTrans(1.0, x, -1.0, d);
  }
  /* m = inv(this)*m */
  void multInvLeft(GeneralMatrix &m) const;
  /* m = inv(this')*m */
  void multInvLeftTrans(GeneralMatrix &m) const;
  /* d = inv(this)*d */
  void multInvLeft(Vector &d) const;
  /* d = inv(this')*d */
  void multInvLeftTrans(Vector &d) const;

  bool isFinite() const;
  /** Returns true of the matrix is exactly zero. */
  bool isZero() const;

  virtual void print() const;
protected:
  void multInvLeft(const char *trans, int mrows, int mcols, int mld, double *d) const;
};

class GeneralMatrix
{
  friend class ConstGeneralMatrix;
protected:
  Vector data;
  int rows;
  int cols;
  int ld;
public:
  GeneralMatrix(int m, int n)
    : data(m*n), rows(m), cols(n), ld(m)
  {
  }
  GeneralMatrix(const double *d, int m, int n)
    : data(d, m*n), rows(m), cols(n), ld(m)
  {
  }
  GeneralMatrix(double *d, int m, int n)
    : data(d, m*n), rows(m), cols(n), ld(m)
  {
  }
  GeneralMatrix(const GeneralMatrix &m);
  GeneralMatrix(const ConstGeneralMatrix &m);
  GeneralMatrix(const GeneralMatrix &m, const char *dummy); // transpose
  GeneralMatrix(const ConstGeneralMatrix &m, const char *dummy); // transpose
  GeneralMatrix(const GeneralMatrix &m, int i, int j, int nrows, int ncols);
  GeneralMatrix(GeneralMatrix &m, int i, int j, int nrows, int ncols);
  /* this = a*b */
  GeneralMatrix(const GeneralMatrix &a, const GeneralMatrix &b);
  /* this = a*b' */
  GeneralMatrix(const GeneralMatrix &a, const GeneralMatrix &b, const char *dum);
  /* this = a'*b */
  GeneralMatrix(const GeneralMatrix &a, const char *dum, const GeneralMatrix &b);
  /* this = a'*b */
  GeneralMatrix(const GeneralMatrix &a, const char *dum1,
                const GeneralMatrix &b, const char *dum2);

  virtual
  ~GeneralMatrix();
  const GeneralMatrix &
  operator=(const GeneralMatrix &m)
  {
    data = m.data; rows = m.rows; cols = m.cols; ld = m.ld; return *this;
  }

  const double &
  get(int i, int j) const
  {
    return data[j*ld+i];
  }
  double &
  get(int i, int j)
  {
    return data[j*ld+i];
  }
  int
  numRows() const
  {
    return rows;
  }
  int
  numCols() const
  {
    return cols;
  }
  int
  getLD() const
  {
    return ld;
  }
  double *
  base()
  {
    return data.base();
  }
  const double *
  base() const
  {
    return data.base();
  }
  Vector &
  getData()
  {
    return data;
  }
  const Vector &
  getData() const
  {
    return data;
  }

  double
  getNormInf() const
  {
    return ConstGeneralMatrix(*this).getNormInf();
  }
  double
  getNorm1() const
  {
    return ConstGeneralMatrix(*this).getNorm1();
  }

  /* place matrix m to the position (i,j) */
  void place(const ConstGeneralMatrix &m, int i, int j);
  void
  place(const GeneralMatrix &m, int i, int j)
  {
    place(ConstGeneralMatrix(m), i, j);
  }

  /* this = a*b */
  void mult(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b);
  void
  mult(const GeneralMatrix &a, const GeneralMatrix &b)
  {
    mult(ConstGeneralMatrix(a), ConstGeneralMatrix(b));
  }

  /* this = this + scalar*a*b */
  void multAndAdd(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b,
                  double mult = 1.0);
  void
  multAndAdd(const GeneralMatrix &a, const GeneralMatrix &b,
             double mult = 1.0)
  {
    multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), mult);
  }

  /* this = this + scalar*a*b' */
  void multAndAdd(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b,
                  const char *dum, double mult = 1.0);
  void
  multAndAdd(const GeneralMatrix &a, const GeneralMatrix &b,
             const char *dum, double mult = 1.0)
  {
    multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), dum, mult);
  }

  /* this = this + scalar*a'*b */
  void multAndAdd(const ConstGeneralMatrix &a, const char *dum, const ConstGeneralMatrix &b,
                  double mult = 1.0);
  void
  multAndAdd(const GeneralMatrix &a, const char *dum, const GeneralMatrix &b,
             double mult = 1.0)
  {
    multAndAdd(ConstGeneralMatrix(a), dum, ConstGeneralMatrix(b), mult);
  }

  /* this = this + scalar*a'*b' */
  void multAndAdd(const ConstGeneralMatrix &a, const char *dum1,
                  const ConstGeneralMatrix &b, const char *dum2, double mult = 1.0);
  void
  multAndAdd(const GeneralMatrix &a, const char *dum1,
             const GeneralMatrix &b, const char *dum2, double mult = 1.0)
  {
    multAndAdd(ConstGeneralMatrix(a), dum1, ConstGeneralMatrix(b), dum2, mult);
  }

  /* this = this + scalar*a*a' */
  void addOuter(const ConstVector &a, double mult = 1.0);
  void
  addOuter(const Vector &a, double mult = 1.0)
  {
    addOuter(ConstVector(a), mult);
  }

  /* this = this * m */
  void multRight(const ConstGeneralMatrix &m);
  void
  multRight(const GeneralMatrix &m)
  {
    multRight(ConstGeneralMatrix(m));
  }

  /* this = m * this */
  void multLeft(const ConstGeneralMatrix &m);
  void
  multLeft(const GeneralMatrix &m)
  {
    multLeft(ConstGeneralMatrix(m));
  }

  /* this = this * m' */
  void multRightTrans(const ConstGeneralMatrix &m);
  void
  multRightTrans(const GeneralMatrix &m)
  {
    multRightTrans(ConstGeneralMatrix(m));
  }

  /* this = m' * this */
  void multLeftTrans(const ConstGeneralMatrix &m);
  void
  multLeftTrans(const GeneralMatrix &m)
  {
    multLeftTrans(ConstGeneralMatrix(m));
  }

  /* x = scalar(a)*x + scalar(b)*this*d */
  void
  multVec(double a, Vector &x, double b, const ConstVector &d) const
  {
    ConstGeneralMatrix(*this).multVec(a, x, b, d);
  }

  /* x = scalar(a)*x + scalar(b)*this'*d */
  void
  multVecTrans(double a, Vector &x, double b, const ConstVector &d) const
  {
    ConstGeneralMatrix(*this).multVecTrans(a, x, b, d);
  }

  /* x = x + this*d */
  void
  multaVec(Vector &x, const ConstVector &d) const
  {
    ConstGeneralMatrix(*this).multaVec(x, d);
  }

  /* x = x + this'*d */
  void
  multaVecTrans(Vector &x, const ConstVector &d) const
  {
    ConstGeneralMatrix(*this).multaVecTrans(x, d);
  }

  /* x = x - this*d */
  void
  multsVec(Vector &x, const ConstVector &d) const
  {
    ConstGeneralMatrix(*this).multsVec(x, d);
  }

  /* x = x - this'*d */
  void
  multsVecTrans(Vector &x, const ConstVector &d) const
  {
    ConstGeneralMatrix(*this).multsVecTrans(x, d);
  }

  /* this = zero */
  void zeros();

  /** this = unit (on main diagonal) */
  void unit();

  /* this = NaN */
  void nans();

  /* this = Inf */
  void infs();

  /* this = scalar*this */
  void mult(double a);

  /* this = this + scalar*m */
  void add(double a, const ConstGeneralMatrix &m);
  void
  add(double a, const GeneralMatrix &m)
  {
    add(a, ConstGeneralMatrix(m));
  }

  /* this = this + scalar*m' */
  void add(double a, const ConstGeneralMatrix &m, const char *dum);
  void
  add(double a, const GeneralMatrix &m, const char *dum)
  {
    add(a, ConstGeneralMatrix(m), dum);
  }

  bool
  isFinite() const
  {
    return (ConstGeneralMatrix(*this)).isFinite();
  }

  bool
  isZero() const
  {
    return (ConstGeneralMatrix(*this)).isZero();
  }

  virtual void
  print() const
  {
    ConstGeneralMatrix(*this).print();
  }
private:
  void copy(const ConstGeneralMatrix &m, int ioff = 0, int joff = 0);
  void
  copy(const GeneralMatrix &m, int ioff = 0, int joff = 0)
  {
    copy(ConstGeneralMatrix(m), ioff, joff);
  }

  void gemm(const char *transa, const ConstGeneralMatrix &a,
            const char *transb, const ConstGeneralMatrix &b,
            double alpha, double beta);
  void
  gemm(const char *transa, const GeneralMatrix &a,
       const char *transb, const GeneralMatrix &b,
       double alpha, double beta)
  {
    gemm(transa, ConstGeneralMatrix(a), transb, ConstGeneralMatrix(b),
         alpha, beta);
  }

  /* this = this * op(m) (without whole copy of this) */
  void gemm_partial_right(const char *trans, const ConstGeneralMatrix &m,
                          double alpha, double beta);
  void
  gemm_partial_right(const char *trans, const GeneralMatrix &m,
                     double alpha, double beta)
  {
    gemm_partial_right(trans, ConstGeneralMatrix(m), alpha, beta);
  }

  /* this = op(m) *this (without whole copy of this) */
  void gemm_partial_left(const char *trans, const ConstGeneralMatrix &m,
                         double alpha, double beta);
  void
  gemm_partial_left(const char *trans, const GeneralMatrix &m,
                    double alpha, double beta)
  {
    gemm_partial_left(trans, ConstGeneralMatrix(m), alpha, beta);
  }

  /* number of rows/columns for copy used in gemm_partial_* */
  static int md_length;
};

class SVDDecomp
{
protected:
  /** Minimum of number of rows and columns of the decomposed
   * matrix. */
  const int minmn;
  /** Singular values. */
  Vector sigma;
  /** Orthogonal matrix U. */
  GeneralMatrix U;
  /** Orthogonal matrix V^T. */
  GeneralMatrix VT;
  /** Convered flag. */
  bool conv;
public:
  SVDDecomp(const GeneralMatrix &A)
    : minmn(std::min<int>(A.numRows(), A.numCols())),
      sigma(minmn),
      U(A.numRows(), A.numRows()),
      VT(A.numCols(), A.numCols()),
      conv(false)
  {
    construct(A);
  }
  const GeneralMatrix &
  getU() const
  {
    return U;
  }
  const GeneralMatrix &
  getVT() const
  {
    return VT;
  }
  void solve(const GeneralMatrix &B, GeneralMatrix &X) const;
  void
  solve(const Vector &b, Vector &x) const
  {
    GeneralMatrix xmat(x.base(), x.length(), 1);
    solve(GeneralMatrix(b.base(), b.length(), 1), xmat);
  }
private:
  void construct(const GeneralMatrix &A);
};

#endif /* GENERAL_MATRIX_H */

// Local Variables:
// mode:C++
// End: