package operaciones;

import javax.swing.*;
import java.io.*;
import java.util.*;

/**
 * <p>Title: Objeto Matriz ver 1.01</p>
 * <p>Description: Define las operaciones basicas con matrices</p>
 * <p>Copyright: Copyright (c) 2005</p>
 * <p>Company: UMSNH</p>
 * @author Dr. Felix Calderon Solorio
 * @version 1.01
 */

public class Matriz
{
  public int nren, ncol;
  public double datos[][];

  /**
   *  Crea una matriz nula
   *
  */

  public Matriz() {
    nren = 0;
    ncol = 0;
  }

  /**
   * Crea una matriz con n renglones y m columnas de puros ceros
   * @param n numero de renglones
   * @param m numero de columnas
   *
   */

  public Matriz(int n, int m) {
    inicializa(n, m);
  }

  /**
   * Crea una matriz con n renglones y m columnas y pone el valor dado
   * @param n numero de renglones
   * @param m numero de columnas
   * @param val valor para inicializar
   *
   */

  public Matriz(int r, int c, double val) {
    nren = r;
    ncol = c;

    datos = new double[nren][ncol];
    for (int i = 0; i < r; i++)
      for (int j = 0; j < c; j++)
        datos[i][j] = val;
  }
  /**
   * Identidad: Crea una matriz de n*n con unos en la diagonal
   * @param r Tamano de la matriz identidad
   */

  public Matriz Identidad(int r)
  {
    inicializa(r, r);
    for(int i=0; i<r; i++)
      this.datos[i][i] = 1;
    return this;
  }

  /**
   * Inicializa las variables para los constructores.
   * @param r int
   * @param c int
   */

  private void inicializa(int r, int c) {
    // Si la matriz es nula reserva memoria.
    if ( (nren == 0 && ncol == 0) || this == null) {
      nren = r;
      ncol = c;

      datos = new double[nren][ncol];
      for (int i = 0; i < r; i++)
        for (int j = 0; j < c; j++)
          datos[i][j] = 0;
    }
  }

  /**
   * Crea una matriz a partir de informacion en un archivo.
   * @param archivo String
   */

  private void Dimensiones(String archivo) {
    StringTokenizer st;
    String linea;
    int num_ren = 0, num_col = 0;

    try {
      RandomAccessFile DIS = new RandomAccessFile(archivo, "r");

      while ( (linea = DIS.readLine()) != null) {
        st = new StringTokenizer(linea);
        if(st.countTokens() != 0) num_ren ++;

        if(num_col == 0)
        {
          while (st.hasMoreTokens())
          {
            num_col++;
            st.nextToken();
          }
        }
      }
      this.nren = num_ren;
      this.ncol = num_col;
      DIS.close();
    }
    catch (IOException e) {
      JOptionPane.showMessageDialog(null, "Error" + e.toString(), "ERROR",
                                    JOptionPane.ERROR_MESSAGE);
    }
  }

  /**
   * Contructor de la clase matriz que resibe como dato la ruta de un archivo
   * texto
   * @param archivo Nombre y ruta del archivo de texto
   */

  public Matriz(String archivo) {
    String linea = "", dato = "";
    int i = 0, j = 0;
    StringTokenizer st;
    double x;

    Dimensiones(archivo) ;
    datos = new double[nren][ncol];

    try {

      RandomAccessFile DIS = new RandomAccessFile(archivo, "r");

      i = 0;
      while ( ( (linea = DIS.readLine()) != null) && i < nren) {
        st = new StringTokenizer(linea);
        if(st.countTokens() != 0)
        {
          st = new StringTokenizer(linea);
          for (j = 0; j < ncol; j++) {
            dato = st.nextToken();
            x = Double.valueOf(dato).doubleValue();
            this.datos[i][j] = x;
          }
          i++;
        }
      }
      DIS.close();
    }
    catch (IOException e) {
      JOptionPane.showMessageDialog(null, "Error" + e.toString(), "ERROR",
                                    JOptionPane.ERROR_MESSAGE);
    }
  }

