A C++ Matrix Template Class

Listing 1: Matrix template class

#ifndef matrix_h
#define matrix_h

#include <stdexcept>
#include <vector>
#include <function>

// undefine to disable range checking
#define RANGE_CHECK

class overdetermined : public std::domain_error
{
   public:
     overdetermined() 
        : std::domain_error("solution is over-determined")
     {}
};

class underdetermined : public std::domain_error
{
   public:
      underdetermined() 
         : std::domain_error("solution is under-determined")
      {}
};

template<class T>
class kahn_sum {     // implements Kahn Summation method
   public:
      kahn_sum() : sum(0.0), cor(0.0) {}
      kahn_sum<T>& operator+=( const T& val ) {
         T old_sum = sum;
         T next = val-cor;
         cor = ( (sum += next) - old_sum ) - next;
         return *this;
      }
      kahn_sum<T>& operator-=( const T& val ) {
         T old_sum = sum;
         T next = val+cor;
         cor = ( (sum -= val) - old_sum ) + next;
         return *this;
      }
      operator T&() { return sum; }
   private:
      T sum;  // running sum
      T cor;  // correction term
};

template<class T>
class matrix {
   private:
      std::vector<T> elements; // array of elements
   public:
      const unsigned rows;  // number of rows
      const unsigned cols;  // number of columns
      
   protected:
      // range check function for matrix access
      void range_check( unsigned i, unsigned j ) const;
   public:
      T& operator()( unsigned i, unsigned j ) {
         #ifdef RANGE_CHECK
         range_check(i,j);
         #endif
         return elements[i*cols+j];
      }
      const T& operator()( unsigned i, unsigned j ) const {
         #ifdef RANGE_CHECK
         range_check(i,j);
         #endif
         return elements[i*cols+j];
      }
      const T& element(unsigned i, unsigned j) const {
         #ifdef RANGE_CHECK
         range_check(i,j);
         #endif
         return elements[i*cols+j];
      }
      T& element(unsigned i, unsigned j) {
         #ifdef RANGE_CHECK
         range_check(i,j);
         #endif
         return elements[i*cols+j];
      }
   public:
      // constructors
      matrix( unsigned rows, unsigned cols, const T* elements = 0 );
      matrix( const matrix<T>& );
      // destructor
      ~matrix();

      // assignment
      matrix<T>& operator=( const matrix<T>& );

      // comparison
      bool operator==( const matrix<T>& ) const;
      bool iszero() const;
      bool operator!() const {
        return iszero();
      }

      // scalar multiplication/division
      matrix<T>& operator*=( const T& a );
      matrix<T> operator*( const T& a ) const {
         return matrix<T>(*this).operator*=(a);
      }
      matrix<T>& operator/=( const T& a );
      matrix<T> operator/( const T& a ) {
        return matrix<T>(*this).operator/=(a);
      }
      matrix<T> operator-() const;
      matrix<T> operator+() const;

      // addition/subtraction
      matrix<T>& operator+=( const matrix<T>& );
      matrix<T>& operator-=( const matrix<T>& );
      matrix<T> operator+( const matrix<T>& M ) const {
         return matrix<T>(*this).operator+=(M);
      }
      matrix<T> operator-( const matrix<T>& M ) const {
         return matrix<T>(*this).operator-=(M);
      }

      // matrix multiplication
      matrix<T> operator*( const matrix<T>& ) const;
      matrix<T>& operator*=( const matrix<T>& M ) {
         return *this = *this * M;
      }
      // matrix division
      matrix<T> leftdiv( const matrix<T>& ) const;
      matrix<T> rightdiv( const matrix<T>& D ) const {
         return transpose().leftdiv(D.transpose()).transpose();
      }
      matrix<T> operator/( const matrix<T>& D ) const {
         return rightdiv(D);
      }
      matrix<T>& operator/=( const matrix<T>& M ) {
         return *this = *this/M;
      }

      // determinants
      matrix<T> minor( unsigned i, unsigned j ) const;
      T det() const;
      T minor_det( unsigned i, unsigned j ) const;

      // these member functions are only valid for squares
      matrix<T> inverse() const;
      matrix<T> pow( int exp ) const;
      matrix<T> identity() const;
      bool isidentity() const;

      // vector operations
      matrix<T> getrow( unsigned j ) const;
      matrix<T> getcol( unsigned i ) const;
      matrix<T>& setcol( unsigned j, const matrix<T>& C );
      matrix<T>& setrow( unsigned i, const matrix<T>& R );
      matrix<T> delrow( unsigned i ) const;
      matrix<T> delcol( unsigned j ) const;

      matrix<T> transpose() const;
      matrix<T> operator~() const {
         return transpose();
      }
};

template<class T>
matrix<T>::matrix( unsigned rows, unsigned cols, const T* elements = 0 )
  : rows(rows), cols(cols), elements(rows*cols,T(0.0))
{
   if( rows == 0 || cols == 0 )
      throw std::range_error("attempt to create a degenerate matrix");
   // initialze from array
   if(elements)
      for(unsigned i=0;i<rows*cols;i++)
         this->elements[i] = elements[i];
};

