Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!usc!snorkelwacker.mit.edu!bloom-picayune.mit.edu!athena.mit.edu!bjaspan
From: bjaspan@athena.mit.edu (Barr3y Jaspan)
Newsgroups: comp.lang.c++
Subject: Re: Matrix Class Libraries
Message-ID: <1991Mar18.025539.18998@athena.mit.edu>
Date: 18 Mar 91 02:55:39 GMT
References: <13505@darkstar.ucsc.edu>
Sender: news@athena.mit.edu (News system)
Organization: Massachusetts Institute of Technology
Lines: 446

I've had several requests for this, and it isn't very long, so I guess I
should just post it here.  I have found it very useful and reliable; your
results may vary.  Naturally, no docs are provided (UTFL). :-)  I don't claim
that is the most efficient or well designed or, in fact, anything.  

Please send any extensions/bugfixes back to me.

-- 
Barr3y Jaspan, bjaspan@mit.edu
Watchmaker Computing

---- snip Matrix.H snip ----

#ifndef __MATRIX_H__
#define __MATRIX_H__

/* $Id: Matrix.H,v 1.1 90/12/03 03:23:20 marc Exp Locker: marc $ */

#include <assert.h>  /* for abort() */
#include <stream.h>

class Matrix {
public:
     /* Constructors */
     Matrix(int, int);
     Matrix(int, int, int *);
     Matrix(int, int, double *);
     Matrix(int, int, int ...);
     Matrix(int, int, double ...);
     Matrix(const Matrix&);
     Matrix(const Matrix&,int,int,int,int);
     ~Matrix();

     /* Member functions */
     inline int Rows() const;
     inline int Cols() const;
     inline int SameSize(const Matrix&) const;
     inline double& operator()(int, int) const;
     inline double *Row(int) const;
     void Print() const;

     /* Overloaded operators */
     friend Matrix operator+(const Matrix&, const Matrix&);
     friend Matrix operator-(const Matrix&, const Matrix&);
     friend Matrix operator*(const Matrix&, const Matrix&);
     friend Matrix operator*(double, const Matrix&);
     Matrix& operator+=(const Matrix&);
     Matrix& operator-=(const Matrix&);
     inline friend Matrix& operator*=(Matrix&, const Matrix&);
     inline friend Matrix operator/(const Matrix&, double);
     friend int operator==(const Matrix&, const Matrix&);
     friend Matrix operator|(const Matrix&, const Matrix&);
     friend Matrix operator&(const Matrix&, const Matrix&);
     Matrix& operator=(const Matrix&);
     friend ostream& operator<<(ostream&, const Matrix&);

     /* Special matrix creators */
     static inline const Matrix& NAM();

     /* Matrix mutators */
     void Identify();

private:
     void Init(int, int);

     static Matrix *M_NAM;
     static Matrix *MakeNAM();
     static double RowDotColumn(const Matrix&, int i, const Matrix&, int j);
     
     int M_m, M_n, nam;
     double *d;
};

inline int Matrix::Rows() const
{
     return M_m;
}

inline int Matrix::Cols() const
{
     return M_n;
}

inline Matrix& operator*=(Matrix& m1, const Matrix& m)
{
   return(m1 = m1*m);
}

inline Matrix operator/(const Matrix& m, double k)
{
   return((1.0/k)*m);
}

inline double& Matrix::operator()(int m, int n) const
{
     /* My kingdom for exception handling ... */
     if (m > M_m || n > M_n || m < 1 || n < 1)
	  abort();
     
     return d[M_n*(m-1) + n - 1];
}

inline double *Matrix::Row(int r) const
{
   return(d+(r-1)*M_n);
}

inline const Matrix& Matrix::NAM()
{
     return *Matrix::M_NAM;
}

inline int Matrix::SameSize(const Matrix& m) const
{
     return (M_m == m.M_m && M_n == m.M_n);
}

/* DO NOT ADD ANYTHING AFTER THIS #endif */
#endif /* __MATRIX_H__ */

---- snip Matrix.C snip ----

/*
 * C++ class for MxN matrices.
 *
 * $Id: Matrix.C,v 1.2 90/12/04 01:38:15 bjaspan Exp Locker: marc $
 *
 */

/* An m x n matrix has m rows and n columns -- ie height x width */
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include "Matrix.H"

