matrix.h

复数矩阵的运算

        _m = m._m;
        return *this;
    }

    //  reallocation method
    MAT_TEMPLATE inline void 
        matrixT::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(T);
        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
    MAT_TEMPLATE inline void
        matrixT::SetSize (size_t row, size_t col) _NO_THROW
    {
        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] = T(0);

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

        return;
    }

    // subscript operator to get/set individual elements
    MAT_TEMPLATE inline T&
        matrixT::operator [][] (size_t row, size_t col) _THROW_MATRIX_ERROR
    {
        if (row >= _m->Row || col >= _m->Col)
            REPORT_ERROR( "matrixT::operator(): Index out of range!");
        if (_m->Refcnt > 1) clone();
        return _m->Val[row][col];
    }

    // subscript operator to get/set individual elements
    MAT_TEMPLATE inline T
        matrixT::operator [][] (size_t row, size_t col) const _THROW_MATRIX_ERROR
    {
        if (row >= _m->Row || col >= _m->Col)
            REPORT_ERROR( "matrixT::operator(): Index out of range!");
        return _m->Val[row][col];
    }

    // input stream function
    MAT_TEMPLATE inline istream&
        operator >> (istream& istrm, matrixT& m)
    {
        for (size_t i=0; i < m.RowNo(); i++)
            for (size_t j=0; j < m.ColNo(); j++)
            {
                T x;
                istrm >> x;
                m(i,j) = x;
            }
            return istrm;
    }

    // output stream function
    MAT_TEMPLATE inline ostream&
        operator << (ostream& ostrm, const matrixT& m)
    {
        for (size_t i=0; i < m.RowNo(); i++)
        {
            for (size_t j=0; j < m.ColNo(); j++)
            {
                T x = m(i,j);
                ostrm << x << '\t';
            }
            ostrm << endl;
        }
        return ostrm;
    }


    //// logical equal-to operator
    //MAT_TEMPLATE inline bool
    //    operator == (const matrixT& m1, const matrixT& m2) _NO_THROW
    //{
    //    if (m1.RowNo() != m2.RowNo() || m1.ColNo() != m2.ColNo())
    //        return false;

    //    for (size_t i=0; i < m1.RowNo(); i++)
    //        for (size_t j=0; j < m1.ColNo(); j++)
    //            if (m1(i,j) != m2(i,j))
    //                return false;

    //    return true;
    //}

    //// logical no-equal-to operator
    //MAT_TEMPLATE inline bool
    //    operator != (const matrixT& m1, const matrixT& m2) _NO_THROW
    //{
    //    return (m1 == m2) ? false : true;
    //}

    // combined addition and assignment operator
    MAT_TEMPLATE inline matrixT&
        matrixT::operator += (const matrixT& m) _THROW_MATRIX_ERROR
    {
        if (_m->Row != m._m->Row || _m->Col != m._m->Col)
            REPORT_ERROR( "matrixT::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
    MAT_TEMPLATE inline matrixT&
        matrixT::operator -= (const matrixT& m) _THROW_MATRIX_ERROR
    {
        if (_m->Row != m._m->Row || _m->Col != m._m->Col)
            REPORT_ERROR( "matrixT::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
    MAT_TEMPLATE inline matrixT&
        matrixT::operator *= (const T& c) _NO_THROW
    {
        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
    MAT_TEMPLATE inline matrixT&
        matrixT::operator *= (const matrixT& m) _THROW_MATRIX_ERROR
    {
        if (_m->Col != m._m->Row)
            REPORT_ERROR( "matrixT::operator*= : Inconsistent matrix sizes in multiplication!");

        matrixT 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] = T(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
    MAT_TEMPLATE inline matrixT&
        matrixT::operator /= (const T& c) _NO_THROW
    {
        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
    MAT_TEMPLATE inline matrixT&
        matrixT::operator ^= (const size_t& pow) _THROW_MATRIX_ERROR
    {
        matrixT temp(*this);

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

        return *this;
    }

    // unary negation operator
    // 单目运算：取负
    MAT_TEMPLATE inline matrixT
        matrixT::operator - () _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._m->Val[i][j] = - _m->Val[i][j];

        return temp;
    }

    // binary addition operator
    MAT_TEMPLATE inline matrixT
        operator + (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
    {
        matrixT temp = m1;
        temp += m2;
        return temp;
    }

    // binary subtraction operator
    MAT_TEMPLATE inline matrixT
        operator - (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
    {
        matrixT temp = m1;
        temp -= m2;
        return temp;
    }

    // binary scalar multiplication operator
    MAT_TEMPLATE inline matrixT
        operator * (const matrixT& m, const T& no) _NO_THROW
    {
        matrixT temp = m;
        temp *= no;
        return temp;
    }


    // binary scalar multiplication operator
    MAT_TEMPLATE inline matrixT
        operator * (const T& no, const matrixT& m) _NO_THROW
    {
        return (m * no);
    }

    // binary matrix multiplication operator
    MAT_TEMPLATE inline matrixT
        operator * (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
    {
        matrixT temp = m1;
        temp *= m2;
        return temp;
    }

    // binary scalar division operator
    MAT_TEMPLATE inline matrixT
        operator / (const matrixT& m, const T& no) _NO_THROW
    {
        return (m * (T(1) / no));
    }


    // binary scalar division operator
    MAT_TEMPLATE inline matrixT
        operator / (const T& no, const matrixT& m) _THROW_MATRIX_ERROR
    {
        return (!m * no);
    }

    // binary matrix division operator
    MAT_TEMPLATE inline matrixT
        operator / (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
    {
        return (m1 * !m2);
    }

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

    // unary transpose operator
    // 矩阵转置
    MAT_TEMPLATE inline matrixT
        operator ~ (const matrixT& m) _NO_THROW
    {
        matrixT temp(m.ColNo(),m.RowNo());

        for (size_t i=0; i < m.RowNo(); i++)
            for (size_t j=0; j < m.ColNo(); j++)
            {
                T x = m(i,j);
                temp(j,i) = x;
            }
            return temp;
    }

    // unary inversion operator
    // 矩阵求逆
    MAT_TEMPLATE inline matrixT
        operator ! (const matrixT m) _THROW_MATRIX_ERROR
    {
        matrixT temp = m;
        return temp.Inv();
    }

    // inversion function
    // 矩阵求逆
    MAT_TEMPLATE inline matrixT
        matrixT::Inv () _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);
        if (_m->Refcnt > 1) clone();


        temp.Unit();
        for (k=0; k < _m->Row; k++)
        {
            int indx = pivot(k);
            if (indx == -1)
                REPORT_ERROR( "matrixT::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->Row; 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->Row; j++)