smatrix.h

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

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

// $Id: smatrix.h,v 1.26 2002/01/15 01:57:23 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

#ifndef IG_MATRIX_H
#include "matrix.h"
#endif

namespace CwMtx
{
  // prefix smat
  template < class T = double >
  class CWTSquareMatrix : public CWTMatrix<T>
  {
  public:
    CWTSquareMatrix(): CWTMatrix<T>() {};
    CWTSquareMatrix(unsigned crowInit): CWTMatrix<T>( crowInit, crowInit) {};
    CWTSquareMatrix(const CWTMatrix<T> &mat): CWTMatrix<T>(mat) {};
    CWTSquareMatrix(const CWTSquareMatrix &smat): CWTMatrix<T>(smat) {};
    CWTSquareMatrix(const CWTMatrix<T> &, unsigned, unsigned, unsigned);
    CWTSquareMatrix(const CWTSquareMatrix<T> &, unsigned, unsigned, unsigned);

    ~CWTSquareMatrix() {};

    void Dimension(unsigned crowInit)
    {
      CWTMatrix<T>::Dimension(crowInit, crowInit);
    }

    void MapInto(const CWTSquareMatrix &, unsigned, unsigned, unsigned);

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

    // not inherited
    CWTSquareMatrix & operator =(const CWTSquareMatrix &smat);
    CWTSquareMatrix & operator +=(const CWTSquareMatrix &smat);
    CWTSquareMatrix & operator -=(const CWTSquareMatrix &smat);
    CWTSquareMatrix & operator *=(const T &value);
    CWTSquareMatrix & operator *=(const CWTSquareMatrix &);
    CWTSquareMatrix & operator /=(const T &value);
    CWTSquareMatrix & operator /=(const CWTSquareMatrix &);

    // stores the adjoint of argument in this
    void StoreAdjoint(const CWTSquareMatrix &);
    // stores the inverse of argument in this
    void StoreInverse(const CWTSquareMatrix &);
    // makes this its own adjoint
    void MakeAdjoint();
    // makes this its own inverse
    void MakeInverse();
    // makes this a unity matrix
    void MakeUnity();
  };

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

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

  // Unity square matrix
  template <class T, unsigned crow>
  class CWTUnity< CWTSMat<CWTSquareMatrix<T>, crow> >:
    public CWTSMat<CWTSquareMatrix<T>, crow>
  {
  public:
    CWTUnity() { MakeUnity(); }
  };

  // Zero matrix.
  template <class T, unsigned crow>
  class CWTZero< CWTSMat<CWTSquareMatrix<T>, crow> >:
    public CWTSMat<CWTSquareMatrix<T>, crow>
  {
  public:
    CWTZero() { Fill(CWTZero<T>()); }
  };

  //
  // Constructors
  //

  template < class T >
  inline CWTSquareMatrix<T>::CWTSquareMatrix(const CWTMatrix<T> &mat,
					     unsigned irowStart,
					     unsigned icolStart,
					     unsigned irowEnd)
    :
    CWTMatrix<T>(mat,
		 irowStart,
		 icolStart,
		 irowEnd,
		 icolStart + irowEnd - irowStart)
  {
  }

  template < class T >
  inline CWTSquareMatrix<T>::CWTSquareMatrix(const CWTSquareMatrix<T> &smat,
					     unsigned irowStart,
					     unsigned icolStart,
					     unsigned irowEnd)
    :
    CWTMatrix<T>(smat,
		 irowStart,
		 icolStart,
		 irowEnd,
		 icolStart + irowEnd - irowStart)
  {
  }

  //
  // User Methods
  //

  template < class T >
  inline void CWTSquareMatrix<T>::MapInto(const CWTSquareMatrix<T> &smat,
					  unsigned irowStart,
					  unsigned icolStart,
					  unsigned irowEnd)
  {
    CWTMatrix<T>::MapInto(smat,
			  irowStart,
			  icolStart,
			  irowEnd,
			  icolStart + irowEnd - irowStart);
  }

  // not inherited
  template < class T >
  inline CWTSquareMatrix<T> &
  CWTSquareMatrix<T>::operator =(const CWTSquareMatrix<T> &smat)
  {
    return static_cast<CWTSquareMatrix<T> &>(CWTMatrix<T>::operator=(smat));
  }

  template < class T >
  inline CWTSquareMatrix<T> &
  CWTSquareMatrix<T>::operator +=(const CWTSquareMatrix<T> &smat)
  {
    return static_cast<CWTSquareMatrix<T> &>(CWTMatrix<T>::operator+=(smat));
  }

  template < class T >
  inline CWTSquareMatrix<T> &
  CWTSquareMatrix<T>::operator -=(const CWTSquareMatrix &smat)
  {
    return static_cast<CWTSquareMatrix<T> &>(CWTMatrix<T>::operator-=(smat));
  }

