matrix.h

用C++写的矩阵和矢量运算库

// -*-c++-*-
#ifndef IG_MATRIX_H
#define IG_MATRIX_H

// $Id: matrix.h,v 1.37 2003/01/13 18:08:48 hkuiper Exp $

// CwMtx matrix and vector math library
// Copyright (C) 1999-2001  Harry Kuiper
// Copyright (C) 2000  Will DeVore(template conversion)

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA

#include <iostream>
#include <iomanip>

// CWMatrix class

// This library was designed to mirror as closely as possible the
// notation used in mathematical writing. A matrix is indexed:
// matMatrixName[row][col].

// CAUTION!!!

// This matrix library was implemented with emphasis on
// speed. Consequently no attempts were made to trap and report
// errors. It is left entirely to the user to write code that does not
// cause errors within the matrix routines.

namespace CwMtx
{
  using std::ostream;

  // matrix status values
  enum { N_NOTALLOCATED, N_ALLOCATED, N_MAPPED };

  // Classes that create unity and zero "objects".  The ones directly
  // below work only for template arguments that are basic (numerical)
  // types that can be initialised by a literal 0 or 1.
  template <class T> class CWTUnity
  {
  public:
    operator T() { return 1; }
  };

  template <class T> class CWTZero
  {
  public:
    operator T() { return 0; }
  };

  // This template defaults to double. Most of the time this template
  // will be working with math functions that only work with
  // doubles. For example, the transcendental function sin(x) takes
  // and returns a double which would force the compiler to convert
  // back and forth from some other data type to a double.

  // prefix mat
  template < class T = double >
  class CWTMatrix
  {
  public:
    typedef T element;

    // creates a matrix, does NOT allocate rows and columns
    CWTMatrix();
    // creates a matrix, allocates rows and columns
    CWTMatrix(unsigned, unsigned);
    CWTMatrix(const CWTMatrix &);
    // sub-matrix mapped into another
    CWTMatrix(const CWTMatrix &, unsigned, unsigned, unsigned, unsigned);
	
    // removes matrix elements from free store
    ~CWTMatrix() { Deallocate(); };

    // allocates rows and colums
    void Dimension(unsigned, unsigned);
    // maps matrix into another
    void MapInto(const CWTMatrix&, unsigned, unsigned, unsigned, unsigned);
    // reverses the effect of Dimension() and MapInto()
    void Deallocate();

    int GetStatus() const { return m_nMatStatus; };
    unsigned GetRows() const { return m_crow; };
    unsigned GetCols() const { return m_ccol; };

    // basic matrix operations

    // returns a row of modifyable elements
    T* operator [](unsigned irow) { return m_rgrow[irow]; };
    // returns a row of non-modifyable elements
    const T* operator [](unsigned irow) const { return m_rgrow[irow]; };

    CWTMatrix operator +(const CWTMatrix &) const;
    CWTMatrix operator -(const CWTMatrix &) const;
    CWTMatrix operator -() const;
    CWTMatrix operator *(const T &) const;
    CWTMatrix operator *(const CWTMatrix &) const;

    // Interesting note here. Because we are defining the "/" operator
    // we can't expect to use it inside the definition unless we
    // safely restrict it to operating on constants. Without the
    // parens the operator does the "* 1" first then does the
    // "/value" which then leads to calling the "/" etc...  With the
    // parens, the intended scalar operation, "1/value", occurs
    // first then the Matrix/Scalar operations occur.  If the parens
    // are missing in the operator below, an ifinite loop
    // occurs. -----------------------------------------------V----------V
    CWTMatrix operator /(const T &value) const
    {
      return (*this)*static_cast<const T &>(CWTUnity<T>()/value);
    }

    // not inherited
    CWTMatrix & operator =(const CWTMatrix &);
    CWTMatrix & operator +=(const CWTMatrix &);
    CWTMatrix & operator -=(const CWTMatrix &);
    CWTMatrix & operator *=(const T &);
    CWTMatrix & operator /=(const T &value)
    {
      return (*this) *= static_cast<const T &>(CWTUnity<T>()/value);
    }

    int operator ==(const CWTMatrix &) const;
    int operator !=(const CWTMatrix &mat) const { return !( (*this) == mat ); }

    // stores CWMatrix + CWMatrix in this
    void StoreSum(const CWTMatrix &, const CWTMatrix &);
    // stores CWMatrix*CWMatrix in this
    void StoreProduct(const CWTMatrix &, const CWTMatrix &);
    // stores transpose of CWMatrix in this
    void StoreTranspose(const CWTMatrix &);
    // stores CWMatrix at indicated position in this
    void StoreAtPosition(unsigned, unsigned, const CWTMatrix &);
    // fills the whole array with a value.
    void Fill(const T &);

    void InterchangeRows(unsigned, unsigned);
    void AddRowToRow(unsigned, unsigned, const T & = CWTUnity<T>());
    void MultiplyRow(unsigned, const T &);

  private:
    // initializes data members
    void Initialize();

    // we keep the data structures used for CWMatrix implementation
    // private

    // row count
    unsigned m_crow;
    // column count
    unsigned m_ccol;
    // an array of rows (stored on free store)
    T **m_rgrow;
    // matrix status
    int m_nMatStatus;
  };