Matrix *Matrix::M_NAM = Matrix::MakeNAM();

/* Constructors */

void Matrix::Init(int m, int n)
{
     M_m = m;
     M_n = n;
     nam = 0;
     d = new double[m*n];
}

Matrix::Matrix(int m, int n)
{
     Init(m, n);
}

Matrix::Matrix(int m, int n, double first ...)
{
     int i;
     va_list ap;
     
     Init(m, n);

     d[0] = first;
     va_start(ap, first);
     for (i=1; i < m*n; i += 1)
	  d[i] = va_arg(ap, double);
}

Matrix::Matrix(int m, int n, int first ...)
{
     int i;
     va_list ap;
     
     Init(m, n);

     d[0] = double(first);
     va_start(ap, first);
     for (i=1; i < m*n; i += 1)
	  d[i] = double(va_arg(ap, int));
}

Matrix::Matrix(int m, int n, int *array)
{
     int i;

     Init(m, n);
     for (i=0; i < (m*n); i += 1)
	  d[i] = double(array[i]);
}

Matrix::Matrix(int m, int n, double *array)
{
     int i;

     Init(m, n);
     for (i=0; i < (m*n); i += 1)
	  d[i] = array[i];
}

Matrix::Matrix(const Matrix& m1)
{
     int i;

     Init(m1.M_m, m1.M_n);
     for (i=0; i < (M_m*M_n); i += 1)
	  d[i] = m1.d[i];
}

Matrix::Matrix(const Matrix& m1, int r, int c, int m, int n)
{
     int i,j;

     if (r<1 || c<1 || m1.Rows() < r-1+m || m1.Cols() < c-1+n)
	  nam = 1;
     else {
	  Init(m,n);

	  for (i=1;i<=m;i++)
	       for (j=1;j<=n;j++)
		    (*this)(i,j) = m1(r+i-1,c+j-1);
     }
}

Matrix::~Matrix()
{
     delete [M_m*M_n] d;
}

/* Overloaded operators */

int operator==(const Matrix& m1, const Matrix& m2)
{
     int i;

     if (m1.nam != m2.nam)
	  return 0;
     if (m1.nam == 1)
	  return 1;
     if (! m1.SameSize(m2))
	  return 0;

     for (i = 0; i < m1.M_m*m1.M_n; i++)
	  if (m1.d[i] != m2.d[i])
	       return 0;

     return 1;
}

Matrix operator+(const Matrix& m1, const Matrix& m2)
{
     if (m1.M_m != m2.M_m || m1.M_n != m2.M_n ||
	 m1 == Matrix::NAM() || m2 == Matrix::NAM())
	  return Matrix::NAM();

     Matrix m(m1);

     for (int i=0; i < (m1.M_m*m1.M_n); i++)
	  m.d[i] += m2.d[i];

     return m;
}

Matrix operator-(const Matrix& m1, const Matrix& m2)
{
     if (m1.M_m != m2.M_m || m1.M_n != m2.M_n ||
	 m1 == Matrix::NAM() || m2 == Matrix::NAM())
	  return Matrix::NAM();

     Matrix m(m1);

     for (int i=0; i < (m1.M_m*m1.M_n); i++)
	  m.d[i] -= m2.d[i];

     return m;
}

Matrix operator*(const Matrix& m1, const Matrix& m2)
{
     if (m1.M_n != m2.M_m || m1 == Matrix::NAM() || m2 == Matrix::NAM())
	  return Matrix::NAM();

     Matrix m(m1.M_m, m2.M_n);

     for (int i=1; i <= m1.M_m; i += 1)
	  for (int j=1; j <= m2.M_n; j += 1)
	       for (int k=1; k <= m1.M_n; k += 1)
		    m(i, j) += m1(i, k) * m2(k, j);
     return m;
}

Matrix& Matrix::operator+=(const Matrix& m)
{
     if (M_m != m.M_m || M_n != m.M_n || *this == Matrix::NAM() ||
	 m == Matrix::NAM())
	  nam = 1;
     else
	  for (int i=0; i < (M_m*M_n); i++)
	       d[i] += m.d[i];

     return *this;
}

