Logo Search packages:      
Sourcecode: sailcut version File versions  Download package

CMatrix CMatrix::gaussjordan ( bool *  is_inv = NULL,
CMatrix inv = NULL,
soltype_t *  soltype = NULL,
CVector bb = NULL,
CMatrix tkern = NULL 
) const

Diagonalises matrix by Gauss-Jordan method with full pivoting optionally returns the inverse matrix.

Definition at line 203 of file matrix.cpp.

References CMatrix(), col(), crop(), CVector::getdim(), id(), m_ncol, m_nrow, CVector::matrix(), swap_col(), and swap_row().

Referenced by diag(), inv(), kern(), and solve().

{
    // check dimensions
    CMatrix b;
    if (bb!=NULL)
    {
        if (bb->getdim() != m_nrow)
            throw CException("CMatrix::solve : matrix <=> right-hand side dimensions incompatible");

        b = bb->matrix();
    }
    else
    {
        b = CMatrix(m_nrow,1);
    }
    //cout << "b" << endl<< b << endl;
    //cout << "m" << endl<< *this << endl;

    real p_value;
    real p_avalue, avalue;
    unsigned int p_i, p_j;
    unsigned int i,j,k;
    unsigned int n;
    // copy over the matrix
    CMatrix m = *this;
    CMatrix t;
    // left-hand operations (row ops)
    CMatrix lops = CMatrix::id(m_nrow);
    //  CMatrix rswap = CMatrix::id(getnrow());
    // right-hand operations (col ops)
    CMatrix cswap = CMatrix::id(m_ncol);

    // determine minimum dimension
    if (m_ncol < m_nrow)
        n = m_ncol;
    else
        n = m_nrow;


    ////////////////////////////////
    //     DIAGONALISATION        //
    ////////////////////////////////

    for (k=0; k < n; k++)
    {
        // find pivot
        p_value = p_avalue = 0;
        p_i = p_j = 0;
        for (i = k; i < m_nrow; i++)
        {
            for (j = k; j < m_ncol; j++)
            {
                avalue = fabs(m(i,j));
                if (avalue > p_avalue)
                {
                    p_value = m(i,j);
                    p_avalue = avalue;
                    p_i = i;
                    p_j = j;
                }
            }
        }
        // if pivot is smaller than epsillon we cannot continue diagonalisation
        if ( p_avalue < EPS )
            break;

        // exchange rows k and p_i
        // swapping
        m.swap_row(k,p_i);
        lops.swap_row(k,p_i);
        //rswap.swap_row(k,p_i);
        b.swap_row(k,p_i);
        // swapping cols k and p_j
        m.swap_col(k,p_j);
        cswap.swap_col(k,p_j);
        // modify rows
        t = CMatrix::id(m_nrow);
        for ( i =0; i < m_nrow; i++)
        {
            if ( i !=k )
                t(i,k) = -m(i,k) / p_value;
            else
                t(k,k) =  1 / p_value;
        }
        m = t * m;
        b = t * b;
        lops = t * lops;
    }


    /////////////////////////////////////
    //  HANDLE REQUEST FOR INVERSION   //
    /////////////////////////////////////

    if ((is_inv != NULL) && (inv != NULL))
    {
        if ((m_nrow == m_ncol)&&(m_ncol == k))
        {
            // matrix is square and invertible
            *inv = cswap*lops;
            *is_inv = true;
        }
    }


    /////////////////////////////////////
    //  HANDLE REQUEST FOR SOLUTION    //
    /////////////////////////////////////

    if ((bb != NULL)&&(soltype != NULL))
    {
        *soltype = ONE;
        for (i = k; i < b.m_nrow; i++)
            if (b(i,0) > EPS)
            {
                *soltype = NONE;
                break;
            }

        if(*soltype==NONE)
        {
            // incompatible system
            *bb = CVector(0);
        }
        else
        {
            // system compatible one or infinity of solutions
            *bb = (cswap * b.crop(m_ncol,1)).col(0);
            if (m_ncol == k)
            {
                // complete diagonalisation, one solution
                *soltype=ONE;
            }
            else
            {
                *soltype=INF;
            }
        }
    }

    /////////////////////////////////////
    //  HANDLE REQUEST FOR KERNEL     //
    /////////////////////////////////////

    if (tkern != NULL)
    {
        if (k < m_ncol)
        {
            /*
                  cout << "has non-zero kernel" << endl;
                  if ((getnrow()==k)) {
                    // no redundant equations
                    cout << "getnrow()==k==" << k << endl;
                  } else {
                    cout << "getnrow()==" << getnrow() << " k==" << k << endl;
                  }
            */
            CMatrix pk = - m.crop(m_ncol-k+1,m_ncol-k,0,k);
            pk = pk.crop(m_ncol,m_ncol-k);
            for (i=0; i + k < m_ncol; i++)
                pk(i+k,i) = 1;
            *tkern = cswap* pk;
        }
        else
        {
            // matrix is inversible
            //cout << "m_ncol==k" << endl;
            *tkern = CMatrix();
        }
    }
    // clean matrix

    for ( i = 0; i < m_nrow; i++ )
        for ( j = 0; j < m_ncol; j++)
            if (fabs(m(i,j)) < EPS)
                m(i,j) = 0;

    //  return m;

    return m.crop(k,m_ncol);
}


Generated by  Doxygen 1.6.0   Back to index