  // Templates to create self-dimensioning CWTMatrix classes - or one
  // of its derived classes - using the syntax of a default
  // constructor.  This facility is required for using matrices as
  // elements of matrices since these will always be created by a call
  // to the default constructor.

  template <class T, unsigned crow, unsigned ccol>
  class CWTMat: public T
  {
  public:
    CWTMat(): T(crow, ccol) {}

    T & operator =(const T &mtx) { return T::operator=(mtx); }
  };

  // NOTE: There exists no unity matrix for a non-square matrix!

  // Zero matrix.  NOTE: A zero matrix can only be constructed for a
  // matrix of known dimensions.  Hence the use of CWTMat<T,n,m>.
  template <class T, unsigned crow, unsigned ccol>
  class CWTZero< CWTMat<CWTMatrix<T>, crow, ccol> >:
    public CWTMat<CWTMatrix<T>, crow, ccol>
  {
  public:
    CWTZero() { Fill(CWTZero<T>()); }
  };

  //
  // Constructors
  //

  template < class T >
  inline CWTMatrix<T>::CWTMatrix()
  {
    Initialize();
  };

  // creates a matrix, allocates rows and columns
  template < class T >
  inline CWTMatrix<T>::CWTMatrix(unsigned crow, unsigned ccol)
  {
    Initialize();
    Dimension(crow, ccol);
  }

  template < class T >
  inline CWTMatrix<T>::CWTMatrix(const CWTMatrix<T> &mat)
  {
    Initialize();

    if (mat.m_nMatStatus == N_NOTALLOCATED)
      {
	// input matrix not allocated, so there's nothing to copy
	return;
      }
    else
      {
	// copy contents of input matrix
	(*this) = mat;
      }
  }

  // Mapped matrix constructor
  template < class T >
  inline CWTMatrix<T>::CWTMatrix(const CWTMatrix<T> &mat,
				 unsigned irowStart,
				 unsigned icolStart,
				 unsigned irowEnd,
				 unsigned icolEnd)
  {
    Initialize();
    MapInto(mat, irowStart, icolStart, irowEnd, icolEnd);
  }

  //
  // Private Methods
  //

  // initial values for matrix attributes
  template < class T >
  inline void CWTMatrix<T>::Initialize()
  {
    m_crow       = 0;
    m_ccol       = 0;
    m_rgrow      = NULL;
    m_nMatStatus = N_NOTALLOCATED;
  }

  //
  // User Methods
  //

  template < class T >
  void CWTMatrix<T>::Dimension(unsigned crowInit, unsigned ccolInit)
  {
    if (m_nMatStatus != N_NOTALLOCATED)
      {
	Deallocate();
      }

    m_crow = crowInit;
    m_ccol = ccolInit;

#ifdef CC_CWTMTX_ASSUME_BASIC_TYPES
    // NOTE: CWTMatrix normally uses the standard C++ new() operator
    // to allocate memory for matrix elements.  Using malloc(3) can
    // save time because it can allocate all required memory in a
    // single call.  However, malloc(3) only works for template
    // arguments that are C++ basic types because basic types are no
    // classes and thus need not be initialised by a constructor.  To
    // use malloc(3) #define CC_CWTMTX_ASSUME_BASIC_TYPES in all
    // source files that use CWTMatrix templates before you #include
    // them.

    // Allocate space for row pointers and the rows themselves using
    // ANSI C malloc(3) function.
    m_rgrow = reinterpret_cast<T **>(malloc(m_crow*sizeof(T *)
					    + m_crow*m_ccol*sizeof(T)));
    T *ptTmp = reinterpret_cast<T *>(&(m_rgrow[m_crow]));
#else
    // Allocate space for row pointers and the rows themselves using
    // ANSI C++ new() operator.
    m_rgrow = new T*[m_crow];
    T *ptTmp = new T[m_crow*m_ccol];
#endif

    // make row pointers point to start of each row
    for (unsigned irow = 0; irow < m_crow; ++irow)
      {
	m_rgrow[irow] = &(ptTmp[irow*m_ccol]);
      }

    m_nMatStatus = N_ALLOCATED;
  }

  // maps a matrix into another matrix, deallocates first if neccesary
  // allocates space on the free store for a matrix, deallocates first
  // if neccesary
  template < class T >
  void CWTMatrix<T>::MapInto(const CWTMatrix<T> &mat,
			     unsigned irowStart,
			     unsigned icolStart,
			     unsigned irowEnd,
			     unsigned icolEnd )
  {
    if (m_nMatStatus != N_NOTALLOCATED)
      {
	Deallocate();
      }

    // calculate columns
    m_crow = irowEnd - irowStart + 1;

    // calculate rows
    m_ccol = icolEnd - icolStart + 1;

    // allocate space for row pointers
#ifdef CC_CWTMTX_ASSUME_BASIC_TYPES
    m_rgrow = reinterpret_cast<T **>(malloc(m_crow*sizeof(T *)));
#else
    m_rgrow = new T*[m_crow];
#endif

    for (unsigned irow = 0; irow < m_crow; ++irow)
      {
	// get values for row pointers
	m_rgrow[irow] = &mat.m_rgrow[irow + irowStart][icolStart];