Matrix& Matrix::operator-=(const Matrix& m)
{
     if (M_m != m.M_m || M_n != m.M_n || *this == Matrix::NAM() ||
	 m == Matrix::NAM())
	  nam = 1;
     else
	  for (int i=0; i < (M_m*M_n); i++)
	       d[i] -= m.d[i];

     return *this;
}

Matrix operator*(double k, const Matrix& m)
{
   Matrix n(m);
   int i;

   for (i=0;i<n.M_m*n.M_n;i++) n.d[i]*=k;
   return(n);
}

Matrix& Matrix::operator=(const Matrix& m)
{
     if (M_m == m.M_m && M_n == m.M_n)
	  bcopy(m.d, d, M_m*M_n*sizeof(double));
     else if (M_m*M_n == m.M_m*m.M_n) {
	  bcopy(m.d, d, M_m*M_n*sizeof(double));
	  M_m = m.M_m;
	  M_n = m.M_n;
     } else {
	  delete [M_m*M_n] d;
	  
	  M_m = m.M_m;
	  M_n = m.M_n;
	  
	  d = new double[M_m*M_n];
	  bcopy(m.d, d, M_m*M_n*sizeof(double));
     }

     nam = m.nam;

     return(*this);
}

Matrix operator|(const Matrix& m1, const Matrix& m2)
   // Concatenates two matrices horizontally.  Requires that m1 and m2
   // have the same number of rows.
   // [ A ] | [ B ] --> [ A | B ]
{
     if (m1.Rows() != m2.Rows() || m1 == Matrix::NAM() || m2 == Matrix::NAM())
	  return Matrix::NAM();

     Matrix m(m1.Rows(), m1.Cols() + m2.Cols());
     int r, c;

     for (r = 1; r <= m1.Rows(); r += 1)
	  for (c = 1; c <= m1.Cols(); c += 1)
	       m(r, c) = m1(r, c);

     for (r = 1; r <= m1.Rows(); r += 1)
	  for (c = 1; c <= m2.Cols(); c += 1)
	       m(r, c + m1.Cols()) = m2(r, c);

     return m;
}

Matrix operator&(const Matrix& m1, const Matrix& m2)
   // Concatenates two matrices vertically.  Requires that m1 and m2
   // have the same number of columns.
   //                   [ A ]
   // [ A ] & [ B ] --> [ - ]
   //                   [ B ]
{
     if (m1.Cols() != m2.Cols() || m1 == Matrix::NAM() || m2 == Matrix::NAM())
	  return Matrix::NAM();

     Matrix m(m1.Rows()+m2.Rows(), m1.Cols());
     int r, c;

     for (r = 1; r <= m1.Rows(); r += 1)
	  for (c = 1; c <= m1.Cols(); c += 1)
	       m(r, c) = m1(r, c);
     for (r = 1; r <= m2.Rows(); r += 1)
	  for (c = 1; c <= m2.Cols(); c += 1)
	       m(r + m1.Rows(), c) = m2(r, c);
     
#if 0
     // This doesn't work, and I don't know why.
     bcopy(m1.d, m.d, m1.Rows()*m1.Cols()*sizeof(double));
     bcopy(m2.d, m.d + m1.Rows()*m1.Cols()*sizeof(double),
	   m2.Rows()*m2.Cols()*sizeof(double));
#endif

     return m;
}

ostream& operator<<(ostream& s, const Matrix& m)
{
     int i,j;
     
     if (m.nam) {
	  s << "NAM\n";
	  return(s);
     }
     
     for (i = 1; i <= m.M_m; i++) {
	  for (j = 1; j <= m.M_n; j++)
#ifdef __GNUG__
	       s.form("%6.2f ", m(i, j));
#else
	       s << form("%6.2f ", m(i, j));
#endif
	  s << "\n";
     }
     
     return(s);
}

/* Member functions */

void Matrix::Print() const
{
   cout << *this;
}

/* Special matrix constructors */
Matrix *Matrix::MakeNAM()
{
     Matrix *m = new Matrix(1,1);
     m->nam = 1;
     return m;
}

/* Matrix mutators */
void Matrix::Identify()	
{
     int i;
     
     if (M_m != M_n) {
	  nam = 1;
	  return;
     }

     bzero((char *) d, sizeof(double)*M_m*M_n);
     for (i = 1; i <= M_m; i++)
	  (*this)(i, i) = 1;
}