matrix.h

一个方便进行矩阵运算的C＋＋函数包

   return *this;
}

// unary negation operator
MAT_TEMPLATE inline matrixT
matrixT::operator - () _NO_THROW
{
   matrixT temp(Row,Col);

   for (size_t i=0; i < Row; i++)
      for (size_t j=0; j < Col; j++)
         temp.Val[i][j] = - Val[i][j];

   return temp;
}

// binary addition operator
MAT_TEMPLATE inline matrixT
operator + (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
{
   if (m1.Row != m2.Row || m1.Col != m2.Col)
      REPORT_ERROR( "matrixT::operator+: Inconsistent matrix size in addition!");

   matrixT temp(m1.Row,m1.Col);

   for (size_t i=0; i < m1.Row; i++)
      for (size_t j=0; j < m1.Col; j++)
         temp.Val[i][j] = m1.Val[i][j] + m2.Val[i][j];

   return temp;
}

// binary subtraction operator
MAT_TEMPLATE inline matrixT
operator - (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
{
   if (m1.Row != m2.Row || m1.Col != m2.Col)
     REPORT_ERROR( "matrixT::operator-: Inconsistent matrix size in subtraction!");

   matrixT temp(m1.Row,m1.Col);

   for (size_t i=0; i < m1.Row; i++)
      for (size_t j=0; j < m1.Col; j++)
         temp.Val[i][j] = m1.Val[i][j] - m2.Val[i][j];

   return temp;
}


// binary scalar multiplication operator
MAT_TEMPLATE inline matrixT
operator * (const matrixT& m, const T& no) _NO_THROW
{
   matrixT temp(m.Row,m.Col);

   for (size_t i=0; i < m.Row; i++)
      for (size_t j=0; j < m.Col; j++)
         temp.Val[i][j] = no * m.Val[i][j];

   return temp;
}


// binary matrix multiplication operator
MAT_TEMPLATE inline matrixT
operator * (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
{
   if (m1.Col != m2.Row)
      REPORT_ERROR( "matrixT::operator*: Inconsistent matrix size in multiplication!");

   matrixT temp(m1.Row,m2.Col);

   for (size_t i=0; i < m1.Row; i++)
      for (size_t j=0; j < m2.Col; j++)
      {
         temp.Val[i][j] = T(0);
         for (size_t k=0; k < m1.Col; k++)
            temp.Val[i][j] += m1.Val[i][k] * m2.Val[k][j];
      }
   return temp;
}

// binary power operator
MAT_TEMPLATE inline matrixT
operator ^ (const matrixT& m, const size_t& pow) _THROW_MATRIX_ERROR
{
   matrixT temp(m);

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

   return temp;
}


// unary transpose operator
MAT_TEMPLATE inline matrixT
operator  ~ (const matrixT& m) _NO_THROW
{
   matrixT temp(m.Col,m.Row);

   for (size_t i=0; i < m.Row; i++)
      for (size_t j=0; j < m.Col; j++)
	 temp.Val[j][i] = m.Val[i][j];

   return temp;
}


// unary inversion operator
MAT_TEMPLATE inline matrixT
operator ! (matrixT m) _THROW_MATRIX_ERROR
{
   size_t i,j,k;
   T a1,a2,*rowptr;

   if (m.Row != m.Col)
      REPORT_ERROR( "matrixT::operator!: Inversion of a non-square matrix");

   matrixT temp(m.Row,m.Col);

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

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


// solve simultaneous equation
MAT_TEMPLATE inline matrixT
matrixT::Solve (const matrixT& v) const _THROW_MATRIX_ERROR
{
   size_t i,j,k;
   T a1;

   if (!(Row == Col && Col == v.Row))
      REPORT_ERROR( "matrixT::Solve():Inconsistent matrices!");

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

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

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


// set zero to all elements of this matrix
MAT_TEMPLATE inline void
matrixT::Null (const size_t& row, const size_t& col) _NO_THROW
{
   if (row != Row || col != Col)
      realloc( row,col);

   for (size_t i=0; i < Row; i++)
      for (size_t j=0; j < Col; j++)
         Val[i][j] = T(0);
   return;
}

// set zero to all elements of this matrix
MAT_TEMPLATE inline void
matrixT::Null() _NO_THROW
{
   for (size_t i=0; i < Row; i++)
      for (size_t j=0; j < Col; j++)
 	 Val[i][j] = T(0);
   return;
}

// set this matrix to unity
MAT_TEMPLATE inline void
matrixT::Unit (const size_t& row) _NO_THROW
{
   if (row != Row || row != Col)
      realloc( row, row);

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

// set this matrix to unity
MAT_TEMPLATE inline void
matrixT::Unit () _NO_THROW
{
   size_t row = min(Row,Col);
   Row = Col = row;

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


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

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


// calculate the determinant of a matrix
MAT_TEMPLATE inline T
matrixT::Det () _THROW_MATRIX_ERROR
{
   size_t i,j,k;
   T piv,detVal = T(1);

   if (Row != Col)
      REPORT_ERROR( "matrixT::Det(): Determinant a non-square matrix!");

   matrixT temp(*this);

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

// calculate the norm of a matrix
MAT_TEMPLATE inline T
matrixT::Norm () _NO_THROW
{
   T retVal = T(0);

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

   return retVal;
}


// calculate the condition number of a matrix
MAT_TEMPLATE inline T
matrixT::Cond () _NO_THROW
{
   matrixT inv(Row,Col);

   inv = ! (*this);
   T retVal = Norm() * inv.Norm();

   return retVal;
}


// calculate the cofactor of a matrix for a given element
MAT_TEMPLATE inline T
matrixT::Cofact (size_t row, size_t col) _THROW_MATRIX_ERROR
{
   size_t i,i1,j,j1;

   if (Row != Col)
      REPORT_ERROR( "matrixT::Cofact(): Cofactor of a non-square matrix!");

   if (row > Row || col > Col)
      REPORT_ERROR( "matrixT::Cofact(): Index out of range!");

   matrixT temp (Row-1,Col-1);

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

   return cof;
}


// calculate adjoin of a matrix
MAT_TEMPLATE inline matrixT
matrixT::Adj () _THROW_MATRIX_ERROR
{
   if (Row != Col)
      REPORT_ERROR( "matrixT::Adj(): Adjoin of a non-square matrix.");

   matrixT temp(Row,Col);

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

   temp = ~temp;
   return temp;
}

// Determine if the matrix is singular
MAT_TEMPLATE inline bool
matrixT::IsSingular () _NO_THROW
{
   if (Row != Col)
      return false;
   return (Det() == T(0));
}

// Determine if the matrix is diagonal
MAT_TEMPLATE inline bool
matrixT::IsDiagonal () _NO_THROW
{
   if (Row != Col)
      return false;
   for (size_t i=0; i < Row; i++)
     for (size_t j=0; j < Col; j++)
        if (i != j && Val[i][j] != T(0))
          return false;
   return true;
}

// Determine if the matrix is scalar
MAT_TEMPLATE inline bool
matrixT::IsScalar () _NO_THROW
{
   if (!IsDiagonal())
     return false;
   T v = Val[0][0];
   for (size_t i=1; i < Row; i++)
     if (Val[i][i] != v)
        return false;
   return true;
}

// Determine if the matrix is a unit matrix
MAT_TEMPLATE inline bool
matrixT::IsUnit () _NO_THROW
{
   if (IsScalar() && Val[0][0] == T(1))
     return true;
   return false;
}

// Determine if this is a null matrix
MAT_TEMPLATE inline bool
matrixT::IsNull () _NO_THROW
{
   for (size_t i=0; i < Row; i++)
      for (size_t j=0; j < Col; j++)
         if (Val[i][j] != T(0))
            return false;
   return true;
}

// Determine if the matrix is symmetric
MAT_TEMPLATE inline bool
matrixT::IsSymmetric () _NO_THROW
{
   if (Row != Col)
      return false;
   for (size_t i=0; i < Row; i++)
      for (size_t j=0; j < Col; j++)
         if (Val[i][j] != Val[j][i])
            return false;
   return true;
}
           
// Determine if the matrix is skew-symmetric
MAT_TEMPLATE inline bool
matrixT::IsSkewSymmetric () _NO_THROW
{
   if (Row != Col)
      return false;
   for (size_t i=0; i < Row; i++)
      for (size_t j=0; j < Col; j++)
         if (Val[i][j] != -Val[j][i])
            return false;
   return true;
}
   
// Determine if the matrix is upper triangular
MAT_TEMPLATE inline bool
matrixT::IsUpperTiangular () _NO_THROW
{
   if (Row != Col)
      return false;
   for (size_t i=1; i < Row; i++)
      for (size_t j=0; j < i-1; j++)
         if (Val[i][j] != T(0))
            return false;
   return true;
}

// Determine if the matrix is lower triangular
MAT_TEMPLATE inline bool
matrixT::IsLowerTiangular () _NO_THROW
{
   if (Row != Col)
      return false;

   for (size_t j=1; j < Col; j++)
      for (size_t i=0; i < j-1; i++)
         if (Val[i][j] != T(0))
            return false;

   return true;
}

#ifndef _NO_NAMESPACE
} 
#endif

#endif //__STD_MATRIX_H