matrix.cpp

快速fft变换的c实现

// Matrix.cpp: implementation of the CMatrix class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "FastRBF.h"
#include "Matrix.h"
#include <math.h>

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CMatrix::CMatrix(size_t row, size_t col)
{
	_m = new base_mat( row, col, 0);
}

CMatrix::~CMatrix()
{
	if (--_m->Refcnt == 0)
		delete _m;
}

// copy constructor
//从另一个矩阵构造一新矩阵
CMatrix::CMatrix(const CMatrix &m)
{
	_m = m._m;
    _m->Refcnt++;
}

// Internal copy constructor
void CMatrix::clone ()
{
    _m->Refcnt--;
    _m = new base_mat( _m->Row, _m->Col, _m->Val);
}

// assignment operator
//赋值操作
CMatrix& CMatrix::operator = (const CMatrix& m)
{
    m._m->Refcnt++;
    if (--_m->Refcnt == 0) delete _m;
    _m = m._m;
    return *this;
}

//  reallocation method
void CMatrix::realloc (size_t row, size_t col)
{
   if (row == _m->RowSiz && col == _m->ColSiz)
   {
      _m->Row = _m->RowSiz;
      _m->Col = _m->ColSiz;
      return;
   }

   base_mat *m1 = new base_mat( row, col, NULL);
   size_t colSize = min(_m->Col,col) * sizeof(double);
   size_t minRow = min(_m->Row,row);

   for (size_t i=0; i < minRow; i++)
      memcpy( m1->Val[i], _m->Val[i], colSize);

   if (--_m->Refcnt == 0) 
       delete _m;
   _m = m1;

   return;
}

// public method for resizing matrix
void CMatrix::SetSize (size_t row, size_t col) 
{
   size_t i,j;
   size_t oldRow = _m->Row;
   size_t oldCol = _m->Col;

   if (row != _m->RowSiz || col != _m->ColSiz)
      realloc( row, col);

   for (i=oldRow; i < row; i++)
      for (j=0; j < col; j++)
	 _m->Val[i][j] = double(0);

   for (i=0; i < row; i++)                      
      for (j=oldCol; j < col; j++)
	 _m->Val[i][j] = double(0);

   return;
}

// subscript operator to get/set individual elements
//下标引用

double& CMatrix::operator () (size_t row, size_t col)
{
//   if (row >= _m->Row || col >= _m->Col)
//      REPORT_ERROR( "CMatrix::operator(): Index out of range!");
   if (_m->Refcnt > 1) clone();
   return _m->Val[row][col];
}

// subscript operator to get/set individual elements
//下标引用
double CMatrix::operator () (size_t row, size_t col) const
{
 //  if (row >= _m->Row || col >= _m->Col)
//      REPORT_ERROR( "CMatrix::operator(): Index out of range!");
   return _m->Val[row][col];
}

// input stream function
//输入


// output stream function
//输出



// logical equal-to operator
bool CMatrix::operator == (const CMatrix& m) const
{
   if (RowNo() != m.RowNo() || ColNo() != m.ColNo())
      return false;

   for (size_t i=0; i < RowNo(); i++)
      for (size_t j=0; j < ColNo(); j++)
	      if ((i,j) != m(i,j))
	         return false;

   return true;
}


// logical no-equal-to operator
bool CMatrix::operator != (const CMatrix& m) const
{
    return (*this == m) ? false : true;
}

// combined addition and assignment operator
CMatrix& CMatrix::operator += (const CMatrix& m) 
{
   if (_m->Row != m._m->Row || _m->Col != m._m->Col)
//      REPORT_ERROR( "CMatrix::operator+= : Inconsistent matrix sizes in addition!");
   if (_m->Refcnt > 1) clone();
   for (size_t i=0; i < m._m->Row; i++)
      for (size_t j=0; j < m._m->Col; j++)
	 _m->Val[i][j] += m._m->Val[i][j];
   return *this;
}

// combined subtraction and assignment operator
CMatrix& CMatrix::operator -= (const CMatrix& m)
{
   if (_m->Row != m._m->Row || _m->Col != m._m->Col)
//      REPORT_ERROR( "CMatrix::operator-= : Inconsistent matrix sizes in subtraction!");
   
   if (_m->Refcnt > 1) clone();
   
   for (size_t i=0; i < m._m->Row; i++)
      for (size_t j=0; j < m._m->Col; j++)
	 _m->Val[i][j] -= m._m->Val[i][j];
   return *this;
}