  template < class T >
  inline CWTSquareMatrix<T> &
  CWTSquareMatrix<T>::operator *=(const T &value)
  {
    return static_cast<CWTSquareMatrix<T> &>(CWTMatrix<T>::operator*=(value));
  }

  template < class T >
  inline CWTSquareMatrix<T> &
  CWTSquareMatrix<T>::operator *=(const CWTSquareMatrix<T> &smat)
  {
    return (*this) = (*this)*smat;
  }

  template < class T >
  inline CWTSquareMatrix<T> &
  CWTSquareMatrix<T>::operator /=(const T &value)
  {
    return (*this) *= static_cast<const T &>(CWTUnity<T>()/value);
  }

  template < class T >
  inline CWTSquareMatrix<T> &
  CWTSquareMatrix<T>::operator /=(const CWTSquareMatrix<T> &smat)
  {
    return (*this) = (*this)/smat;
  }

  template < class T >
  CWTSquareMatrix<T>
  CWTSquareMatrix<T>::operator +(const CWTSquareMatrix<T> &smat) const
  {
    return CWTSquareMatrix(*this) += smat;
  }

  template < class T >
  CWTSquareMatrix<T>
  CWTSquareMatrix<T>::operator -(const CWTSquareMatrix<T> &smat) const
  {
    return CWTSquareMatrix(*this) -= smat;
  }

  template < class T >
  CWTSquareMatrix<T> CWTSquareMatrix<T>::operator -() const
  {
    return (*this)*static_cast<const T &>(CWTZero<T>() - CWTUnity<T>());
  }

  template < class T >
  CWTSquareMatrix<T> CWTSquareMatrix<T>::operator *(const T &value) const
  {
    return CWTSquareMatrix(*this) *= value;
  }

  template < class T >
  CWTSquareMatrix<T>
  CWTSquareMatrix<T>::operator *(const CWTSquareMatrix<T> &smat) const
  {
    CWTSquareMatrix smatResult(GetRows());
    smatResult.StoreProduct(*this, smat);
    return smatResult;
  }

  template < class T >
  CWTSquareMatrix<T>
  CWTSquareMatrix<T>::operator /(const CWTSquareMatrix<T> &smat) const
  {
    return (*this)*inv(smat);
  }

  // stores inverse of argument in this
  template < class T >
  void CWTSquareMatrix<T>::StoreInverse(const CWTSquareMatrix<T> &smat)
  {
    // copy input matrix in this
    (*this) = smat;
    MakeInverse();
  }

  // makes this a unity matrix
  template < class T >
  void CWTSquareMatrix<T>::MakeUnity()
  {
    unsigned crow(GetRows()), ccol(GetCols());

    for (unsigned irow = 0; irow < crow; ++irow)
      {
	for (unsigned icol = 0; icol < ccol; ++icol)
	  {
	    if (irow == icol)
	      {
		(*this)[irow][icol] = CWTUnity<T>();
	      }
	    else
	      {
		(*this)[irow][icol] = CWTZero<T>();
	      }
	  }
      }
  }

  // makes this its own adjoint
  template < class T >
  void CWTSquareMatrix<T>::MakeAdjoint()
  {
    // we need a copy of this
    CWTSquareMatrix smatCopy(*this);
    // for easier access to crows
    unsigned crowCopy = smatCopy.GetRows();
    T elemTmp;
    // will eventually contain det(smatCopy)
    T elemDet = CWTUnity<T>();

    // make this a unity matrix
    MakeUnity();

    // start row reduction
    for (unsigned irow = 0; irow < crowCopy; ++irow)
      {
	// if element [irow][irow] is zero, add a row with non-zero
	// element at column irow
	if (smatCopy[irow][irow] == CWTZero<T>())
	  {
	    for (unsigned irowTmp = irow; irowTmp < crowCopy; ++irowTmp)
	      {
		if (smatCopy[irowTmp][irow] != CWTZero<T>())
		  {
		    // has no effect on adj(smat)
		    smatCopy.AddRowToRow(irowTmp, irow);
		    // repeat action on this
		    AddRowToRow(irowTmp, irow);
		    break;
		  };
	      };
	  };

	elemTmp = smatCopy[irow][irow];
	smatCopy.MultiplyRow(irow, CWTUnity<T>()/elemTmp);
	// repeat action on this
	MultiplyRow(irow, CWTUnity<T>()/elemTmp);
	elemDet *= elemTmp;

	for (unsigned irowToAddTo = 0; irowToAddTo < crowCopy; ++irowToAddTo)
	  {
	    if (irowToAddTo != irow)
	      {
		elemTmp = -smatCopy[irowToAddTo][irow];
		smatCopy.AddRowToRow(irow, irowToAddTo, elemTmp);
		// repeat action on this
		AddRowToRow(irow, irowToAddTo, elemTmp);
	      };
	  };
      };

    // this now contains its adjoint
    (*this) *= elemDet;
  }

  // stores adjoint of input matrix in this
  template < class T >
  void CWTSquareMatrix<T>::StoreAdjoint(const CWTSquareMatrix<T> &smat)
  {
    // copy input matrix in this
    (*this) = smat;
    MakeAdjoint();
  }

