package ppmrmn;

public class Jacobi extends RealMat {

  private int p;
  private int q;
  private RealMat matR = new RealMat();
  private RealMat matA = new RealMat();

  public Jacobi() {
  }

  public RealMat getMatR() {
    return matR;
  }

  public int sign (double x) {
    if ( x < 0) return -1;
    return 1;
  }

  //Final matrix matA is diagonal matrix of eiganvalues
  //Final matrix matR has eigenvectors in rows
  public RealMat jacRot(RealMat matADiagonaliser) {
    for (int i = 0; i < n; i++)
      for (int j = 0; j < n; j++)
        matA.re[i][j] = matADiagonaliser.re[i][j];
    int iter = 0;

    for (int i = 0; i < n; i++)
      for (int z = 0; z < n; z++)
        matR.re[i][z] = 0.0;   // RAZ

    for (int i = 0; i < n; i++)
      matR.re[i][i] = 1.0;   // unit matrix

    String etat = "ITERATING";
    double tolerance = 0.00001;
    int maxIter = 100;

    while (etat == "ITERATING") {
      double max = maxOffDiag();
      if (max > tolerance) {
        rotate();
        if (++iter == maxIter) {
          etat = "WONSTOP";
          System.out.println("maximum iteration atteinte");
        }
      }
      else {
        etat = "SUCCESS";
//        matR = RealMat.transpose(matR); //eigenvectors in rows
      }
    }
    return matA;
  }

  private double maxOffDiag() {
    double max = 0.0;
    for (int i = 0; i < n - 1; i++) {
      for (int j = i + 1; j < n; j++) {
        double aij = Math.abs(matA.re[i][j]);
        if (aij > max) {
          max = aij;
          p = i;
          q = j;
        }
      }
    }
    return max;
  }

  //Rotates matrix matA through theta in pq-plane to set matA.re[p][q] = 0
  //Rotation stored in matrix matR whose columns are eigenvectors of matA
  private void rotate(){
    // d = cot 2*theta, t = tan theta, c = cos theta, s = sin theta
    double d = (matA.re[p][p] - matA.re[q][q])/(2.0*matA.re[p][q]);
    double t = sign(d)/(Math.abs(d) + Math.sqrt(d*d + 1));
    double c = 1.0/Math.sqrt(t*t + 1);
    double s = t*c;
    matA.re[p][p] += t*matA.re[p][q];
    matA.re[q][q] -= t*matA.re[p][q];
    matA.re[p][q] = matA.re[q][p] = 0.0;
    for (int k = 0; k < n; k++) {  // Transform matA
      if (k != p && k != q) {
        double akp = c*matA.re[k][p] + s*matA.re[k][q];
        double akq = -s*matA.re[k][p] + c*matA.re[k][q];
        matA.re[k][p] = matA.re[p][k] = akp;
        matA.re[k][q] = matA.re[q][k] = akq;
      }
    }
    for (int k = 0; k < n; k++) {  // Store matR.re
        double rkp = c*matR.re[k][p] + s*matR.re[k][q];
        double rkq = -s*matR.re[k][p] + c*matR.re[k][q];
        matR.re[k][p] = rkp;
        matR.re[k][q] = rkq;
    }
  }
}