  /**
   * Inicializa la matriz con un arreglo bidimensional
   *
   * @param d arreglo de dobles con la informacion a almacenar en la
   * matriz
   *
   */

  public Matriz(double d[][]) {
    // crea una matriz con n renglones y m columnas

    nren = d.length;
    ncol = d[0].length;

    datos = new double[nren][ncol];
    for (int i = 0; i < nren; i++)
      for (int j = 0; j < ncol; j++)
        datos[i][j] = d[i][j];
  }

  /**
   * Inicializa la matriz con un arreglo unidimensional
   *
   * @param d arreglo de dobles con la informacion a almacenar en la
   * matriz
   *
   */

  public Matriz(double d[]) {
    // crea una matriz con n renglones y 1 columna

    nren = d.length;
    ncol = 1;

    datos = new double[nren][ncol];
    for (int i = 0; i < nren; i++)
        datos[i][0] = d[i];
  }

  /**
   * Inicializa la matriz con un arreglo unidimensional
   *
   * @param r renglones de la matriz
   * @param c columna de la matriz
   * @param d arreglo de dobles con la informacion a almacenar en la
   * matriz
   *
   */

  public Matriz(int r, int c, double d[]) {
    // crea una matriz con n renglones y m columnas dado una lista de datos

    nren = r;
    ncol = c;
    int k = 0;

    datos = new double[nren][ncol];
    for (int i = 0; i < nren; i++)
      for(int j = 0; j< ncol; j++)
        datos[i][j] = d[k++];
  }


  /**
   * Imprime el contenido en una matriz
   */

  public void imprime() {
    int i, j;
    double aux;
    if ( (nren == 0) && (ncol == 0)) {
      System.out.println("No tiene informacion la matriz");
      return;
    }

    for (i = 0; i < this.nren; i++) {
      for (j = 0; j < this.ncol; j++)
        System.out.print(redondea(datos[i][j], 7) + " \t");
      System.out.println(" ");
    }
  }

  /**
   * Rutina para redondear un numero
   * @param dato valor a redondear
   * @param max tamano del numero a imprimir
   */

  String redondea(double dato, int max) {
    String salida = "", aux = "";
    float a;

    a = ( Math.round(dato * 10000)) / 10000.0f;

    aux += a;
    max -= aux.length();
    for (int i = 0; i < max; i++)
      salida += " ";
    salida += aux;
    return salida;
  }

  /**
   * Algoritmo para suma de matrices
   *
   * @param b Matriz a sumar
   * @return resul suma de la matriz que hace el llamado y  b
   */

  public Matriz mas(Matriz b) {
    if ( (this.nren != b.nren) || (this.ncol != b.ncol))
      return null;

    Matriz resul = new Matriz(b.nren, b.ncol);

    for (int i = 0; i < b.nren; i++)
      for (int j = 0; j < b.ncol; j++)
        resul.datos[i][j] = this.datos[i][j] + b.datos[i][j];

    return resul;
  }

  /**
   * Algoritmo para sumar una constante a una matriz
   *
   * @param b conatante a sumar
   * @return resul suma
   */

  public Matriz mas(double b) {

    Matriz resul = new Matriz(this.nren, this.ncol);

    for (int i = 0; i < resul.nren; i++)
      for (int j = 0; j < resul.ncol; j++)
        resul.datos[i][j] = this.datos[i][j] + b;

    return resul;
  }

  /**
   * Algoritmo para resta de matrices
   * @param b Matriz sustraendo
   * @return resul resta de la matriz que hace el llamado y b
   */

  public Matriz menos(Matriz b) {
    if ( (this.nren != b.nren) || (this.ncol != b.ncol))
      return null;

    Matriz resul = new Matriz(b.nren, b.ncol);

    for (int i = 0; i < b.nren; i++)
      for (int j = 0; j < b.ncol; j++)
        resul.datos[i][j] = this.datos[i][j] - b.datos[i][j];

    return resul;
  }

