/*
               Solution of  Linear Equations by LU decomposition

                                 version KIJ

FILE: lu-kij.c           date: 26.09.1994.Mo
*/

#include <stdio.h>
#include <sys/types.h>
#include <sys/time.h>

#define SIZEM 512
#define NO_RECORD

outtext( nome, beginning, ending)
    char *nome;
    time_t beginning, ending;
{
  printf("%5s.c :%4ld-%4ld=%4ld (%02dm%02ds)\n" ,nome,
                                                 ending,
                                                 beginning,
                                                 ending-beginning,
                                                 (ending-beginning)/60,
                                                 (ending-beginning)%60);
}

main(argc, argv)
     int argc;
     char **argv;
{
 FILE *ArqA,
      *ArqL,
      *ArqU;

 int i, j, k, s, n,
     rows, columns;

 float A[SIZEM][SIZEM],
       L[SIZEM][SIZEM],
       U[SIZEM][SIZEM];

 time_t beginning, ending, delta;

/*---------------------------------------------------------------------------*/
/*                         Loading matrix A                                  */
/*---------------------------------------------------------------------------*/

 ArqA = fopen ("A","rb");

 fread (&rows, sizeof(int), 1, ArqA);
 fread (&columns, sizeof(int), 1, ArqA);

 for (i=0; i<rows; i++)
     for (j=0; j<columns; j++)
         {
	  fread(&A[i][j], sizeof(float), 1, ArqA);
         }

 fclose(ArqA);

 n = rows;

/*---------------------------------------------------------------------------*/
/*                    Initialization of matrices                             */
/*---------------------------------------------------------------------------*/

 for (i=0; i<n; i++)
     {
      for (j=i+1; j<n; j++)
	  L[i][j] = 0;

/*      for (j=0; j<i; j++)
	  U[i][j] = 0;*/

      L[i][i] = 1;
     }

/*---------------------------------------------------------------------------*/
/*               Beginning of LU decomposition routine                       */
/*---------------------------------------------------------------------------*/

 beginning = time (NULL);

 for (k=0; k<n-1; k++)
     {
      for (i=k+1; i<n; i++)
	  {
	    L[i][k] = A[i][k]/A[k][k];
 	    for (j=k+1; j<n; j++)
              {
 	        A[i][j] -= L[i][k] * A[k][j];
              }
          }
     }

 ending = time (NULL);

 outtext (argv[0], beginning, ending);

/*---------------------------------------------------------------------------*/
/*                     Saving results in files                               */
/*---------------------------------------------------------------------------*/

#ifndef NO_RECORD

 ArqU = fopen ("U.calc","wb");
 fwrite (&n, sizeof(int), 1, ArqU);
 fwrite (&n, sizeof(int), 1, ArqU);

 for (i=0; i<n; i++)
     for (j=0; j<n; j++)
         {
          fwrite(&A[i][j], sizeof(float), 1, ArqU);
         }

 fclose (ArqU);

 ArqL = fopen ("L.calc","wb");
 fwrite (&n, sizeof(int), 1, ArqL);
 fwrite (&n, sizeof(int), 1, ArqL);

 for (i=0; i<n; i++)
     for (j=0; j<n; j++)
         {
          fwrite(&L[i][j], sizeof(float), 1, ArqL);
         }

 fclose (ArqL);

#endif

}