lineal_iterativo.cpp

   1 /* lineal_iterativo.cpp
   2  **********************
   3  * Análisis Numérico I
   4  *---------------------
   5  * Para direccionar la salida a un archivo, ejecutelo con:
   6  *
   7  *  lineal_iterativo > [archivo]
   8  *
   9  *  Para compilar en *nix:
  10  *
  11  *  gcc -lm -Wall -O3 -o lineal_iterativo lineal_iterativo.cpp
  12  *
  13  *  Para debug:
  14  *  gcc -lm -Wall -g3 -o lineal_iterativo lineal_iterativo.cpp
  15  *
  16  */
  17
  18 #include <time.h>
  19 #include "lib.h"
  20
  21 // Constantes de error
  22 const int MATRIZ_SINGULAR    = 1;
  23 const int MAXIMO_ITERACIONES = 2;
  24
  25 // Declaración de funciones
  26 Numero* cargar_matriz( const char* );
  27
  28 // Variables
  29 Matriz A;
  30 Vector b;
  31 VectorPermutaciones p;
  32
  33 int main( int argc, char *argv[] ) {
  34
  35     // Reinicio los numeros aleatorios
  36     //srand( time( NULL ) );
  37
  38     // Genero matriz al azar
  39     generar_matriz( A );
  40     // Genero vector al azar
  41     generar_vector( b );
  42     // Genero vector de permutaciones
  43     generar_vector_permutaciones( p );
  44
  45     // Imprimo matriz
  46     printf( "Matriz A:\n=========\n" );
  47     imprimir_matriz( A );
  48     printf( "\n" );
  49     // Imprimo vector
  50     printf( "Vector b:\n=========\n" );
  51     imprimir_vector( b );
  52     printf( "\n\n" );
  53
  54     // Recorro por columna (o diagonal).
  55     for ( Indice j = 0; j < N; j++ ) {
  56         // Busco máximo en la columna para usar de pivote
  57         Indice maxi = j;
  58         Numero max  = fabs( A[j][j] );
  59         for ( Indice i = j + 1; i < N; i++ )
  60             if ( fabs( A[i][j] ) > max ) {
  61                 maxi = i;
  62                 max  = fabs( A[i][j] );
  63             }
  64
  65         // Nos fijamos si la matriz es singular
  66         if ( max == 0 )
  67             return MATRIZ_SINGULAR;
  68
  69         // Nos fijamos si hay permutación
  70         if ( maxi != j ) {
  71             // Intercambio vector de permutaciones
  72             Indice aux  = p[maxi];
  73             p[maxi]     = p[j];
  74             p[j]        = aux;
  75         }
  76
  77         // Reduzco las filas
  78         for ( Indice i = j + 1; i < N; i++ ) {
  79             // Calculo el m, almacenandolo en A
  80             A[p[i]][j] /= A[p[j]][j];
  81             // Recorro columna por columna, operando
  82             for ( Indice k = j + 1; i < N; i++ )
  83                 A[p[i]][k] -= A[p[i]][j] * A[p[j]][k]; // Aij = Aij - m * Ajk
  84         }
  85     }
  86
  87     // Imprimo la matriz modificada
  88     printf( "Matriz A modificada:\n====================\n" );
  89     imprimir_matriz_permutada( A, p );
  90     printf( "\n" );
  91     // Imprimo la matriz L
  92     printf( "Matriz L:\n========\n" );
  93     imprimir_matriz_L( A, p );
  94     printf( "\n" );
  95     // Imprimo la matriz U
  96     printf( "Matriz U:\n========\n" );
  97     imprimir_matriz_U( A, p );
  98     printf( "\n" );
  99
 100     // Calcula solucion de Ly=b
 101     sustitucion_directa( A, b, p );
 102
 103     // Imprimo el vector solucion y
 104     printf( "y (Ly=b):\n=========\n" );
 105     imprimir_vector_permutado( b, p );
 106     printf( "\n" );
 107
 108     // Calcula solucion de Ux=y
 109     sustitucion_inversa( A, b, p );
 110
 111     // Imprimo el vector solucion x
 112     printf( "x (Ux=y):\n=========\n" );
 113     imprimir_vector_permutado( b, p );
 114     printf( "\n" );
 115
 116     return EXIT_SUCCESS;
 117
 118 }
 119
 120 Numero* cargar_matriz( const char* archivo ) {
 121
 122     Indice n;
 123
 124     FILE* f;
 125     f = fopen( archivo, "r" );
 126
 127     char s[4096];
 128     n = Indice( atol( fgets( s, 4095, f ) ) );
 129     for ( Indice i = 0; !feof( f ) || ( i > n ); i++ ) {
 130         fgets( s, 4095, f );
 131