  /**
   * Algoritmo para resta una constante a una matriz
   * @param b sustraendo
   * @return resul resta de la matriz que hace el llamado y b
   */

  public Matriz menos(double b) {

    Matriz resul = new Matriz(this.nren, this.ncol);

    for (int i = 0; i < resul.nren; i++)
      for (int j = 0; j < resul.ncol; j++)
        resul.datos[i][j] = this.datos[i][j] - b;

    return resul;
  }


  /**
   * Algoritmo para multiplicacion de matrices
   * @param b Matriz por la que se multiplica
   * @return resul producto de la matriz que hace el llamado y la matriz b
   */

  public Matriz por(Matriz b) {
    int i, j;

    double suma;

    if (this.ncol == b.nren) {
      Matriz resul = new Matriz(this.nren, b.ncol);
      for (i = 0; i < resul.nren; i++)
        for (j = 0; j < resul.ncol; j++) {
          suma = 0;
          for (int k = 0; k < this.ncol; k++)
            suma += this.datos[i][k] * b.datos[k][j];
          resul.datos[i][j] = suma;

        }
      return resul;
    }
    else
      System.out.println("Matrices con diferente tamano");
    return null;
  }

  /**
   * Algoritmo para multiplicacion de una matricez por una constante
   * @param b constante
   * @return resul producto de la matriz que hace el llamado y la cte b
   */

  public Matriz por(double b) {

    Matriz resul = new Matriz(this.nren, this.ncol);

    for (int i = 0; i < resul.nren; i++)
      for (int j = 0; j < resul.ncol; j++)
        resul.datos[i][j] = this.datos[i][j] * b;

    return resul;
  }

  /**
   * Eleva una matriz a un exponente dado
   * @param b double Exponente
   * @return Matriz
   */


  public Matriz pow(double b) {

    Matriz resul = new Matriz(this.nren, this.ncol);

    for (int i = 0; i < resul.nren; i++)
      for (int j = 0; j < resul.ncol; j++)
        resul.datos[i][j] = Math.pow(this.datos[i][j],  b);

    return resul;
  }


  /**
   *  Algoritmo para encontrar la transpuesta de un matriz
   * @return resul matriz transpuesta
   */

  public Matriz T() {
    Matriz resul = new Matriz(this.ncol, this.nren);

    for (int i = 0; i < resul.nren; i++)
      for (int j = 0; j < resul.ncol; j++)
        resul.datos[i][j] = this.datos[j][i];
    return resul;
  }

  /**
   * Algoritmo de Sustitucion hacia atras
   * @param a Matriz cuadrada triangular superior
   * @param b Vector de terminos independientes
   * @return resul solucion del sistema triangulas superior ax = b
   */

  public Matriz sustitucion_hacia_atras(Matriz a, Matriz b) {
    if ( (a.nren == a.ncol) && (b.ncol == 1)) {
      Matriz resul = new Matriz(b.nren, b.ncol);

      double suma;
      int k, j;
      int n = this.nren;

      for (k = n - 1; k >= 0; k--) {
        suma = 0;
        for (j = k + 1; j < n; j++)
          suma += (a.datos[k][j]) * (resul.datos[j][0]);
        resul.datos[k][0] = ( (b.datos[k][0]) - (suma)) / (a.datos[k][k]);
      }
      return resul;
    }
    else
      System.out.println("Las dimenciones de las matrices no son adecuadas");
    return null;
  }

  /*******************************************************
    Algoritmo de Sustitucion hacia atras
    para resolver sistemas L'DL o L'L.
   ********************************************************/

  public void sustitucion_hacia_atras2(Matriz a, Matriz b) {
    if ( (a.nren == a.ncol) && (b.ncol == 1)) {
      this.inicializa(a.ncol, 1);

      double suma;
      int k, j;
      int n = this.nren;

      for (k = n - 1; k >= 0; k--) {
        suma = 0;
        for (j = k + 1; j < n; j++)
          suma += (a.datos[j][k]) * (this.datos[j][0]);
        this.datos[k][0] = b.datos[k][0] / a.datos[k][k] - suma;
      }
    }
    else
      System.out.println("Las dimenciones de las matrices no son adecuadas");
  }

