/*****************************************************************************/
/*                                                                           */
/*     Project: Homework 2 - Exercise 2.12                                   */
/*   Copyright: Copyright © 2003. All Rights Reserved                        */
/*    Compiler: Microsoft Visual C++ .NET, Linux gcc 2.95.3                  */
/*     Company: RT Mnemonic <mnemonic@eunet.at>                              */
/*      Author: Rainer Trummer                                               */
/*        Date: November 9, 2003                                             */
/*                                                                           */
/*****************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <time.h>


inline void randomize( void )
{
    srand( (unsigned int) time( NULL ) );
}

inline double random_double( void )
{
    return( (double) RAND_MAX / (double)(rand( ) + 1) );
}


int main( void )
{
    const int     N = 1000;

    double        **L, *x, *b;
    unsigned long i, j, k = 0, hN = N >> 1;

    try
    {
        L = new double *[N];

        for( i = 0; i < N; ++i )
        {
            L[i] = new double[N];
        }

        x = new double[N];
        b = new double[N];
    }
    catch( ... )
    {
        printf( "Not enough memory!\n\n" );
        return( 1 );
    }
        
    randomize( );

    // Create lower triangular matrix L and vector b
    //
    for( i = 0; i < N; ++i )
    {
        for( j = 0; j < N; ++j )
        {
            L[i][j] = ( i < j ) ? 0.0 : random_double( );
        }

        b[i] = random_double( );
    }

    // Print out the initialized linear system
    //
    if( N < 7 )
    {
        for( i = 0; i < N; ++i )
        {
            printf( "|" );

            for( j = 0; j < N; ++j )
            {
                printf( " %6.2lf", L[i][j] );
            }

            if( i == hN )
            {
                printf( " |  *  | x%d |  =  | %6.2lf |\n", i, b[i] );
            }
            else
            {
                printf( " |     | x%d |     | %6.2lf |\n", i, b[i] );
            }
        }

        printf( "\n\n" );
    }

    // Apply forward-substitution to the system
    //
    for( j = 0; j < N; ++j )
    {
        if( L[j][j] == 0.0 )
        {
            printf( "Singular Matrix!\n\n" );
            break;
        }

        x[j] = b[j] / L[j][j];
        ++k;

        for( i = j + 1; i < N; ++i )
        {
            b[i] -= L[i][j] * x[j];
            ++k;
        }
    }

    // Print out the solved linear system
    //
    if( N < 7 )
    {
        for( i = 0; i < N; ++i )
        {
            printf( "|" );

            for( j = 0; j < N; ++j )
            {
                printf( " %6.2lf", L[i][j] );
            }

            if( i == hN )
            {
                printf( " |  *  | %6.2lf |  =  | %6.2lf |\n", x[i], b[i] );
            }
            else
            {
                printf( " |     | %6.2lf |     | %6.2lf |\n", x[i], b[i] );
            }
        }

        printf( "\n\n" );
    }

    printf( "Operation count: %d\n\n"
            "( N^2 + N ) / 2: %d\n\n", k, (N * N + N) >> 1 );

    for( i = 0; i < N; ++i )
    {
        delete [] L[i];
    }

    delete [] L;
    delete [] x;
    delete [] b;
    return( 0 );
}

// End of file.