Gaussian_Elimination.h

Ir a la documentación de este archivo.

// Gaussian_Elimination.h (c) 2006 adolfo@di-mare.com

#ifdef Spanish_dox
/// \file   Gaussian_Elimination.h
/// \brief  Implementación de la eliminación Gaussiana.
/// \author Adolfo Di Mare <adolfo@di-mare.com>
/// \date   2009
#endif
#ifdef English_dox
/// \file   Gaussian_Elimination.h
/// \brief  Gaussian elimination implementation.
/// \author Adolfo Di Mare <adolfo@di-mare.com>
/// \date   2009
#endif

#ifndef  Gaussian_Elimination_h
#define  Gaussian_Elimination_h

#ifdef Spanish_dox
    /// Documentación en español.
    #define Spanish_dox "Documentación en español"
   /// \def Spanish_dox   ///< \c "Doxygen: Documentación en español"
   /// \def Gaussian_Elimination_h ///< Evita la inclusión múltiple.
#endif
#ifdef English_dox
    /// Doxygen English documentation.
    #define English_dox   "Doxygen English documentation"
   /// \def English_dox   ///< "Doxygen: English documentation"
   /// \def Gaussian_Elimination_h ///< Avoids multiple inclusion.
#endif

// #include "Matrix_BASE.h" // es extraño, pero no es necesario

template <class E>
inline E abs(const E& r) {
    return E( r>=E(0) ? r : -r );
}

// Resuelve M*X = B
template <class Matrix>
typename Matrix::value_type Gaussian_Elimination(
    const Matrix & M,
          Matrix & X,
    const Matrix & B
) {
    int i, j, k, maxrow;
    typename Matrix::value_type tmp;

    {   // verifica que "M" y "B" tengan las dimensiones correctas
        bool M_isSquare    = (M.rows() == M.cols());
        bool M_B_same_dim  = (M.cols() == B.rows());
        bool M_notCero     = (M.rows() != 0);
        bool M_B_dim       = (B.cols() >= 1);
        if ( ! ( M_isSquare && M_B_same_dim && M_notCero && M_B_dim) ) {
            // El sistema no tiene solución
            return typename Matrix::value_type(0);
        }
    }

    int n = M.rows();
    Matrix A;
    // La última columna tiene una copia de "B"
    A.reSize( M.rows(), M.cols()+1 );

    // Copia "M" en "A" y pone "B" como última columna de "A"
    for (i = 0; i<n; i++) {
        for (j = 0; j<n; j++) {
            A(i,j) = M(i,j);
        }
        int r = B.cols();
        for (j = 0; j<r; j++) {
            A(i,n+j) = B(i,j);
        }
    }

    // Algoritmo de eliminación gaussiana
    for (j = 0; j<n; j++) {
        // "maxrow" es la fila con el mayor valor de primero
        maxrow = j;
        for (i = j+1; i<n; i++) {
            if ( abs(A(i,j)) > abs(A(maxrow,j)) ) {
                maxrow = i;
            }
        }

        // Intercambia la fila "maxrow" con la j-ésima
        for (k=j; k<n+1; k++) {
            tmp = A(j,k);
            A(j,k) = A(maxrow,k);
            A(maxrow,k) = tmp;
        }

        // ¿Es la matriz singular?
        if ( A(j,j) == typename Matrix::value_type(0) ) {
            return typename Matrix::value_type(0); // el determinante es cero
        }

        // Anula el j-ésimo elemento en la columna "i"
        for (i=j+1; i<n; i++) {
            for (k=n; k>=j; k--) {
                A(i,k) -= A(j,k) * A(i,j) / A(j,j);
            }
        }
    }

    // Hace la sustitución hacia atrás
    X.reSize( M.rows(), B.cols() );
    for (i=n-1; i>=0; i--) { // OJO: "unsigned i" es un error
        // aquí puede ocurrir (i<0)
        tmp = typename Matrix::value_type(0);
        for (k=i+1; k<n; k++) {
            tmp += A(i,k) * X(k,0);
        }
        X(i,0) = (A(i,n) - tmp) / A(i,i);
    }

    // Calcula y retorna el determinante
    tmp = 1;
    for (i=0; i<n; i++) {
        tmp *= A(i,i);
    }
    return tmp;
}

#ifdef Spanish_dox
/** \fn template <class E> \
        inline E abs(const E& r);
    \brief Retorna el valor absoluto de \c r.
*/
/** \fn template <class Matrix> \
    typename Matrix::value_type Gaussian_Elimination( \
        const Matrix & M, \
              Matrix & X, \
        const Matrix & B \
    );
    \brief Eliminación Gaussiana.
    - Resuelve un sistema de \c "n" ecuaciones lineales con \c "n"
      incógnitas usando la "Eliminación Gaussiana".
    - El sistema tiene la forma matricial <code>MxX = B</code>:
      \code
      M0,0    M1,0    M2,0    ....  Mn-1,0  ==   B0,0
      M0,1    M1,1    M2,1    ....  Mn-1,1  ==   B0,1
      M0,2    M1,2    M2,2    ....  Mn-1,2  ==   B0,2
      :       :       :             :            :
      :       :       :             :            :
      M0,n-1  M1,n-1  M2,n-1  ....  Mn-1,n-1 ==  B0,n-1
      \endcode

    - Si la matriz \c "M" no es singular, \c Gaussian_Elimination()
      retorna su determinante.
    - Si la matriz \c "M" sí es singular, \c Gaussian_Elimination()
      retorna \c 0 (cero).
    - El resultado calculado es devuelto como valor del vector \c "X".

    \par Reconocimientos
    Adaptado de de código Internet escrito por:
    - Paul Bourke, 1997
      - http://astronomy.swin.edu.au/pbourke/analysis/gausselim/
      - http://www.gilsinan.com/thesis/gaussian.c
    - Jim Gilsinan IV
      - http://www.fas.harvard.edu/~gilsinan/
*/
#endif
#ifdef English_dox
/** \fn template <class Matrix> \
    typename Matrix::value_type Gaussian_Elimination( \
        const Matrix & M, \
              Matrix & X, \
        const Matrix & B \
    );
    \brief Gaussianan Elimination.
    - Solves a system of \c "n" linear equations with \c "n"
      unknowns ussign "Gaussianan Elimination".
    - The system has the matricial form <code>MxX = B</code>:
      \code
      M0,0    M1,0    M2,0    ....  Mn-1,0  ==   B0,0
      M0,1    M1,1    M2,1    ....  Mn-1,1  ==   B0,1
      M0,2    M1,2    M2,2    ....  Mn-1,2  ==   B0,2
      :       :       :             :            :
      :       :       :             :            :
      M0,n-1  M1,n-1  M2,n-1  ....  Mn-1,n-1 ==  B0,n-1
      \endcode

    - If matrix \c "M" is not singular, \c Gaussian_Elimination()
      returns its determinant.
    - If matrix \c "M" is singular, \c Gaussian_Elimination()
      returns \c 0 (cero).
    - The calculated result is returned as the value of vector \c "X".

    \par Acknowledgments
    Adapted from Internet code written by:
    - Paul Bourke, 1997
      - http://astronomy.swin.edu.au/pbourke/analysis/gausselim/
      - http://www.gilsinan.com/thesis/gaussian.c
    - Jim Gilsinan IV
      - http://www.fas.harvard.edu/~gilsinan/
*/
/** \fn    template <class E> inline E abs(const E& r);
    \brief Returns the absolute value for \c r.
*/
#endif

#endif // Gaussian_Elimination_h

// EOF: Gaussian_Elimination.h