template<class T>
matrix<T>::matrix( const matrix<T>& cp )
  : rows(cp.rows), cols(cp.cols), elements(cp.elements)
{
}

template<class T>
matrix<T>::~matrix()
{
}

template<class T>
matrix<T>& matrix<T>::operator=( const matrix<T>& cp )
{
   if(cp.rows != rows && cp.cols != cols )
      throw std::domain_error("matrix op= not of same order");
   for(unsigned i=0;i<rows*cols;i++)
      elements[i] = cp.elements[i];
   return *this;
}
template<class T>
void matrix<T>::range_check( unsigned i, unsigned j ) const
{
   if( rows <= i )
      throw std::range_error("matrix access row out of range");
   if( cols <= j )
      throw std::range_error("matrix access col out of range");
}

//not shown: operator==(const matrix<T>& A) const, iszero() const,
//operator*=(const T& a), operator/=(const T& a), operator-() const,
//operator+() const, operator*(const T& a, const matrix<T>& M),
//operator+=(const matrix<T>& M), operator-=(const matrix<T>& M),
//minor(unsigned i, unsigned j) const,
//minor_det(unsigned i, unsigned j) const, det() const -- mb
//operator*( const matrix<T>& B) const ... mb

template<class T>
matrix<T> matrix<T>::leftdiv( const matrix<T>& D ) const
{
   const matrix<T>& N = *this;

   if( N.rows != D.rows )
      throw std::domain_error("matrix divide: incompatible orders");

   matrix<T> Q(D.cols,N.cols);  // quotient matrix

   if( N.cols > 1 ) {
       // break this down for each column of the numerator
       for(unsigned j=0;j<Q.cols;j++)
          Q.setcol( j, N.getcol(j).leftdiv(D) );  // note: recursive
       return Q;
   }

   // from here on, N.col == 1

   if( D.rows < D.cols )
      throw underdetermined();

   if( D.rows > D.cols ) {
      bool solution = false;
      for(unsigned i=0;i<D.rows;i++) {
         matrix<T> D2 = D.delrow(i);  // delete a row from the matrix
         matrix<T> N2 = N.delrow(i);
         matrix<T> Q2(Q);  
         try {
            Q2 = N2.leftdiv(D2);
         } 
         catch( underdetermined x ) {
            continue;  // try again with next row
         }
         if( !solution ) {
            // this is our possible solution
            solution = true;
            Q = Q2;
         } else {
            // do the solutions agree?
            if( Q != Q2 )
               throw overdetermined();
         }
      }
      if( !solution )
         throw underdetermined();
      return Q;
   }

   // D.rows == D.cols && N.cols == 1
   // use Kramer's Rule
   //
   // D is a square matrix of order N x N
   // N is a matrix of order N x 1

   const T T0(0.0); // additive identity

   if( D.cols<=3 ) {
      T ddet = D.det();
      if( ddet == T0 )
         throw underdetermined();
      for(unsigned j=0;j<D.cols;j++) {
         matrix<T> A(D); // make a copy of the D matrix
         // replace column with numerator vector
         A.setcol(j,N);
         Q(j,0) = A.det()/ddet;
      }
   } else {
      // this method optimizes the determinant calculations
      // by saving a cofactors used in calculating the
      // denominator determinant.
      kahn_sum<T> sum;
      vector<T> cofactors(D.cols);  // save cofactors
      for(unsigned j=0;j<D.cols;j++) {
         T c = D.minor_det(0,j);
         cofactors[j] = c;
         T a = D(0,j);
         if( a != T0 ) {
            a *= c;
            if(j%2)
               sum -= a;
            else
               sum += a;
         }
      }
      T ddet = sum;
      if( ddet == T0 )
         throw underdetermined();
      for(unsigned j = 0;j<D.cols;j++) {
         matrix<T> A(D);
         A.setcol(j,N);
         kahn_sum<T> ndet;
         for(unsigned k=0;k<D.cols;k++) {
            T a = A(0,k);
            if( a != T0 ) {
               if(k==j)
                  a *= cofactors[k];  // use previously calculated
                                      // cofactor
               else
                  a *= A.minor_det(0,k); // calculate minor's determinant
               if(k%2)
                  ndet -= a;
               else
                  ndet += a;
            }
         }
         Q(j,0) = T(ndet)/ddet;
      }
   }
   return Q;
}

//not shown: inverse() const, getrow(unsigned i) const,
//getcol(unsigned j) const, delcol(unsigned j) const,
//delrow(unsigned i) const, setcol(unsigned j, const matrix<T>& C)
//setrow(unsigned i, const matrix<T>& R), identity() const,
//isidentity() const, pow(int exp) const,
//pow(const matrix<T>& M, int exp) --mb
// ...
#endif
/* End of File */