matrix.h

复数矩阵的运算

                    {
                        _m->Val[i][j] -= a2 * _m->Val[k][j];
                        temp._m->Val[i][j] -= a2 * temp._m->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 (!(_m->Row == _m->Col && _m->Col == v._m->Row))
            REPORT_ERROR( "matrixT::Solve():Inconsistent matrices!");

        matrixT 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( "matrixT::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];
            }
        }
        matrixT 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
    MAT_TEMPLATE inline void
        matrixT::Null (const size_t& row, const size_t& col) _NO_THROW
    {
        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] = T(0);
        return;
    }

    // set zero to all elements of this matrix
    MAT_TEMPLATE inline void
        matrixT::Null() _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] = T(0);
        return;
    }

    // set this matrix to unity
    MAT_TEMPLATE inline void
        matrixT::Unit (const size_t& row) _NO_THROW
    {
        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 ? T(1) : T(0);
        return;
    }

    // set this matrix to unity
    MAT_TEMPLATE inline void
        matrixT::Unit () _NO_THROW
    {
        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 ? 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 < _m->Row; i++)
            if ( (temp = abs( _m->Val[i][row])) > amax && temp != 0.0)
            {
                amax = temp;
                k = i;
            }
            if (_m->Val[k][row] == T(0))
                return -1;
            if (k != int(row))
            {
                T* rowptr = _m->Val[k];
                _m->Val[k] = _m->Val[row];
                _m->Val[row] = rowptr;
                return k;
            }
            return 0;
    }

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

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

        matrixT 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
    MAT_TEMPLATE inline T
        matrixT::Norm () _NO_THROW
    {
        T retVal = T(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
    MAT_TEMPLATE inline T
        matrixT::Cond () _NO_THROW
    {
        matrixT inv = ! (*this);
        return (Norm() * inv.Norm());
    }

    // 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 (_m->Row != _m->Col)
            REPORT_ERROR( "matrixT::Cofact(): Cofactor of a non-square matrix!");

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

        matrixT 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++;
        }
        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 (_m->Row != _m->Col)
            REPORT_ERROR( "matrixT::Adj(): Adjoin of a non-square matrix.");

        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[j][i] = Cofact(i,j);
        return temp;
    }

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

    // Determine if the matrix is diagonal
    MAT_TEMPLATE inline bool
        matrixT::IsDiagonal () _NO_THROW
    {
        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] != 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 = _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
    MAT_TEMPLATE inline bool
        matrixT::IsUnit () _NO_THROW
    {
        if (IsScalar() && _m->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 < _m->Row; i++)
            for (size_t j=0; j < _m->Col; j++)
                if (_m->Val[i][j] != T(0))
                    return false;
        return true;
    }

    // Determine if the matrix is symmetric
    MAT_TEMPLATE inline bool
        matrixT::IsSymmetric () _NO_THROW
    {
        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 (_m->Val[i][j] != _m->Val[j][i])
                    return false;
        return true;
    }

    // Determine if the matrix is skew-symmetric
    MAT_TEMPLATE inline bool
        matrixT::IsSkewSymmetric () _NO_THROW
    {
        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 (_m->Val[i][j] != -_m->Val[j][i])
                    return false;
        return true;
    }

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

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

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

        return true;
    }

#ifndef _NO_NAMESPACE
} 
#endif

#endif //__STD_MATRIX_H