package ejemplos;

/**
 * <p>Title: Sistemas lineales</p>
 * <p>Description: Resuelve sistemas de ecuaciones Ax = b</p>
 * <p>Copyright: Copyright (c) 2004</p>
 * <p>Company: UMSNH</p>
 * @author Dr. Felix Calderon Solorio
 * @version 1.0
 */

public class slineal {

  static public void SolucionLU(double A[][], double x[], double b[]) {
    LU(A);
    sustitucion_hacia_delante(A, x, b);
    sustitucion_hacia_atras(A, b, x);
    for(int i=0; i<x.length; i++)
      x[i] = b[i];
  }

  static public void Solucion(double A[][], double x[], double b[]) {
    EliminacionGaussiana(A, x, b);
    sustitucion_hacia_atras(A, x, b);
  }

  static public void EliminacionGaussiana(double a[][], double x[], double b[]) {
    int n = a.length;

    for (int k = 0; k <= n - 2; k++) {
      for (int i = k + 1; i <= (n - 1); i++) {
        b[i] -= a[i][k] * b[k] / a[k][k];
        for (int j = n - 1; j >= k; j--)
          a[i][j] -= a[i][k] * a[k][j] / a[k][k];
      }
      /*
      for(int i=0; i<n; i++)
      {
        for (int j = 0; j < n; j++)
          System.out.print(a[i][j] + " ");

        System.out.println(" " + b[i]);
      }*/

    }
  }

  static public void LU(double a[][]) {
    int n = a.length;
    int i, j, k;
    double suma;

    for (i = 0; i < n; i++) {

      for (j = 0; j <= i - 1; j++) {
        suma = 0;
        for (k = 0; k <= j - 1; k++)
          suma += a[i][k] * a[k][j];
        a[i][j] = (a[i][j] - suma) / a[j][j];
      }

      for (j = i; j < n; j++) {
        suma = 0;
        for (k = 0; k <= i - 1; k++)
          suma += a[i][k] * a[k][j];
        a[i][j] = a[i][j] - suma;
      }
    }
  }

  static public void sustitucion_hacia_atras(double a[][], double x[], double b[]) {
    double suma;
    int k, j;
    int n = x.length;

    for (k = n - 1; k >= 0; k--) {
      suma = 0;
      for (j = k + 1; j < n; j++)
        suma += a[k][j] * x[j];
      x[k] = (b[k] - suma) / a[k][k];
    }
  }

  static public void sustitucion_hacia_delante(double A[][], double y[],
                                               double b[]) {
    int k, i, j, n;
    double suma;

    n = A.length;
    for (k = 0; k < n; k++) {
      suma = 0;
      for (j = 0; j < k; j++) {
        suma += A[k][j] * y[j];
      }
      y[k] = b[k] - suma;
    }
  }
}