  // makes this its own inverse
  template < class T >
  void CWTSquareMatrix<T>::MakeInverse()
  {
    // we need a copy of this
    CWTSquareMatrix smatCopy(*this);
    // for easier access to crows
    unsigned crowCopy = smatCopy.GetRows();
    T elemTmp;

    // make this a unity matrix
    MakeUnity();

    // start row reduction
    for (unsigned irow = 0; irow < crowCopy; ++irow)
      {
	// if element [irow][irow] is zero, add a row with non-zero
	// element at column irow
	if (smatCopy[irow][irow] == CWTZero<T>())
	  {
	    for (unsigned irowTmp = irow; irowTmp < crowCopy; ++irowTmp)
	      {
		// has no effect on inv(smat)
		if (smatCopy[irowTmp][irow] != CWTZero<T>())
		  {
		    smatCopy.AddRowToRow(irowTmp, irow);
		    // repeat action on this
		    AddRowToRow(irowTmp, irow);
		    break;
		  };
	      };
	  };

	elemTmp = smatCopy[irow][irow];
	smatCopy.MultiplyRow(irow,
			     static_cast<const T &>(CWTUnity<T>()/elemTmp));
	MultiplyRow(irow, static_cast<const T &>(CWTUnity<T>()/elemTmp));

	for (unsigned irowToAddTo = 0; irowToAddTo < crowCopy; ++irowToAddTo)
	  {
	    if (irowToAddTo != irow)
	      {
		elemTmp = -smatCopy[irowToAddTo][irow];
		smatCopy.AddRowToRow(irow, irowToAddTo, elemTmp);
		// repeat action on this
		AddRowToRow(irow, irowToAddTo, elemTmp);
	      };
	  };
      };

    // this now contains its inverse
  }

  //
  // template functions designed work with the Square Matrix class.
  //

  // SCLR*CWTSquareMatrix
  template < class T >
  inline CWTSquareMatrix<T> operator *(const T &value,
				       const CWTSquareMatrix<T> &smat)
  {
    return smat*value;
  }

  template < class T >
  CWTSquareMatrix<T> transpose(const CWTSquareMatrix<T> &smat)
  {
    CWTSquareMatrix<T> smatResult(smat.GetRows());
    smatResult.StoreTranspose(smat);
    return smatResult;
  }

  // returns the adjoint of a square matrix
  template < class T >
  CWTSquareMatrix<T> adj(const CWTSquareMatrix<T> &smat)
  {
    CWTSquareMatrix<T> smatResult(smat);   // copy input matrix
    smatResult.MakeAdjoint();
    return smatResult;
  }

  // calculates the inverse of a square matrix
  template < class T >
  CWTSquareMatrix<T> inv(const CWTSquareMatrix<T> &smat)
  {
    // copy input matrix
    CWTSquareMatrix<T> smatResult(smat);
    smatResult.MakeInverse();
    return smatResult;
  }

  // calculates the determinant of a square matrix
  template < class T >
  T det(const CWTSquareMatrix<T> &smat)
  {
    // a copy of the input matrix
    CWTSquareMatrix<T> smatCopy(smat);
    unsigned crowCopy = smatCopy.GetRows();

    // start row reduction
    T elemTmp;
    // will eventually contain det(smatCopy)
    T elemDet = CWTUnity<T>();

    for (unsigned irow = 0; irow < crowCopy; ++irow)
      {
	// if element [irow][irow] is zero, add a row with non-zero
	// element at column irow
	if (smatCopy[irow][irow] == CWTZero<T>())
	  {
	    for (unsigned irowTmp = irow; irowTmp < crowCopy; ++irowTmp)
	      {
		// has no effect on inv(smat)
		if (smatCopy[irowTmp][irow] != CWTZero<T>())
		  {
		    smatCopy.AddRowToRow(irowTmp, irow);
		    break;
		  };
	      };
	  };

	elemTmp =  smatCopy[irow][irow];
	elemDet *= elemTmp;

	if (elemDet == CWTZero<T>())
	  {
	    // once elemDet is zero it will stay zero
	    return elemDet;
	  }

	smatCopy.MultiplyRow(irow, CWTUnity<T>()/elemTmp);

	for (unsigned irowToAddTo = 0; irowToAddTo < crowCopy; ++irowToAddTo)
	  {
	    if (irowToAddTo != irow)
	      {
		smatCopy.AddRowToRow(irow,
				     irowToAddTo,
				     -smatCopy[irowToAddTo][irow]);
	      };
	  };
      };

    return elemDet;
  }

  // trace
  template < class T >
  T tr(const CWTSquareMatrix<T> &smat)
  {
    T elemTmp = CWTZero<T>();

    for (unsigned c = 0; c < smat.GetCols(); c++)
      {
	elemTmp += smat[c][c];
      }

    return elemTmp;
  }

}

#endif // IG_SMATRIX_H