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