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

/*
  esercitazione3.c
  Copyright (C) 2008 by Giovanni.Organtini@roma1.infn.it

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

#define NMAX 100

main() {
  /* algoritmo di crout per la decomposizione di una matrice */

  float min, max, a, b;
  float m[NMAX][NMAX];
  float alpha[NMAX][NMAX] = {0.};
  float beta[NMAX][NMAX] = {0.};
  float res[NMAX][NMAX];
  int i, j, k, N;

  srand(time(NULL));

  /* presentazione */
  printf("Questo programma decompone una matrice quadrata in due matrici\n"
	 "triangolari usando l'algoritmo di Crout.\n\n");
  do {
    /* notate l'uso della costante nel messaggio */
    printf("Inserisci l'ordine della matrice [2,%d): ", NMAX);
    scanf("%d", &N);
  } while ((N < 2) || (N > NMAX));
  printf("Inserisci i limiti entro cui possono variare gli elementi\n"
	 "della matrice, separati da uno spazio: ");
  do {
    scanf("%f %f", &a, &b);
    if (a >= b) {
      printf("ERRORE: il primo valore deve essere minore del secondo: ");
    }
  } while (a >= b);

  /* genera la matrice random */
  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      m[i][j] = (b - a)*rand()/RAND_MAX + a;
    }
  }

  /* stampa */
  printf("M\n=============================\n");
  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      printf("%+f ", m[i][j]);
    }
    printf("\n");
  }
  printf("\n");

  /* step 1: assegna i valori della diagonale */
  for (i = 0; i < N; i++) {
    alpha[i][i] = 1.;
  }

  /* step 2: calcola gli elementi fuori della diagonale */
  for (j = 0; j < N; j++) {
    for (i = 0; i <= j; i++) {
      /* notate come si ottimizza il codice scrivendo nell'ordine
	 giusto le istruzioni */
      beta[i][j] = m[i][j];
      for (k = 0; k <= i - 1; k++) {
	beta[i][j] -= alpha[i][k]*beta[k][j];
      }
    }
    for (i = j + 1; i < N; i++) {
      alpha[i][j] = m[i][j];
      for (k = 0; k <= j - 1; k++) {
	alpha[i][j] -= alpha[i][k] * beta[k][j];
      }
      alpha[i][j] /= beta[j][j];
    }
  }

  /* stampa */
  printf("L\n=============================\n");
  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      printf("%+f ", alpha[i][j]);
    }
    printf("\n");
  }
  printf("\n");

  printf("U\n=============================\n");
  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      printf("%+f ", beta[i][j]);
    }
    printf("\n");
  }
  printf("\n");


  /* verifica */
  printf("I\n============================\n");
  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      res[i][j] = 0;
      for (k = 0; k < N; k++) {
	res[i][j] += alpha[i][k] * beta[k][j];
      }
      printf("%+f ", res[i][j] - m[i][j]);
    }
    printf("\n");
  }

}