  public void sustitucion_hacia_atras3(Matriz a, Matriz b) {
    if ( (a.nren == a.ncol) && (b.ncol == 1)) {
      this.inicializa(a.ncol, 1);

      double suma;
      int k, j;
      int n = this.nren;

      for (k = n - 1; k >= 0; k--) {
        suma = 0;
        for (j = k + 1; j < n; j++)
          suma += (a.datos[j][k]) * (this.datos[j][0]);
        this.datos[k][0] = (b.datos[k][0] - suma) / a.datos[k][k];
      }
    }
    else
      System.out.println("Las dimenciones de las matrices no son adecuadas");
  }

  /**************************************************************
   *         Sustitucion hacia adelante
   *
   ***********************************************************/

  public void sustitucion_hacia_adelante(Matriz A, Matriz b) {
    int k, i, j, n;
    double suma;

    if ( (A.nren == A.ncol) && (b.ncol == 1)) {
      this.inicializa(A.ncol, 1);

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

  public void sustitucion_hacia_adelante3(Matriz A, Matriz b) {
    int k, i, j, n;
    double suma;

    if ( (A.nren == A.ncol) && (b.ncol == 1)) {
      this.inicializa(A.ncol, 1);

      n = A.nren;
      for (k = 0; k < n; k++) {
        suma = 0;
        for (j = 0; j < k; j++) {
          suma += A.datos[k][j] * this.datos[j][0];
        }
        this.datos[k][0] = (b.datos[k][0] - suma) / A.datos[k][k];
      }
    }
  }

  public Matriz LU()
  {
      if(this.nren != this.ncol) return null;

      int i, j, k, N;

      Matriz res = new Matriz(this.nren, this.ncol);
      float suma;

      N = this.nren;

      for(i=0; i<N; i++)
          for(j=0; j<N; j++)
              res.datos[i][j] = this.datos[i][j];


      for (i = 0; i<N; i++) {
          for (j = 0; j <= i - 1; j++) {
              suma = 0;
              for (k = 0; k <= j - 1; k++)
                  suma += res.datos[i][k] * res.datos[k][j];
              res.datos[i][j] = (res.datos[i][j] - suma) / res.datos[j][j];
          }
          for (j = i; j < N; j++) {
              suma = 0;
              for (k = 0; k <= i - 1; k++)
                  suma += res.datos[i][k] * res.datos[k][j];
              res.datos[i][j] = res.datos[i][j] - suma;
          }
      }
      return res;
  }


  /**
   * Soluci�n de un sistema de ecuaciones Ax=b
   * utilizando Eliminacion Gaussiana y Sustitucion hacia atras
   * @param A matriz cuadrada
   * @return resul soluci�n del sistema
   */

  public Matriz entre(Matriz A)
  {
    Matriz Ar=null,  br = null;

    Ar = Matriz.igual_a(A);
    br = Matriz.igual_a(this);

    Matriz.eliminacion_gaussiana(Ar, br);

    Matriz resul = sustitucion_hacia_atras(Ar, br);

    return resul;
  }

  /**
   * Algoritmo de Eliminacion Gaussiana
   * Entrega una matriz diagonal superior
   * @param a Matriz cuadrada
   * @param b Vector de terminos independientes
   */

  public static void eliminacion_gaussiana(Matriz a, Matriz b) {
    if ( (a.nren == a.ncol) && (b.ncol == 1) && (a.nren == b.nren)) {
      int n = a.nren;
      for (int k = 0; k <= n - 2; k++) {
        for (int i = k + 1; i <= (n - 1); i++) {
          b.datos[i][0] -= ( (a.datos[i][k] * b.datos[k][0]) / a.datos[k][k]);
          for (int j = n - 1; j >= k; j--)
            a.datos[i][j] -= ( (a.datos[i][k] * a.datos[k][j]) / a.datos[k][k]);
        }
      }
    }
  }

  static public double determinante(Matriz A)
  {
      Matriz Ar=null;
      double det=1;

      Ar = Matriz.igual_a(A);
      int n = A.nren;
      if(A.nren != A.ncol)
      {
          System.out.println("No se puede calcular el determinante");
          return 0;
      }

      for (int k = 0; k <= n - 2; k++) {
          for (int i = k + 1; i <= (n - 1); i++) {
              for (int j = n - 1; j >= k; j--)
                  Ar.datos[i][j] -=
                          ((Ar.datos[i][k] * Ar.datos[k][j]) / Ar.datos[k][k]);
          }
      }

      for(int k=0; k <n; k++)
          det *= Ar.datos[k][k];

      return det;
  }

  /**
   * Algoritmo de la Inversa de Shiplay para inviertir una matriz
   * return resul A-1
   */

  public Matriz inversa() {
    Matriz resul = new Matriz(this.datos);
    if (resul.nren == resul.ncol && this != null) {
      int n = resul.nren;
      int k, j, i;
      for (k = 0; k < n; k++) {
        for (i = 0; i < n; i++)
          for (j = 0; j < n; j++) {
            if ( (i != k) && (j != k))
              resul.datos[i][j] -= (resul.datos[i][k] * resul.datos[k][j]) /
                  resul.datos[k][k];
          }

        for (j = 0; j < n; j++) {
          if (j != k)
            resul.datos[k][j] = -resul.datos[k][j] / resul.datos[k][k];
        }

        for (i = 0; i < n; i++) {
          if (i != k)
            resul.datos[i][k] = resul.datos[i][k] / resul.datos[k][k];
        }

        resul.datos[k][k] = 1 / resul.datos[k][k];
      }
      return resul;
    }
    else
      System.out.println("no se pudo");
    return null;
  }

  /**
   * Metodo para obtener un dato
   * @param i renglon
   * @param j columna
   * @return dato en la posicion <i,j>
   */

  public double obten(int i, int j)
  {
    return (this.datos[i][j]);
  }

  /**
   * Metodo para poner un dato en la matriz
   * @param i renglon
   * @param j columna
   * @param d dato
   */

  public void inserta(int i, int j, double d)
  {
    this.datos[i][j] = d;
  }

  /**
   * Factorizacion de Cholesky
   */

  public Matriz Cholesky() {
    int i, j, k, n;
    double fact;

    if (this.nren != this.ncol || this == null) {
      System.out.println("Matriz mal condicionada");
      return null;
    }

    n = this.nren;
    for (k = 0; k < n - 1; k++) {
      for (i = k + 1; i < n; i++) {
        fact = this.datos[i][k] / this.datos[k][k];
        this.datos[i][k] = fact;
        for (j = k + 1; j < n; j++)
          this.datos[i][j] = this.datos[i][j] - fact * this.datos[k][j];
      }
      for (j = k + 1; j < n; j++)
        this.datos[k][j] = this.datos[k][j] / this.datos[k][k];
    }
    return this;
  }

  /**
   * Algoritmo de Factorizacion Cholesky L'DL
   */

  public Matriz Cholesky2() {
    int i, j, k, n, s;
    double fact, suma = 0;

    if (this.nren != this.ncol || this == null) {
      System.out.println("Matriz mal condicionada");
      return null;
    }

    n = this.nren;
    for (j = 0; j < n; j++) {
      suma = 0;
      for (s = 0; s <= j - 1; s++)
        suma += this.datos[s][s] * this.datos[j][s] * this.datos[j][s];

      this.datos[j][j] = this.datos[j][j] - suma;
      for (i = j + 1; i < n; i++) {
        suma = 0;
        for (s = 0; s <= j - 1; s++)
          suma += this.datos[s][s] * this.datos[i][s] * this.datos[j][s];
        this.datos[i][j] = this.datos[i][j] - suma;
        this.datos[i][j] = this.datos[i][j] / this.datos[j][j];
      }
    }
    return this;
  }

  public Matriz Cholesky3() {
    int i, j, k, n, s;
    double fact, suma = 0;

    if (this.nren != this.ncol || this == null) {
      System.out.println("Matriz mal condicionada");
      return null;
    }
    n = this.nren;
    for (k = 0; k < n; k++) {
      for (i = 0; i <= k - 1; i++) {
        suma = 0;
        for (j = 0; j <= i - 1; j++)
          suma += this.datos[i][j] * this.datos[k][j];

        this.datos[k][i] = (this.datos[k][i] - suma) / this.datos[i][i];
      }
      suma = 0;
      for (j = 0; j <= k - 1; j++)
        suma += this.datos[k][j] * this.datos[k][j];
      this.datos[k][k] = Math.sqrt(this.datos[k][k] - suma);
    }
    return this;
  }

  /**
   * Metodo de Cholesky incompleto
   */

  public Matriz Cholesky_Incompleto() {
    int i, j, k, n, s;
    double fact, suma = 0;

    if (this.nren != this.ncol || this == null) {
      System.out.println("Matriz mal condicionada");
      return null;
    }
    

    n = this.nren;
    for (k = 0; k < n; k++) {
      for (i = 0; i <= k - 1; i++) {
        if(this.datos[k][i] != 0) {
          suma = 0;
          for (j = 0; j <= i - 1; j++)
            suma += this.datos[i][j] * this.datos[k][j];

          this.datos[k][i] = (this.datos[k][i] - suma) / this.datos[i][i];
        }
      }
      suma = 0;
      for (j = 0; j <= k - 1; j++)
        suma += this.datos[k][j] * this.datos[k][j];
      this.datos[k][k] = Math.sqrt(this.datos[k][k] - suma);
    }
    return this;
  }

  /**
   * Copia una matriz en otra
   * @param A matriz a copiar en la matriz que hace el llamado
   */

  public static Matriz igual_a(Matriz A)
  {
    Matriz resul  = new Matriz(A.nren, A.ncol);
    for(int i=0; i<A.nren; i++)
      for(int j=0; j<A.ncol; j++)
        resul.datos[i][j] = A.datos[i][j];

    return resul;
  }

  /**
   * Gatdiente conjugado precondicionado
   * @param A Matriz
   * @param M Matriz
   * @param x Matriz
   * @param b Matriz
   * @return Matriz
   */

  public Matriz Gradiente_conjugado_pre(Matriz A, Matriz M, Matriz x, Matriz b)
  {
    double alpha, beta, rho, rho1;

    Matriz rr     = new Matriz(1,1);
    Matriz pAp    = new Matriz(1,1);

    Matriz Ap     = new Matriz(A.nren, 1);
    Matriz r      = new Matriz(A.nren, 1);
    Matriz p      = new Matriz(A.nren, 1);
    Matriz aux1   = new Matriz(1, A.nren);
    Matriz aux2   = new Matriz(A.nren, 1);

    Matriz y      = new Matriz(A.nren, 1);

    rho  = 1;
    beta = 0;

    Ap = A.por(x);
    r = b.menos(Ap);

    System.out.println("r = ");
    r.imprime();

    for(int it =1; it <= 16; it++)
    {
      aux2.sustitucion_hacia_adelante3(M,r);
      y.sustitucion_hacia_atras3(M,aux2);

      rho1 = rho;

      rr = (r.T()).por(y);
      rho = rr.datos[0][0];

      if(it ==1) p.igual_a(y);
      else beta = rho/rho1;

      p = y.mas(p.por(beta));

      Ap = A.por(p);
      pAp = (p.T()).por(Ap);

      alpha = rho/pAp.datos[0][0];

      x = x.mas(p.por(alpha));
      r = r.menos(Ap.por(alpha));
      System.out.println("iteracion " + it);
      x.imprime();
    }
    return x;
  }
}