matrix.cpp

矩阵类

      MatrixError( "Column", "Unable to copy the source matrix column." );
    }
    return A;
  }

  // Return the row matrix specified by the column index. Returns (ncols x 1).
  Matrix  Matrix::Row(const unsigned row)
  {
    Matrix A;
    if( !MTX_CopyRow(&m_Matrix, row, &A.m_Matrix ) )
    {
      MatrixError( "Column", "Unable to copy the source matrix column." );
    }
    return A;
  }

  // Return the tranpose of the matrix.
  Matrix  Matrix::Transpose()
  {
    Matrix A;
    if( !MTX_Transpose( &m_Matrix, &A.m_Matrix ) )
    {
      MatrixError( "Column", "Unable to transpose the source matrix." );
    }
    return A;
  }

  // Return the tranpose of the matrix.
  Matrix  Matrix::T()
  {
    return Transpose();
  }

  // Return the diagonal of the matrix as a vector.
  Matrix  Matrix::Diagonal()
  {
    Matrix A;
    if( !MTX_Diagonal( &m_Matrix, &A.m_Matrix ) )
    {
      MatrixError( "Diagonal", "Unable to get the diagonal of the source matrix." );
    }
    return A;
  }

  // Return the inverse of the matrix.
  Matrix  Matrix::Inverse()
  {
    Matrix A;
    if( !MTX_Copy( &m_Matrix, &A.m_Matrix ) )
    {
      MatrixError( "Inverse", "Unable to copy the source matrix." );
    }
    if( !MTX_InvertInPlace( &A.m_Matrix ) )
    {
      MatrixError( "Inverse", "Unable to invert the matrix." );
    }
    return A;
  }

  // Return the inverse of the matrix.
  Matrix  Matrix::Inv()// short version
  {
    return Inverse();
  }

  // Return the Fourier Transform of each column of the matrix. Power of two uses FFT, otherwise fast DFT.
  Matrix  Matrix::FFT()
  {
    Matrix A;
    if( !MTX_FFT( &m_Matrix, &A.m_Matrix ) )
    {
      MatrixError( "FFT", "Unable to perform the FFT." );
    }
    return A;
  }

  // Return the inverse Fourier Transform of each column of the matrix. Power of two uses IFFT, otherwise fast IDFT.
  Matrix  Matrix::IFFT()
  {
    Matrix A;
    if( !MTX_IFFT( &m_Matrix, &A.m_Matrix ) )
    {
      MatrixError( "IFFT", "Unable to perform the IFFT." );
    }
    return A;
  }



  /// Get a reference to an element in the matrix to set its value.
  Matrix::Element& Matrix::operator() (unsigned row, unsigned col)
  {
    if( IndexCheck(row,col) )
    {
      m_MatrixElement.m_row = row;
      m_MatrixElement.m_col = col;
    }
    else
    {
      // This code should not be reached!
      m_MatrixElement.m_row = 0;
      m_MatrixElement.m_col = 0;
    }
    return m_MatrixElement; 
  }


  /// Get a reference to an element in the matrix as a column or row vector to set its value.
  Matrix::Element& Matrix::operator() (unsigned index)
  {
    if( IndexCheck( index ) )
    {
      if( m_Matrix.ncols == 1 ) // column array
      {
        m_MatrixElement.m_row = index;
        m_MatrixElement.m_col = 0;
      }
      else if( m_Matrix.nrows == 1 ) // row array
      {
        m_MatrixElement.m_row = 0;
        m_MatrixElement.m_col = index;
      }
      else
      {
        // access the matrix as a singular column array
        m_MatrixElement.m_col = index / m_Matrix.nrows;
        m_MatrixElement.m_row = index - m_MatrixElement.m_col*m_Matrix.nrows;
      }
    }
    else
    {
      // This code should not be reached!
      m_MatrixElement.m_row = 0;
      m_MatrixElement.m_col = 0;
    }

    return m_MatrixElement;
  }


  Matrix::Element::Element(MTX& mtx)
    : m_mtx(mtx), m_row(0), m_col(0)
  {
    // Note that there is no index checking at this level. The 
    // index checking is performed at the Matrix operator() level.
  }

  Matrix::Element::~Element()
  {}

  const double Matrix::Element::real()
  {
    if( m_mtx.isReal )
    {
      return m_mtx.data[m_col][m_row];
    }
    else
    {
      return m_mtx.cplx[m_col][m_row].re;
    }
  }
  
  const double Matrix::Element::imag()
  {
    if( m_mtx.isReal )
    {
      return 0.0;
    }
    else
    {
      return m_mtx.cplx[m_col][m_row].im;
    }
  }

  const Matrix::Element& Matrix::Element::operator= (double v)
  {
    if( m_mtx.isReal )
    {
      m_mtx.data[m_col][m_row] = v;
    }
    else
    {
      m_mtx.cplx[m_col][m_row].re = v;
      m_mtx.cplx[m_col][m_row].im = 0.0;
    }
    return *this;
  }
  
  const Matrix::Element& Matrix::Element::operator= (std::complex<double> v)
  {
    if( m_mtx.isReal )
    {
      if( v.imag() != 0.0 )
      {
        // This is on the fly conversion to a complex matrix!
        if( !MTX_ConvertRealToComplex( &m_mtx ) )
        {
          MTX_Free( &m_mtx );
          Matrix::StaticMatrixError( "Element::operator=", "Unable to convert matrix from real to complex" );
          return *this;
        }
      }
      else
      {
        m_mtx.data[m_col][m_row] = v.real();
        return *this;
      }
    }
    m_mtx.cplx[m_col][m_row].re = v.real();
    m_mtx.cplx[m_col][m_row].im = v.imag();
    return *this;
  }

  const Matrix::Element& Matrix::Element::operator= (Element v)
  {
    if( v.m_mtx.isReal )
    {
      if( m_mtx.isReal )
      {
        m_mtx.data[m_col][m_row] = v.m_mtx.data[v.m_col][v.m_row];
      }
      else
      {
        m_mtx.cplx[m_col][m_row].re = v.m_mtx.data[v.m_col][v.m_row];
        m_mtx.cplx[m_col][m_row].im = 0.0;
      }
    }
    else
    {
      if( m_mtx.isReal )
      {
        // This is on the fly conversion to a complex matrix!
        if( !MTX_ConvertRealToComplex( &m_mtx ) )
        {
          MTX_Free( &m_mtx );
          Matrix::StaticMatrixError( "Element::operator=", "Unable to convert matrix from real to complex" );
          return *this;
        }        
      }
      
      m_mtx.cplx[m_col][m_row].re = v.m_mtx.cplx[v.m_col][v.m_row].re;
      m_mtx.cplx[m_col][m_row].im = v.m_mtx.cplx[v.m_col][v.m_row].im;      
    }
    return *this;
  }


  Matrix::Element::operator const std::complex<double>() const
  {
    if( m_mtx.isReal )
    {
      std::complex<double> v( m_mtx.data[m_col][m_row], 0.0 );
      return v;
    }
    else
    {
      std::complex<double> v( m_mtx.cplx[m_col][m_row].re, m_mtx.cplx[m_col][m_row].im );
      return v;
    }
  }


  void Matrix::Element::operator+= (const double scalar)
  {
    if( m_mtx.isReal )
    {
      m_mtx.data[m_col][m_row] += scalar;
    }
    else
    {
      m_mtx.cplx[m_col][m_row].re += scalar;
    }
  }

  void Matrix::Element::operator+= (const std::complex<double>& v)
  {
    if( m_mtx.isReal )
    {
      if( v.imag() != 0.0 )
      {
        // This is on the fly conversion to a complex matrix!
        if( !MTX_ConvertRealToComplex( &m_mtx ) )
        {
          MTX_Free( &m_mtx );
          Matrix::StaticMatrixError( "Element::operator=", "Unable to convert matrix from real to complex" );
          return;
        }
      }
      else
      {
        m_mtx.data[m_col][m_row] += v.real();
        return;
      }
    }
    
    m_mtx.cplx[m_col][m_row].re += v.real();
    m_mtx.cplx[m_col][m_row].im += v.imag();    
  }
  
  void Matrix::Element::operator+= (const Element& v)
  {
    std::complex<double> cplx = (const std::complex<double>)v;
    *this += cplx;
  }


  void Matrix::Element::operator-= (const double scalar)
  {
    if( m_mtx.isReal )
    {
      m_mtx.data[m_col][m_row] -= scalar;
    }
    else
    {
      m_mtx.cplx[m_col][m_row].re -= scalar;
    }
  }

  void Matrix::Element::operator-= (const std::complex<double>& v)
  {
    if( m_mtx.isReal )
    {
      if( v.imag() != 0.0 )
      {
        // This is on the fly conversion to a complex matrix!
        if( !MTX_ConvertRealToComplex( &m_mtx ) )
        {
          MTX_Free( &m_mtx );
          Matrix::StaticMatrixError( "Element::operator=", "Unable to convert matrix from real to complex" );
          return;
        }
      }
      else
      {
        m_mtx.data[m_col][m_row] -= v.real();
        return;
      }
    }
    
    m_mtx.cplx[m_col][m_row].re -= v.real();
    m_mtx.cplx[m_col][m_row].im -= v.imag();    
  }
  
  void Matrix::Element::operator-= (const Element& v)
  {
    std::complex<double> cplx = (const std::complex<double>)v;
    *this -= cplx;
  }


  void Matrix::Element::operator*= (const double scalar)
  {
    if( m_mtx.isReal )
    {
      m_mtx.data[m_col][m_row] *= scalar;
    }
    else
    {
      m_mtx.cplx[m_col][m_row].re *= scalar;
      m_mtx.cplx[m_col][m_row].im *= scalar;
    }
  }

  void Matrix::Element::operator*= (const std::complex<double>& v)
  {
    if( m_mtx.isReal )
    {
      if( v.imag() != 0.0 )
      {
        // This is on the fly conversion to a complex matrix!
        if( !MTX_ConvertRealToComplex( &m_mtx ) )
        {
          MTX_Free( &m_mtx );
          Matrix::StaticMatrixError( "Element::operator=", "Unable to convert matrix from real to complex" );
          return;
        }
      }
      else
      {
        m_mtx.data[m_col][m_row] *= v.real();
        return;
      }
    }

    double re = m_mtx.cplx[m_col][m_row].re;
    double im = m_mtx.cplx[m_col][m_row].im;
    
    m_mtx.cplx[m_col][m_row].re = re*v.real() - im*v.imag();
    m_mtx.cplx[m_col][m_row].im = re*v.imag() + im*v.real();    
  }
  
  void Matrix::Element::operator*= (const Element& v)
  {
    std::complex<double> cplx = (const std::complex<double>)v;
    *this *= cplx;
  }


  void Matrix::Element::operator/= (const double scalar)
  {
    //if( scalar == 0.0 )
    //{
    //  MTX_Free( &m_mtx );
    //  Matrix::StaticMatrixError("Element/=","Divide by zero!");
    //}
    //else
    //{
      if( m_mtx.isReal )
      {
        m_mtx.data[m_col][m_row] /= scalar;
      }
      else
      {
        m_mtx.cplx[m_col][m_row].re /= scalar;
        m_mtx.cplx[m_col][m_row].im /= scalar;
      }
    //}
  }

  void Matrix::Element::operator/= (const std::complex<double>& v)
  {
    if( m_mtx.isReal )
    {
      if( v.imag() == 0.0 )
      {
      //  if( v.real() == 0.0 )
      //  {
      //    MTX_Free( &m_mtx );
      //    Matrix::StaticMatrixError( "Element/=", "Divide by zero." );
      //    return;
      //  }
      //  else
      //  {
          m_mtx.data[m_col][m_row] /= v.real();
          return;
      //  }
      }
      else
      {
        // This is on the fly conversion to a complex matrix!
        if( !MTX_ConvertRealToComplex( &m_mtx ) )
        {
          MTX_Free( &m_mtx );
          Matrix::StaticMatrixError( "Element/=", "Unable to convert matrix from real to complex" );
          return;
        }
      }
    }

    double d = v.real()*v.real() + v.imag()*v.imag(); 
    //if( d == 0.0 )
    //{ 
    //  MTX_Free( &m_mtx );
    //  Matrix::StaticMatrixError( "Element/=", "Divide by zero." );
    //}
    //else
    {
      double r = m_mtx.cplx[m_col][m_row].re;
      double i = m_mtx.cplx[m_col][m_row].im;
      m_mtx.cplx[m_col][m_row].re = (r * v.real() + i * v.imag()) / d;
      m_mtx.cplx[m_col][m_row].im = (i * v.real() - r * v.imag()) / d;
    }
  }
  
  void Matrix::Element::operator/= (const Element& v)
  {
    std::complex<double> cplx = (const std::complex<double>)v;
    *this /= cplx;
  }






  //friend 
  const std::complex<double> operator+ (const Matrix::Element& m, double scalar)
  {
    std::complex<double> v;
    v = (const std::complex<double>)m;
    v += scalar;
    return v;
  }

  //friend 
  const std::complex<double> operator+ (const Matrix::Element& a, const Matrix::Element& b)
  {
    std::complex<double> v1;
    std::complex<double> v2;
    v1 = (const std::complex<double>)a;
    v2 = (const std::complex<double>)b;
    v1 += v2;
    return v1;
  }
    
  //friend 
  const std::complex<double> operator+ (const Matrix::Element& a, const std::complex<double>& b)
  {
    std::complex<double> v;
    v = (const std::complex<double>)a;
    v += b;
    return v;
  }

  //friend 
  const std::complex<double> operator+ (double scalar, const Matrix::Element& m)
  {
    return (m+scalar);

  }

  //friend 
  const std::complex<double> operator+ (const std::complex<double>& b, const Matrix::Element& a)
  {
    return (a+b);
  }


  //friend 
  const std::complex<double> operator- (const Matrix::Element& m, double scalar)
  {
    std::complex<double> v;
    v = (const std::complex<double>)m;
    v -= scalar;
    return v;
  }

  //friend 
  const std::complex<double> operator- (const Matrix::Element& a, const Matrix::Element& b)
  {
    std::complex<double> v1;
    std::complex<double> v2;
    v1 = (const std::complex<double>)a;
    v2 = (const std::complex<double>)b;
    v1 -= v2;
    return v1;
  }
    
  //friend 
  const std::complex<double> operator- (const Matrix::Element& a, const std::complex<double>& b)
  {
    std::complex<double> v;
    v = (const std::complex<double>)a;
    v -= b;
    return v;
  }

  //friend 
  const std::complex<double> operator- (double scalar, const Matrix::Element& m)
  {
    return (m+(-1.0*scalar));

  }

  //friend 
  const std::complex<double> operator- (const std::complex<double>& b, const Matrix::Element& a)
  {
    std::complex<double> v = b;
    v -= (const std::complex<double>)a;
    return v;    
  }


  //friend