Abstract non Polymorphyc Matrix:
 All Classes Namespaces Files Functions Variables Typedefs Friends Defines
Gaussian_Elimination.h
Go to the documentation of this file.
00001 // Gaussian_Elimination.h (c) 2006 adolfo@di-mare.com
00002 
00003 #ifdef Spanish_dox
00004 /// \file   Gaussian_Elimination.h
00005 /// \brief  Implementación de la eliminación Gaussiana.
00006 /// \author Adolfo Di Mare <adolfo@di-mare.com>
00007 /// \date   2009
00008 #endif
00009 #ifdef English_dox
00010 /// \file   Gaussian_Elimination.h
00011 /// \brief  Gaussian elimination implementation.
00012 /// \author Adolfo Di Mare <adolfo@di-mare.com>
00013 /// \date   2009
00014 #endif
00015 
00016 #ifndef  Gaussian_Elimination_h
00017 #define  Gaussian_Elimination_h
00018 
00019 #ifdef Spanish_dox
00020     /// Documentación en español.
00021     #define Spanish_dox "Documentación en español"
00022    /// \def Spanish_dox   ///< \c "Doxygen: Documentación en español"
00023    /// \def Gaussian_Elimination_h ///< Evita la inclusión múltiple.
00024 #endif
00025 #ifdef English_dox
00026     /// Doxygen English documentation.
00027     #define English_dox   "Doxygen English documentation"
00028    /// \def English_dox   ///< "Doxygen: English documentation"
00029    /// \def Gaussian_Elimination_h ///< Avoids multiple inclusion.
00030 #endif
00031 
00032 // #include "Matrix_BASE.h" // es extraño, pero no es necesario
00033 
00034 template <class E>
00035 inline E abs(const E& r) {
00036     return E( r>=E(0) ? r : -r );
00037 }
00038 
00039 // Resuelve M*X = B
00040 template <class Matrix>
00041 typename Matrix::value_type Gaussian_Elimination(
00042     const Matrix & M,
00043           Matrix & X,
00044     const Matrix & B
00045 ) {
00046     int i, j, k, maxrow;
00047     typename Matrix::value_type tmp;
00048 
00049     {   // verifica que "M" y "B" tengan las dimensiones correctas
00050         bool M_isSquare    = (M.rows() == M.cols());
00051         bool M_B_same_dim  = (M.cols() == B.rows());
00052         bool M_notCero     = (M.rows() != 0);
00053         bool M_B_dim       = (B.cols() >= 1);
00054         if ( ! ( M_isSquare && M_B_same_dim && M_notCero && M_B_dim) ) {
00055             // El sistema no tiene solución
00056             return typename Matrix::value_type(0);
00057         }
00058     }
00059 
00060     int n = M.rows();
00061     Matrix A;
00062     // La última columna tiene una copia de "B"
00063     A.reSize( M.rows(), M.cols()+1 );
00064 
00065     // Copia "M" en "A" y pone "B" como última columna de "A"
00066     for (i = 0; i<n; i++) {
00067         for (j = 0; j<n; j++) {
00068             A(i,j) = M(i,j);
00069         }
00070         int r = B.cols();
00071         for (j = 0; j<r; j++) {
00072             A(i,n+j) = B(i,j);
00073         }
00074     }
00075 
00076     // Algoritmo de eliminación gaussiana
00077     for (j = 0; j<n; j++) {
00078         // "maxrow" es la fila con el mayor valor de primero
00079         maxrow = j;
00080         for (i = j+1; i<n; i++) {
00081             if ( abs(A(i,j)) > abs(A(maxrow,j)) ) {
00082                 maxrow = i;
00083             }
00084         }
00085 
00086         // Intercambia la fila "maxrow" con la j-ésima
00087         for (k=j; k<n+1; k++) {
00088             tmp = A(j,k);
00089             A(j,k) = A(maxrow,k);
00090             A(maxrow,k) = tmp;
00091         }
00092 
00093         // ¿Es la matriz singular?
00094         if ( A(j,j) == typename Matrix::value_type(0) ) {
00095             return typename Matrix::value_type(0); // el determinante es cero
00096         }
00097 
00098         // Anula el j-ésimo elemento en la columna "i"
00099         for (i=j+1; i<n; i++) {
00100             for (k=n; k>=j; k--) {
00101                 A(i,k) -= A(j,k) * A(i,j) / A(j,j);
00102             }
00103         }
00104     }
00105 
00106     // Hace la sustitución hacia atrás
00107     X.reSize( M.rows(), B.cols() );
00108     for (i=n-1; i>=0; i--) { // OJO: "unsigned i" es un error
00109         // aquí puede ocurrir (i<0)
00110         tmp = typename Matrix::value_type(0);
00111         for (k=i+1; k<n; k++) {
00112             tmp += A(i,k) * X(k,0);
00113         }
00114         X(i,0) = (A(i,n) - tmp) / A(i,i);
00115     }
00116 
00117     // Calcula y retorna el determinante
00118     tmp = 1;
00119     for (i=0; i<n; i++) {
00120         tmp *= A(i,i);
00121     }
00122     return tmp;
00123 }
00124 
00125 #ifdef Spanish_dox
00126 /** \fn template <class E> \
00127         inline E abs(const E& r);
00128     \brief Retorna el valor absoluto de \c r.
00129 */
00130 /** \fn template <class Matrix> \
00131     typename Matrix::value_type Gaussian_Elimination( \
00132         const Matrix & M, \
00133               Matrix & X, \
00134         const Matrix & B \
00135     );
00136     \brief Eliminación Gaussiana.
00137     - Resuelve un sistema de \c "n" ecuaciones lineales con \c "n"
00138       incógnitas usando la "Eliminación Gaussiana".
00139     - El sistema tiene la forma matricial <code>MxX = B</code>:
00140       \code
00141       M0,0    M1,0    M2,0    ....  Mn-1,0  ==   B0,0
00142       M0,1    M1,1    M2,1    ....  Mn-1,1  ==   B0,1
00143       M0,2    M1,2    M2,2    ....  Mn-1,2  ==   B0,2
00144       :       :       :             :            :
00145       :       :       :             :            :
00146       M0,n-1  M1,n-1  M2,n-1  ....  Mn-1,n-1 ==  B0,n-1
00147       \endcode
00148 
00149     - Si la matriz \c "M" no es singular, \c Gaussian_Elimination()
00150       retorna su determinante.
00151     - Si la matriz \c "M" sí es singular, \c Gaussian_Elimination()
00152       retorna \c 0 (cero).
00153     - El resultado calculado es devuelto como valor del vector \c "X".
00154 
00155     \par Reconocimientos
00156     Adaptado de de código Internet escrito por:
00157     - Paul Bourke, 1997
00158       - http://astronomy.swin.edu.au/pbourke/analysis/gausselim/
00159       - http://www.gilsinan.com/thesis/gaussian.c
00160     - Jim Gilsinan IV
00161       - http://www.fas.harvard.edu/~gilsinan/
00162 */
00163 #endif
00164 #ifdef English_dox
00165 /** \fn template <class Matrix> \
00166     typename Matrix::value_type Gaussian_Elimination( \
00167         const Matrix & M, \
00168               Matrix & X, \
00169         const Matrix & B \
00170     );
00171     \brief Gaussianan Elimination.
00172     - Solves a system of \c "n" linear equations with \c "n"
00173       unknowns ussign "Gaussianan Elimination".
00174     - The system has the matricial form <code>MxX = B</code>:
00175       \code
00176       M0,0    M1,0    M2,0    ....  Mn-1,0  ==   B0,0
00177       M0,1    M1,1    M2,1    ....  Mn-1,1  ==   B0,1
00178       M0,2    M1,2    M2,2    ....  Mn-1,2  ==   B0,2
00179       :       :       :             :            :
00180       :       :       :             :            :
00181       M0,n-1  M1,n-1  M2,n-1  ....  Mn-1,n-1 ==  B0,n-1
00182       \endcode
00183 
00184     - If matrix \c "M" is not singular, \c Gaussian_Elimination()
00185       returns its determinant.
00186     - If matrix \c "M" is singular, \c Gaussian_Elimination()
00187       returns \c 0 (cero).
00188     - The calculated result is returned as the value of vector \c "X".
00189 
00190     \par Acknowledgments
00191     Adapted from Internet code written by:
00192     - Paul Bourke, 1997
00193       - http://astronomy.swin.edu.au/pbourke/analysis/gausselim/
00194       - http://www.gilsinan.com/thesis/gaussian.c
00195     - Jim Gilsinan IV
00196       - http://www.fas.harvard.edu/~gilsinan/
00197 */
00198 /** \fn    template <class E> inline E abs(const E& r);
00199     \brief Returns the absolute value for \c r.
00200 */
00201 #endif
00202 
00203 #endif // Gaussian_Elimination_h
00204 
00205 // EOF: Gaussian_Elimination.h