// combined scalar multiplication and assignment operator
CMatrix& CMatrix::operator *= (const double& c)
{
    if (_m->Refcnt > 1) clone();
    for (size_t i=0; i < _m->Row; i++)
	for (size_t j=0; j < _m->Col; j++)
	    _m->Val[i][j] *= c;
    return *this;
}

// combined matrix multiplication and assignment operator
CMatrix& CMatrix::operator *= (const CMatrix& m)
{
//   if (_m->Col != m._m->Row)
//      REPORT_ERROR( "CMatrix::operator*= : Inconsistent matrix sizes in multiplication!");
   CMatrix temp(_m->Row,m._m->Col);
   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < m._m->Col; j++)
      {
         temp._m->Val[i][j] = double(0);
         for (size_t k=0; k < _m->Col; k++)
            temp._m->Val[i][j] += _m->Val[i][k] * m._m->Val[k][j];
      }
   *this = temp;

   return *this;
}

// combined scalar division and assignment operator
CMatrix& CMatrix::operator /= (const double& c)
{
    if (_m->Refcnt > 1) clone();
    for (size_t i=0; i < _m->Row; i++)
	for (size_t j=0; j < _m->Col; j++)
	    _m->Val[i][j] /= c;

    return *this;
}

// combined power and assignment operator
CMatrix& CMatrix::operator ^= (const size_t& pow) 
{
	CMatrix temp(*this);

	for (size_t i=2; i <= pow; i++)
      *this *= temp;

	return *this;
}

CMatrix CMatrix::operator ! ()
{
	return Reverse();
}

CMatrix CMatrix::operator + (CMatrix &m)
{
	if (m.RowNo() != RowNo())
	{
		TRACE("Unmatched row size of matrixes!");
		return *this;
	}
	if (m.ColNo() != ColNo())
	{
		TRACE("Unmatched col size of matrixes!");
		return *this;
	}
	for (size_t i=0; i<RowNo(); i++)
	{
		for (size_t j=0; j<ColNo(); j++)
			_m->Val[i][j] += m(i,j);
	}
	return *this;
}

CMatrix CMatrix::operator - (CMatrix &m)
{
	if (m.RowNo() != RowNo())
	{
		TRACE("Unmatched row size of matrixes!");
		return *this;
	}
	if (m.ColNo() != ColNo())
	{
		TRACE("Unmatched col size of matrixes!");
		return *this;
	}
	for (size_t i=0; i<RowNo(); i++)
	{
		for (size_t j=0; j<ColNo(); j++)
			_m->Val[i][j] -= m(i,j);
	}
	return *this;
}

CMatrix CMatrix::operator ~ ()
{
	CMatrix temp(ColNo(),RowNo());
	for (size_t i=0; i < RowNo(); i++)
	{
		for (size_t j=0; j < ColNo(); j++)
		{
			temp(j,i) = (i,j);
		}
	}
   return temp;
}

// reverse function
CMatrix  CMatrix::Reverse () 
{
   size_t i,j,k;
   double a1,a2,*rowptr;
   CMatrix  tempthis=*this;

 //  if (_m->Row != _m->Col)
 //     REPORT_ERROR( "CMatrix::operator!: Inversion of a non-square matrix");

   CMatrix temp(_m->Row,_m->Col);
   if (_m->Refcnt > 1) clone();


   temp.Unit();
   for (k=0; k < _m->Row; k++)
   {
      int indx = pivot(k);
//     if (indx == -1)
//	      REPORT_ERROR( "CMatrix::operator!: Inversion of a singular matrix");

      if (indx != 0)
      {
	      rowptr = temp._m->Val[k];
	      temp._m->Val[k] = temp._m->Val[indx];
	      temp._m->Val[indx] = rowptr;
      }
      a1 = _m->Val[k][k];
      for (j=0; j < _m->Col; j++)
      {
	      _m->Val[k][j] /= a1;
	      temp._m->Val[k][j] /= a1;
      }
      for (i=0; i < _m->Row; i++)
	   if (i != k)
	   {
	      a2 = _m->Val[i][k];
	      for (j=0; j < _m->Col; j++)
	      {
	         _m->Val[i][j] -= a2 * _m->Val[k][j];
	         temp._m->Val[i][j] -= a2 * temp._m->Val[k][j];
	      }
	   }
   }
   *this=tempthis;
   return temp;
}

