Tools/Math/LA_Matrix.cpp

Go to the documentation of this file.

//------------------------------------------------------------------------------
#include <vector>
#include <iomanip>
#include <sstream>
#include <assert.h>
#include <float.h>
#include <math.h>
#include "LA_Matrix.h"
#include "LA_Lu_decomp.h"
//------------------------------------------------------------------------------
namespace LA
{
//------------------------------------------------------------------------------
Matrix::Matrix() : r(0), c(0), content(NULL)
{
}
//------------------------------------------------------------------------------
Matrix::Matrix(int rows, int columns) : r(0), c(0), content(NULL)
{
  create(rows, columns);
}
//------------------------------------------------------------------------------
Matrix::Matrix(int rows, int columns, const double* M)
: r(0), c(0), content(NULL)
{
  create(rows, columns, M);
}
//------------------------------------------------------------------------------
Matrix::Matrix(const Matrix& M)
: r(0), c(0), content(NULL)
{
  *this = M;
}
//------------------------------------------------------------------------------
Matrix::Matrix(const Matrix3x3<double>& M)
: r(0), c(0), content(NULL)
{
  *this = M;
}
//------------------------------------------------------------------------------
Matrix Matrix::identity(int dim)
{
  Matrix M(dim, dim);
  for (int i = 1; i <= dim; ++i)
    M(i, i) = 1.0;
  return M;
}
//------------------------------------------------------------------------------
Matrix::~Matrix()
{
  destroy();
}
//------------------------------------------------------------------------------
Matrix& Matrix::operator=(const Matrix& M)
{
  create(M.r, M.c);
  for (int i = 0; i < M.r; ++i)
    content[i] = M.content[i];
  return *this;
}
//------------------------------------------------------------------------------
Matrix& Matrix::operator=(const Matrix3x3<double>& M)
{
  create(3, 3);
  for (int i = 0; i < 3; ++i)
  {
    content[0][i] = M.c[i].x;
    content[1][i] = M.c[i].y;
    content[2][i] = M.c[i].z;
  }
  return *this;
}
//------------------------------------------------------------------------------
void Matrix::create(int rows, int columns)
{
  destroy();
  r = rows;
  c = columns;
  if (r != 0)
  {
    content = new Vector[r];
    for (int i = 0; i < r; ++i)
      content[i] = Vector(c);
  }
}
//------------------------------------------------------------------------------
void Matrix::create(int rows, int columns, const double* M)
{
  destroy();
  r = rows;
  c = columns;
  if (r != 0)
  {
    content = new Vector[r];
    const double* pos = M;
    for (int i = 0; i < r; ++i)
    {
      content[i] = Vector(c, pos);
      pos += c;
    }
  }
}
//------------------------------------------------------------------------------
const Vector& Matrix::operator[](int row) const
{
  assert((row >= 0) && (row < r));
  return content[row];
}
//------------------------------------------------------------------------------
Vector& Matrix::operator[](int row)
{
  assert((row >= 0) && (row < r));
  return content[row];
}
//------------------------------------------------------------------------------
const Vector& Matrix::operator()(int row) const
{
  assert((row >= 1) && (row <= r));
  return content[row-1];  
}
//------------------------------------------------------------------------------
Vector& Matrix::operator()(int row)
{
  assert((row >= 1) && (row <= r));
  return content[row-1];  
}
//------------------------------------------------------------------------------
const double& Matrix::operator()(int row, int column) const
{
  assert((row >= 1) && (row <= r));
  assert((column >= 1) && (column <= c));
  return content[row-1].content[column-1];
}
//------------------------------------------------------------------------------
double& Matrix::operator()(int row, int column)
{
  assert((row >= 1) && (row <= r));
  assert((column >= 1) && (column <= c));
  return content[row-1].content[column-1];
}
//------------------------------------------------------------------------------
int Matrix::rows() const
{
  return r;
}
//------------------------------------------------------------------------------
int Matrix::columns() const
{
  return c;
}
//------------------------------------------------------------------------------
Vector Matrix::get_row(int index) const
{
  assert((index >= 1) && (index <= r));
  return content[index-1];
}
//------------------------------------------------------------------------------
void Matrix::set_row(int index, const Vector& v)
{
  assert((index >= 1) && (index <= r));
  assert(v.n == c);
  content[index-1] = v;
}
//------------------------------------------------------------------------------
Vector Matrix::get_column(int index) const
{
  assert((index >= 1) && (index <= c));
  Vector v(c);
  index--;
  for (int i = 0; i < c; ++i)
    v[i] = content[i][index];    
  return v;
}
//------------------------------------------------------------------------------
void Matrix::set_column(int index, const Vector& v)
{
  assert((index >= 1) && (index <= c));
  assert(v.n == r);
  index--;
  for (int i = 0; i < r; ++i)
    content[i].content[index] = v.content[i];
}
//------------------------------------------------------------------------------
Matrix Matrix::get_sub(int rowstart, int rowlength,
                       int columnstart, int columnlength) const
{
  assert((rowstart >= 1) && ((rowstart+rowlength-1) <= r));
  assert((columnstart >= 1) && ((columnstart+columnlength-1) <= c));

  Matrix M(rowlength, columnlength);
  for (int i = 0; i < rowlength; ++i)
    M.content[i] = content[rowstart + i - 1].sub(columnstart, columnlength);

  return M;
}
//------------------------------------------------------------------------------
void Matrix::set_sub(int rowstart, int columnstart, const Matrix& M)
{
  assert((rowstart >= 1) && ((rowstart+M.r-1) <= r));
  assert((columnstart >= 1) && ((columnstart+M.c-1) <= c));

  for (int i = 0; i < M.r; ++i)
    content[rowstart + i - 1].set_sub(columnstart, M.content[i]);
}
//------------------------------------------------------------------------------
Matrix& Matrix::operator+=(const Matrix& M)
{
  assert((r == M.r) && (c == M.c));
  for (int i = 0; i < r; ++i)
    content[i] += M.content[i];
  return *this;
}
//------------------------------------------------------------------------------
Matrix& Matrix::operator-=(const Matrix& M)
{
  assert((r == M.r) && (c == M.c));
  for (int i = 0; i < r; ++i)
    content[i] -= M.content[i];
  return *this;
}
//------------------------------------------------------------------------------
Matrix& Matrix::operator*=(const Matrix& M)
{
  *this = *this * M;
  return *this;
}
//------------------------------------------------------------------------------
Matrix& Matrix::operator*=(const double& s)
{
  for (int i = 0; i < r; ++i)
    content[i] *= s;
  return *this;
}
//------------------------------------------------------------------------------
Matrix& Matrix::operator/=(const Matrix& M)
{
  *this *= LA::invert(M);
  return *this;
}
//------------------------------------------------------------------------------
Matrix& Matrix::operator/=(const double& s)
{
  for (int i = 0; i < r; ++i)
    content[i] /= s;
  return *this;
}
//------------------------------------------------------------------------------
void Matrix::transpose()
{
  Matrix MT(c, r);
  for (int i = 0; i < r; ++i)
  {
    for (int j = 0; j < c; ++j)
    {
      MT.content[j].content[i] = content[i].content[j];
    }
  }
  *this = MT;
}
//------------------------------------------------------------------------------
Matrix Matrix::transposition() const
{
  Matrix M(*this);
  M.transpose();
  return M;
}
//------------------------------------------------------------------------------
Matrix Matrix::T() const
{
  return transposition();
}
//------------------------------------------------------------------------------
void Matrix::invert()
{
  assert(r == c);

  LU lu(*this);
  *this = lu.inverse();
}
//------------------------------------------------------------------------------
Matrix Matrix::inverse() const
{
  return LA::invert(*this);
}
//------------------------------------------------------------------------------
Vector Matrix::solve(const Vector& b)
{
  assert(r == c);
  assert(b.n == r);

  LU lu(*this);
  return lu.solve(b);
}
//------------------------------------------------------------------------------
double Matrix::det() const
{
  assert(r == c);

  LU lu(*this);
  return lu.det();
}
//------------------------------------------------------------------------------
std::string Matrix::asString() const
{
  std::ostringstream os;
  os << *this;
  return os.str();
}
//------------------------------------------------------------------------------
void Matrix::destroy()
{
  if (content != NULL)
  {
    delete [] content;
    content = NULL;
  }
}
//------------------------------------------------------------------------------
Matrix operator+(const Matrix& M1, const Matrix& M2)
{
  Matrix M(M1);
  M += M2;
  return M;
}
//------------------------------------------------------------------------------
Matrix operator-(const Matrix& M1, const Matrix& M2)
{
  Matrix M(M1);
  M -= M2;
  return M;
}
//------------------------------------------------------------------------------
Matrix operator*(const Matrix& M1, const Matrix& M2)
{
  assert(M1.c == M2.r);
  Matrix M(M1.r, M2.c);

  for (int i = 0; i < M1.r; ++i)
  {
    for (int j = 0; j < M2.c; ++j)
    {
      for (int k = 0; k < M1.c; ++k)
      {
        M.content[i].content[j] +=
          M1.content[i].content[k]*M2.content[k].content[j];
      }
    }
  }

  return M;
}
//------------------------------------------------------------------------------
Vector operator*(const Matrix& M, const Vector& v)
{
  assert(M.c == v.n);
  Vector res(M.r);

  for (int i = 0; i < M.r; ++i)
  {
    for (int j = 0; j < M.c; ++j)
    {
      res.content[i] += M.content[i].content[j] * v.content[j];
    }
  }

  return res;
}
//------------------------------------------------------------------------------
Matrix operator*(const Matrix& M, const double& s)
{
  Matrix Mres(M);
  Mres *= s;
  return Mres;
}
//------------------------------------------------------------------------------
Matrix operator*(const double& s, const Matrix& M)
{
  Matrix Mres(M);
  Mres *= s;
  return Mres;
}
//------------------------------------------------------------------------------
Matrix operator/(const Matrix& M, const double& s)
{
  Matrix Mres(M);
  Mres /= s;
  return Mres;
}
//------------------------------------------------------------------------------
Matrix operator/(const Matrix& M1, const Matrix& M2)
{
  return M1 * invert(M2);
}
//------------------------------------------------------------------------------
Matrix transpose(const Matrix& M)
{
  Matrix MT(M);
  MT.transpose();
  return MT;
}
//------------------------------------------------------------------------------
Matrix invert(const Matrix& M)
{
  Matrix Minv(M);
  Minv.invert();
  return Minv;
}
//------------------------------------------------------------------------------
Matrix inv(const Matrix& M)
{
  Matrix Minv(M);
  Minv.invert();
  return Minv;
}
//------------------------------------------------------------------------------
double det(const Matrix& M)
{
  return M.det();
}
//------------------------------------------------------------------------------
std::ostream& operator<<(std::ostream& out, const Matrix& M)
{
  for (int i = 1; i <= M.rows(); ++i)
  {
    out << "| ";
    for (int j = 1; j <= M.columns(); ++j)
    {
      out << std::setw(10) << std::setfill(' ') << M(i, j) << ' ';      
    }
    out << "|" << std::endl;
  }
  return out;  
}
//------------------------------------------------------------------------------
Out& operator<<(Out& stream, const Matrix& M)
{
  stream << M.r;
  stream << M.c;
  for (int i = 0; i < M.r; ++i)
    for (int j = 0; j < M.c; ++j)
      stream << M.content[i].content[j];
  return stream;
}
//------------------------------------------------------------------------------
In& operator>>(In& stream, Matrix& M)
{
  int r, c;
  stream >> r;
  stream >> c;
  M.create(r, c);
  for (int i = 0; i < r; ++i)
    for (int j = 0; j < c; ++j)
      stream >> M.content[i].content[j];
  return stream;
}
//------------------------------------------------------------------------------
}
//------------------------------------------------------------------------------

---