/*
               Solution of  Linear Equations by LU decomposition

                                 version KJI

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

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

#define SIZEM 512
#define NO_RECORD
typedef float matriz[SIZEM][SIZEM];

transM(A, AT, n)
  matriz A, AT;
  int n;
{
  int i,j;
   
  for (j=0; j<n; j++)
      {
        for (i=0; i<n; i++)
            {
              AT[j][i] = A[i][j];
            }
      }
}

printM(A, n, text)
  matriz A;
  int n;
  char * text;
{
  int x,y;
  
  printf("%s\n",text);
  for (y=0; y<n; y++)
      {
        for (x=0; x<n; x++)
            {
              printf("%f ",A[y][x]);
            }
        printf("\n");
      }
}


outtext( nome, b1, e1, b2, e2, b3, e3, total)
    char *nome;
    time_t b1, b2, b3, e1, e2, e3, total;
{
  printf("-- %10s.c :\n",nome);
  printf("-- A,L --> AT,LT:%4ld-%4ld=%3d (%02dm%02ds)\n",e1,
                                                         b1,
                                                         e1-b1,
                                                         (e1-b1)/60,
                                                         (e1-b1)%60);
  printf("-- %12s :%4ld-%4ld=%3d (%02dm%02ds)\n","decompos. LU",e2,
                                                         b2,
                                                         e2-b2,
                                                         (e2-b2)/60,
                                                         (e2-b2)%60);
  printf("-- AT,LT --> A,L:%4ld-%4ld=%3d (%02dm%02ds)\n",e3,
                                                         b3,
                                                         e3-b3,
                                                         (e3-b3)/60,
                                                         (e3-b3)%60);
                                                         
  printf("-- %12s :%4ld-%4ld=%3d (%02dm%02ds)\n","total",e3,
                                                         b1,
                                                         total,
                                                         (total)/60,
                                                         (total)%60);
}

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

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

 matriz A, AT, L, LT;
      
 time_t b1,b2,b3,e1,e2,e3,total;

/*---------------------------------------------------------------------------*/
/*                         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                       */
/*---------------------------------------------------------------------------*/

 total = 0;
 b1 = time (NULL);
   transM(A, AT, n);
   transM(L, LT, n);
 e1 = time (NULL);
 total = e1 - b1;

 b2 = time (NULL);

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

 e2 = time (NULL);
 total += e2 - b2;
 
 b3 = time (NULL);
   transM(AT, A, n);
   transM(LT, L, n);
 e3 = time (NULL);
 total += e3 - b3;
 outtext(argv[0],b1,e1,b2,e2,b3,e3,total);

/*---------------------------------------------------------------------------*/
/*                     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

}