// solve simultaneous equation
CMatrix CMatrix::Solve (const CMatrix& v) const 
{
   size_t i,j,k;
   double a1;

 //  if (!(_m->Row == _m->Col && _m->Col == v._m->Row))
 //     REPORT_ERROR( "CMatrix::Solve():Inconsistent matrices!");

   CMatrix temp(_m->Row,_m->Col+v._m->Col);
   for (i=0; i < _m->Row; i++)
   {
      for (j=0; j < _m->Col; j++)
	 temp._m->Val[i][j] = _m->Val[i][j];
      for (k=0; k < v._m->Col; k++)
	 temp._m->Val[i][_m->Col+k] = v._m->Val[i][k];
   }
   for (k=0; k < _m->Row; k++)
   {
      int indx = temp.pivot(k);
//      if (indx == -1)
//	 REPORT_ERROR( "CMatrix::Solve(): Singular matrix!");

      a1 = temp._m->Val[k][k];
      for (j=k; j < temp._m->Col; j++)
	 temp._m->Val[k][j] /= a1;

      for (i=k+1; i < _m->Row; i++)
      {
	 a1 = temp._m->Val[i][k];
	 for (j=k; j < temp._m->Col; j++)
	   temp._m->Val[i][j] -= a1 * temp._m->Val[k][j];
      }
   }
   CMatrix s(v._m->Row,v._m->Col);
   for (k=0; k < v._m->Col; k++)
      for (int m=int(_m->Row)-1; m >= 0; m--)
      {
	 s._m->Val[m][k] = temp._m->Val[m][_m->Col+k];
	 for (j=m+1; j < _m->Col; j++)
	    s._m->Val[m][k] -= temp._m->Val[m][j] * s._m->Val[j][k];
      }
   return s;
}

// set zero to all elements of this matrix
//置零
void  CMatrix::Null (const size_t& row, const size_t& col) 
{
    if (row != _m->Row || col != _m->Col)
	realloc( row,col);

    if (_m->Refcnt > 1) 
	clone();

    for (size_t i=0; i < _m->Row; i++)
	for (size_t j=0; j < _m->Col; j++)
	    _m->Val[i][j] = double(0);
    return;
}

// set zero to all elements of this matrix
//全部元素置零
void CMatrix::Null() 
{
    if (_m->Refcnt > 1) clone();   
    for (size_t i=0; i < _m->Row; i++)
	for (size_t j=0; j < _m->Col; j++)
		_m->Val[i][j] = double(0);
    return;
}

// set this matrix to unity
//生成单位矩阵
void CMatrix::Unit (const size_t& row) 
{
    if (row != _m->Row || row != _m->Col)
	realloc( row, row);
	
    if (_m->Refcnt > 1) 
	clone();

    for (size_t i=0; i < _m->Row; i++)
	for (size_t j=0; j < _m->Col; j++)
	    _m->Val[i][j] = i == j ? double(1) : double(0);
    return;
}

// set this matrix to unity
void CMatrix::Unit () 
{
    if (_m->Refcnt > 1) clone();   
    size_t row = min(_m->Row,_m->Col);
    _m->Row = _m->Col = row;

    for (size_t i=0; i < _m->Row; i++)
	for (size_t j=0; j < _m->Col; j++)
	    _m->Val[i][j] = i == j ? double(1) : double(0);
    return;
}

// set this matrix to zero matrix
void CMatrix::Zero () 
{
    if (_m->Refcnt > 1) clone();   
   
    for (size_t i=0; i < _m->Row; i++)
	for (size_t j=0; j < _m->Col; j++)
	    _m->Val[i][j] = 0;
    return;
}

// private partial pivoting method
int CMatrix::pivot (size_t row)
{
  int k = int(row);
  double amax,temp;

  amax = -1;
  for (size_t i=row; i < _m->Row; i++)
    if ( (temp = fabs( _m->Val[i][row])) > amax && temp != 0.0)
     {
       amax = temp;
       k = i;
     }
  if (_m->Val[k][row] == double(0))
     return -1;
  if (k != int(row))
  {
     double* rowptr = _m->Val[k];
     _m->Val[k] = _m->Val[row];
     _m->Val[row] = rowptr;
     return k;
  }
  return 0;
}

