todo_inline.cpp

   1 /* ec_lineal.h
   2  ***********************
   3  * Análisis Numérico I *
   4  ***********************
   5  */
   6
   7 // Tipos básicos
   8 typedef unsigned int Indice;
   9 typedef float Numero;
  10
  11 // Cantidad máxima de elementos
  12 const Indice N = 1000;
  13
  14 // Tipos compuestos
  15 typedef Indice VectorPermutaciones[N];
  16 typedef Numero Vector[N];
  17 typedef Numero Matriz[N][N];
  18
  19 #include <stdio.h>
  20 #include <stdlib.h>
  21 #include <math.h>
  22 #include <time.h>
  23
  24 inline void generar_vector( Vector& v ) {
  25     for ( Indice i = 0; i < N; i++ )
  26         v[i] = Numero( rand() ) * 3.0 / RAND_MAX;
  27 }
  28
  29 inline void generar_vector_permutaciones( VectorPermutaciones& v ) {
  30     for ( Indice i = 0; i < N; i++ )
  31         v[i] = i;
  32 }
  33
  34 inline void imprimir_matriz( Matriz& m ) {
  35     for ( Indice i = 0; i < N; i++ ) {
  36         for ( Indice j = 0; j < N; j++ )
  37             printf( "%.02f ", m[i][j] );
  38         printf( "\n" );
  39     }
  40 }
  41
  42 inline void imprimir_vector( Vector& v ) {
  43     for ( Indice i = 0; i < N; i++ )
  44         printf( "%.02f ", v[i] );
  45     printf( "\n" );
  46 }
  47
  48 inline void imprimir_matriz_permutada( Matriz& m, VectorPermutaciones& p ) {
  49     for ( Indice i = 0; i < N; i++ ) {
  50         for ( Indice j = 0; j < N; j++ )
  51             printf( "%.02f ", m[p[i]][j] );
  52         printf( "\n" );
  53     }
  54 }
  55
  56 inline void imprimir_vector_permutado( Vector& v, VectorPermutaciones& p ) {
  57     for ( Indice i = 0; i < N; i++ )
  58         printf( "%.02f ", v[p[i]] );
  59     printf( "\n" );
  60 }
  61
  62 inline void imprimir_matriz_L( Matriz& m, VectorPermutaciones& p ) {
  63     for ( Indice i = 0; i < N; i++ ) {
  64         for ( Indice j = 0; j < N; j++ )
  65             if ( i < j )
  66                 printf( "0.00 " );
  67             else if ( i == j )
  68                 printf( "1.00 " );
  69             else
  70                 printf( "%.02f ", m[p[i]][j] );
  71         printf( "\n" );
  72     }
  73 }
  74
  75 inline void imprimir_matriz_U( Matriz& m, VectorPermutaciones& p ) {
  76     for ( Indice i = 0; i < N; i++ ) {
  77         for ( Indice j = 0; j < N; j++ )
  78             if ( i > j )
  79                 printf( "0.00 " );
  80             else
  81                 printf( "%.02f ", m[p[i]][j] );
  82         printf( "\n" );
  83     }
  84 }
  85
  86 inline void sustitucion_directa( Matriz& A, Vector& b, VectorPermutaciones& p ) {
  87     for ( Indice i = 1; i < N; i++ ) {
  88         Numero sum = 0;
  89         for ( Indice j = 0; j < i - 1; j++ )
  90             sum += A[p[i]][j] * b[p[i]];
  91         b[p[i]] -= sum;
  92     }
  93 }
  94
  95 inline void sustitucion_inversa( Matriz& A, Vector& b, VectorPermutaciones& p ) {
  96     for ( Indice i = N; i > 0; i-- ) {
  97         Numero sum = 0;
  98         for ( Indice j = i; j < N; j++ )
  99             sum += A[p[i-1]][j] * b[p[j]];
 100         b[p[i-1]] = ( b[p[i-1]] - sum ) / A[p[i-1]][i];
 101     }
 102 }
 103 // Constantes de error
 104 const int MATRIZ_SINGULAR = 1;
 105
 106 int main( int argc, char *argv[] ) {
 107
 108     Matriz A;
 109     Vector b;
 110     VectorPermutaciones p;
 111
 112     // Reinicio los numeros aleatorios
 113     //srand( time( NULL ) );
 114
 115     // Genero matriz al azar
 116     generar_matriz( A );
 117     // Genero vector al azar
 118     generar_vector( b );
 119     // Genero vector de permutaciones
 120     generar_vector_permutaciones( p );
 121
 122     // La imprimo
 123     printf( "Matriz A:\n=========\n" );
 124     imprimir_matriz( A );
 125     printf( "\n" );
 126     // La imprimo
 127     printf( "Vector b:\n=========\n" );
 128     imprimir_vector( b );
 129     printf( "\n\n" );
 130
 131     // Recorro por columna (o diagonal).
 132     for ( Indice j = 0; j < N; j++ ) {
 133         // Busco máximo en la columna para usar de pivote
 134         Indice maxi = j;
 135         Numero max  = fabs( A[j][j] );
 136         for ( Indice i = j + 1; i < N; i++ )
 137             if ( fabs( A[i][j] ) > max ) {
 138                 maxi = i;
 139                 max  = fabs( A[i][j] );
 140             }
 141
 142         // Nos fijamos si la matriz es singular
 143         if ( max == 0 )
 144             return MATRIZ_SINGULAR;
 145
 146         // Nos fijamos si hay permutación
 147         if ( maxi != j ) {
 148             // Intercambio vector de permutaciones
 149             Indice aux  = p[maxi];
 150             p[maxi]     = p[j];
 151             p[j]        = aux;
 152         }
 153
 154         // Reduzco las filas
 155         for ( Indice i = j + 1; i < N; i++ ) {
 156             // Calculo el m, almacenandolo en A
 157             A[p[i]][j] /= A[p[j]][j];
 158             // Recorro columna por columna, operando
 159             for ( Indice k = j + 1; i < N; i++ )
 160                 A[p[i]][j] -= A[p[i]][j] * A[p[j]][k]; // Aij = Aij - m * Ajk
 161         }
 162     }
 163
 164     // Imprimo la matriz modificada
 165     printf( "Matriz A modificada:\n====================\n" );
 166     imprimir_matriz_permutada( A, p );
 167     printf( "\n" );
 168     // Imprimo la matriz L
 169     printf( "Matriz L:\n========\n" );
 170     imprimir_matriz_L( A, p );
 171     printf( "\n" );
 172     // Imprimo la matriz U
 173     printf( "Matriz U:\n========\n" );
 174     imprimir_matriz_U( A, p );
 175     printf( "\n" );
 176
 177     // Calcula solucion de Ly=b
 178     sustitucion_directa( A, b, p );
 179
 180     // Imprimo el vector solucion y
 181     printf( "y (Ly=b):\n=========\n" );
 182     imprimir_vector_permutado( b, p );
 183     printf( "\n" );
 184
 185     // Calcula solucion de Ux=y
 186     sustitucion_inversa( A, b, p );
 187
 188     // Imprimo el vector solucion x
 189     printf( "x (Ux=y):\n=========\n" );
 190     imprimir_vector_permutado( b, p );
 191     printf( "\n" );
 192
 193     return EXIT_SUCCESS;
 194
 195 }
 196
 197 inline void generar_matriz( Matriz& m ) {
 198     for ( Indice i = 0; i < N; i++ )
 199         for ( Indice j = 0; j < N; j++ )
 200             m[i][j] = Numero( rand() ) * 3.0 / RAND_MAX;
 201 }
 202
 203