// calculate the determinant of a matrix
double CMatrix::Det () const 
{
   size_t i,j,k;
   double piv,detVal = double(1);

//   if (_m->Row != _m->Col)
//      REPORT_ERROR( "CMatrix::Det(): Determinant a non-square matrix!");
   
   CMatrix temp(*this);
   if (temp._m->Refcnt > 1) temp.clone();

   for (k=0; k < _m->Row; k++)
   {
      int indx = temp.pivot(k);
      if (indx == -1)
	 return 0;
      if (indx != 0)
	 detVal = - detVal;
      detVal = detVal * temp._m->Val[k][k];
      for (i=k+1; i < _m->Row; i++)
      {
	 piv = temp._m->Val[i][k] / temp._m->Val[k][k];
	 for (j=k+1; j < _m->Row; j++)
	    temp._m->Val[i][j] -= piv * temp._m->Val[k][j];
      }
   }
   return detVal;
}

// calculate the norm of a matrix
double CMatrix::Norm () 
{
   double retVal = double(0);

   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < _m->Col; j++)
	 retVal += _m->Val[i][j] * _m->Val[i][j];
   retVal = sqrt( retVal);

   return retVal;
}

// calculate the condition number of a matrix
double CMatrix::Cond () 
{
   CMatrix inv = ! (*this);
   return (Norm() * inv.Norm());
}

// calculate the cofactor of a matrix for a given element
double CMatrix::Cofact (size_t row, size_t col) 
{
   size_t i,i1,j,j1;

//   if (_m->Row != _m->Col)
//      REPORT_ERROR( "CMatrix::Cofact(): Cofactor of a non-square matrix!");

//   if (row > _m->Row || col > _m->Col)
//      REPORT_ERROR( "CMatrix::Cofact(): Index out of range!");

   CMatrix temp (_m->Row-1,_m->Col-1);

   for (i=i1=0; i < _m->Row; i++)
   {
      if (i == row)
	continue;
      for (j=j1=0; j < _m->Col; j++)
      {
	 if (j == col)
	    continue;
	 temp._m->Val[i1][j1] = _m->Val[i][j];
	 j1++;
      }
      i1++;
   }
   double  cof = temp.Det();
   if ((row+col)%2 == 1)
      cof = -cof;

   return cof;
}

// calculate adjoin of a matrix
CMatrix CMatrix::Adj () 
{
//   if (_m->Row != _m->Col)
//      REPORT_ERROR( "CMatrix::Adj(): Adjoin of a non-square matrix.");
   CMatrix temp(_m->Row,_m->Col);

   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < _m->Col; j++)
	  temp._m->Val[j][i] = Cofact(i,j);
   return temp;
}

// Determine if the matrix is singular
bool  CMatrix::IsSingular () 
{
   if (_m->Row != _m->Col)
      return false;
   return (Det() == double(0));
}

// Determine if the matrix is diagonal
bool CMatrix::IsDiagonal () 
{
   if (_m->Row != _m->Col)
      return false;
   for (size_t i=0; i < _m->Row; i++)
     for (size_t j=0; j < _m->Col; j++)
	if (i != j && _m->Val[i][j] != double(0))
	  return false;
   return true;
}

// Determine if the matrix is scalar
 bool CMatrix::IsScalar () 
{
   if (!IsDiagonal())
     return false;
   double v = _m->Val[0][0];
   for (size_t i=1; i < _m->Row; i++)
     if (_m->Val[i][i] != v)
	return false;
   return true;
}

// Determine if the matrix is a unit matrix
 bool CMatrix::IsUnit () 
{
   if (IsScalar() && _m->Val[0][0] == double(1))
     return true;
   return false;
}

// Determine if this is a null matrix
 bool CMatrix::IsNull () 
{
   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < _m->Col; j++)
	 if (_m->Val[i][j] != double(0))
	    return false;
   return true;
}

// Determine if the matrix is symmetric
 bool CMatrix::IsSymmetric () 
{
   if (_m->Row != _m->Col)