// test_qrd.java  QR decomposition using Jama

import Jama.Matrix;
import Jama.QRDecomposition;

import java.text.*;
import java.io.*;

public class test_qrd
{
  public static void main(String[] args)
  { 
    System.out.println("test_qrd.java running ");
    System.out.println("test_qrd.java using Jama QRD ");

    System.out.println("test mat_mul ");
    double tst[][] = {{1.0, 1.0, 1.0},{1.0, 2.0, 1.0},{1.0, 1.0, 2.0}};
    double tstr[][] = {{ -1.732, -2.309, -2.309},
		       {  0.000, -0.816,  0.408},
                       {  0.000,  0.000, -0.707}};
    double tstq[][] = {{-0.577,  0.408,  0.707},
	               {-0.577, -0.816, -0.000},
	               {-0.577,  0.408, -0.707}};

    System.out.println("mat_mul test ");
    print("tst",tst);
    print("tstr",tstr);
    print("tstq",tstq);
    mat_mul(tstq,tstr,tst);
    System.out.println("check tstq * tstr = tst ");
    print("tst",tst);
    System.out.println("end mat_mul test ");
    System.out.println(" ");

	 
    System.out.println("MM = ");
    double[][] vals =
            {{27.7089,   -3.3344,  -23.2704,  -17.2821},
             {27.5579,   -2.1332,  -26.0760,  -17.8489},
             { 3.6486,    0.1602,    0.1226,   -3.2913},
             {28.0106,   -4.5348,  -27.2590,  -15.6983}};
    Matrix MM = new Matrix(vals);
    MM.print(4, 4);

    QRDecomposition QRMM = MM.qr();
    Matrix RMM = QRMM.getR();
    Matrix QMM = QRMM.getQ();

    System.out.println("RMM =");
    RMM.print(4, 4);
    System.out.println("QMM =");
    QMM.print(4, 4);
    Matrix tmp = Matrix.random(4, 4); // over written
    tmp = QMM.times(RMM);
    System.out.println("check MM = QMM*RMM =");
    tmp.print(4, 4);
    System.out.println(" ");
    Matrix QT = Matrix.random(4, 4); // over written
    QT = QMM.transpose();
    System.out.println("QT = transpose(QMM) =");
    QT.print(4,4);
    System.out.println(" ");
    System.out.println("check = QMM*QT =");
    tmp = QMM.times(QT);
    tmp.print(4,4);
    System.out.println("check for I ");
    System.out.println(" ");
    
    Matrix RT = Matrix.random(4, 4); // over written
    tmp = RMM.times(MM);
    tmp = tmp.times(RT);
    System.out.println("check = RMM*MM*RT =");
    tmp.print(4,4);
    System.out.println(" ");
    
    tmp = RMM.times(QMM);
    tmp = tmp.times(RT);
    System.out.println("check = RMM*QMM*RT =");
    tmp.print(4,4);
    System.out.println(" ");
    System.out.println(" ");

    System.out.println("MM=");
    System.out.println("[[ 27.7089+0.j  -3.3344+0.j -23.2704+0.j -17.2821+0.j]");
    System.out.println(" [ 27.5579+0.j  -2.1332+0.j -26.076 +0.j -17.8489+0.j]");
    System.out.println(" [  3.6486+0.j   0.1602+0.j   0.1226+0.j  -3.2913+0.j]");
    System.out.println(" [ 28.0106+0.j  -4.5348+0.j -27.259 +0.j -15.6983+0.j]]");
    System.out.println("eigen values EMM=");
    System.out.println("[1.00073787  4.00011124  2.99911203  2.00003886]");
    System.out.println("eigen vectors are columns ZMM=");
    System.out.println("[[0.53321094  0.71095368   0.57067791  0.71314171]");
    System.out.println(" [0.59282242  0.360302     0.5751235   0.11096673]");
    System.out.println(" [0.08308547  0.59472527   0.09388313  0.31375614]");
    System.out.println(" [0.59778297  0.10501993   0.57857207  0.61698652]]");
    System.out.println("end ZMM");
    System.out.println(" ");
    System.out.println(" ");
 
    System.out.println("MPD = ");
    double[][] mpdv =
            {{ 1.6886,  0.0430,  0.4649, -0.0103},
            { 0.0430,  3.0672,  0.0291, -0.2544},
            { 0.4649,  0.0291,  1.3139, -0.0069},
            {-0.0103, -0.2544, -0.0069,  3.9303}};
    Matrix MPD = new Matrix(mpdv);
    MPD.print(4, 4);

    QRDecomposition QRMPD = MPD.qr();
    Matrix RMPD = QRMPD.getR();
    Matrix QMPD = QRMPD.getQ();

    System.out.println("RMPD =");
    RMPD.print(4, 4);
    System.out.println("QMPD =");
    QMPD.print(4, 4);
    System.out.println(" ");
    System.out.println(" ");


    System.out.println("M = ");
    double[][] m =
            {{ 1.0,  5.0,  6.0,  7.0},
             { 0.0,  2.0,  7.0,  6.0},
             { 0.0,  0.0,  3.0,  5.0},
             {-0.0,  0.0,  0.0,  4.0}};
    Matrix M = new Matrix(m);
    M.print(4, 4);

    QRDecomposition QRM = M.qr();
    Matrix R = QRM.getR();
    Matrix Q = QRM.getQ();

    // System.out.println("QRM =");
    // QRM.print(4, 4);
    System.out.println("R =");
    R.print(4, 4);
    System.out.println("Q =");
    Q.print(4, 4);
    System.out.println(" ");
    System.out.println(" ");

    System.out.println("M3a = ");
    double[][] m3a =
            {{ 1.0,  1.0,  1.0},
             { 1.0,  2.0,  1.0},
             { 1.0,  1.0,  2.0}};
    Matrix M3a = new Matrix(m3a);
    M3a.print(3, 3);
    QRDecomposition QRM3a = M3a.qr();
    Matrix R3a = QRM3a.getR();
    Matrix Q3a = QRM3a.getQ();
    System.out.println("R3a =");
    R3a.print(3, 3);
    System.out.println("Q3a =");
    Q3a.print(3, 3);
    System.out.println(" ");
    Matrix tmp3 = Matrix.random(3, 3); // over written
    tmp3 = Q3a.times(R3a);
    System.out.println("check M3a = Q3a*R3a =");
    tmp3.print(3, 3);
    System.out.println(" ");
    System.out.println(" ");


    System.out.println("M3b = ");
    double[][] m3b =
            {{ 2.0,  1.0,  1.0},
             { 1.0,  2.0,  1.0},
             { 1.0,  1.0,  2.0}};
    Matrix M3b = new Matrix(m3b);
    M3b.print(3, 3);
    QRDecomposition QRM3b = M3b.qr();
    Matrix R3b = QRM3b.getR();
    Matrix Q3b = QRM3b.getQ();
    System.out.println("R3b =");
    R3b.print(3, 3);
    System.out.println("Q3b =");
    Q3b.print(3, 3);
    System.out.println(" ");
    tmp3 = Q3b.times(R3b);
    System.out.println("check M3b = Q3b*R3b =");
    tmp3.print(3, 3);
    System.out.println(" ");
    System.out.println(" ");


    System.out.println("M3c = ");
    double[][] m3c =
          {{ (1.0/3.0),  (-2.0/3.0),  ( 2.0/3.0)},
           { (2.0/3.0),  (-1.0/3.0),  (-2.0/3.0)},
	   { (2.0/3.0),  ( 2.0/3.0),  ( 1.0/3.0)}};
    Matrix M3c = new Matrix(m3c);
    M3c.print(3, 3);
    QRDecomposition QRM3c = M3c.qr();
    Matrix R3c = QRM3c.getR();
    Matrix Q3c = QRM3c.getQ();
    System.out.println("R3c =");
    R3c.print(3, 3);
    System.out.println("Q3c =");
    Q3c.print(3, 3);
    System.out.println(" ");
    tmp3 = Q3c.times(R3c);
    System.out.println("check M3c = Q3c*R3c =");
    tmp3.print(3, 3);
    System.out.println(" ");
    System.out.println(" ");

    // create a symmetric positive definite matrix
    System.out.println("create a symmetric positive definite matrix AA");
    int N = 4;
    Matrix A = Matrix.random(N, N);
    Matrix AA = Matrix.random(N, N); // over written
    System.out.println("random A =");
    A.print(4,4);
    AA = A.transpose().times(A);
    System.out.println("AA = A*A^t =");
    AA.print(4,4);

    System.out.println("A3 = ");
    double[][] vA3 =
            {{12.0,   -51.0,    4.0},
             { 6.0,   167.0,  -68.0},
             {-4.0,    24.0,  -41.0}};
    Matrix A3 = new Matrix(vA3);
    A3.print(3, 3);

    System.out.println("web Q wQ3 = ");
    double [][] wQ =
            {{ 6.0/7.0, -69.0/175.0, -58.0/175.0},
             { 3.0/7.0, 158.0/175.0,   6.0/175.0},
             {-2.0/7.0,   6.0/35.0,  -33.0/35.0 }};
    Matrix wQ3 = new Matrix(wQ);
    wQ3.print(3,3);

    System.out.println("web R wR3 = ");
    double [][] wR =
             {{ 14.0,  21.0, -14.0 },
              {  0.0, 175.0, -70.0 },
              {  0.0,   0.0,  35.0 }};
    Matrix wR3 = new Matrix(wR);
    wR3.print(3,3);

    System.out.println("web A = wQ3*wR3 = ");
    Matrix wA3 = Matrix.random(3, 3); // over written
    wA3 = wQ3.times(wR3);
    wA3.print(3,3);

    System.out.println("given A[][] = | a0  a1  a2 |");
    System.out.println("a0 = { 12,   6,  -4}   column");
    System.out.println("a1 = {-51, 167,  24}   column");
    System.out.println("a2 = {  4, -68, -41}   column");
    System.out.println("vectors aj  uj  ej   all same length ");
    System.out.println("u0 = a0");
    System.out.println("e0 = u0/||u0||");
    System.out.println("u1 = a1 - proj(u0,a1)");
    System.out.println("e1 = u1/||u1||");
    System.out.println("u2 = a2 - proj(u0,a2) - proj(u1,a2");
    System.out.println("e2 = u2/||u2||");
    System.out.println("uk = ak - sumj=0..k(proj(uj,ak)");
    System.out.println("ek = uk/||k2||");
    
    System.out.println("wQ = | e0 e1 e2 |");
    System.out.println("wR = | <e0,a0>  <e0,a1>  <e0,a2> |");
    System.out.println("     | 0        <e1,a1>  <e1,a2> |");
    System.out.println("     | 0        0        <e2,a2> |");
    System.out.println(" ");
    System.out.println("v = proj(u,a) = <u,a>/<u,u> u");
    System.out.println("<u,a> = sumj=0..n(u[j]*a[j])");
    System.out.println(" ");

    System.out.println("compute myQ,myR from A3 myA ");
    double myQ[][] = new double[3][3]; 
    double myR[][] = new double[3][3];
    double myA[][] = {{1.0, 1.0, 1.0},{1.0, 2.0, 1.0},{1.0, 1.0,2.0}}; 
    print("myA",myA);
    qr_from_a(myA, myQ, myR);
    System.out.println("myQ = ");
    print("myQ",myQ);
    System.out.println("myR = ");
    print("myR",myR);
    System.out.println(" ");
    mat_mul(myQ,myR,myA);
    System.out.println("check myA = myQ*myR =");
    print("myA",myA);
    System.out.println("end compute myQ,myR from A3 ");
    System.out.println(" ");
    
    System.out.println("using Jama to get Q R");
    QRDecomposition QRA3 = A3.qr();
    Matrix RA3 = QRA3.getR();
    Matrix QA3 = QRA3.getQ();

    System.out.println("RA3 =");
    RA3.print(3, 3);
    System.out.println("QA3 =");
    QA3.print(3, 3);
    Matrix tA3 = Matrix.random(3, 3); // over written
    tA3 = QA3.times(RA3);
    System.out.println("check A3 = QA3*RA3 =");
    tA3.print(3, 3);
    System.out.println(" ");
    Matrix RA3T = Matrix.random(3, 3); // over written
    RA3T = RA3.transpose();
    System.out.println("RA3T = transpose(RA3) =");
    RA3T.print(3,3);
    System.out.println(" ");
    System.out.println("check = RA3*RA3T =");
    tA3 = RA3.times(RA3T);
    tA3.print(3,3);
    System.out.println(" "); 
    System.out.println("v = |  0.3281   -0.9905    0.2548 |");
    System.out.println("    | -0.9369    0.0872    0.3028 |");
    System.out.println("    | -0.1207    0.1061    0.9184 |");
    System.out.println("e =  156.1367   16.0600  -34.1967  ");
    System.out.println(" "); 
	
    System.out.println("end test_qrd.java ");
  } // end main

  static double [][] transpose(double M[][])
  {
    int row = M.length;
    int col = M[0].length;
    double T[][] = new double[col][row];
    for(int j=0; j<col; j++)
    {
      for(int i=0; i<row; i++)
      {
        T[j][i] = M[i][j];
      } // end i
    } // end j
    return T;
  } // end transpose

  static double  len(double a[]) // ||a||
  {
    int n = a.length;
    double sum = 0.0;
    for(int i=0; i<n; i++)
    {
      sum += a[i]*a[i];
    } // end i
    return Math.sqrt(sum);
  } // end len

  static double  inner_product(double a[], double u[]) // <a,u>
  {
    int n = a.length;
    double ip = 0.0;
    for(int i=0; i<n; i++)
    {
      ip += a[i]*u[i];
    } // end i
    return ip;
  } // end inner_product

  static double [] proj(double a[], double u[]) // <u,a>/<u,u> * u
  {
    int n = a.length;
    double v[] = new double[n];
    double ipuu = inner_product(u,u);
    double ipua = inner_product(u,a);
    double r = ipua/ipuu;
    for(int i=0; i<n; i++)
    {
      v[i] = r*u[i];
    } // end i
    return v;
  } // end proj  projection

  static void  vec_sub(double a[], double b[], double c[])
  {
    int n = a.length;
    for(int i=0; i<n; i++)
    {
      c[i] = a[i]-b[i];
    } // end i
  } // end vec_sub


  static void  col_to_vec(double A[][], int col, double V[])
  {
    int n = A.length;
    for(int i=0; i<n; i++)
    {
      V[i] = A[i][col];
    } // end i
  } // end col_to_vec

  static void  vec_to_col(double V[], int col, double A[][])
  {
    int n = A.length;
    for(int i=0; i<n; i++)
    {
      A[i][col] = V[i];
    } // end i
  } // end vec_to_col

  static void  qr_from_a(double A[][], double Q[][], double R[][])
  {
    int n = A.length;
    double u[] = new double[n];
    double uj[] = new double[n]; // temp
    double U[][] = new double[n][n]; // saved copies
    double e[] = new double[n];
    double E[][] = new double[n][n]; // saved copies
    double a[] = new double[n];
    double ai[] = new double[n]; // temp
    double ulen = 0.0;
    double sub = 0.0;

    System.out.println("qr_from_a running");
    for(int i=0; i<n; i++)
    {
      System.out.println("i="+i);
      System.out.println("compute U");
      if(i==0)
      {
        col_to_vec(A,i,u);
        vec_to_col(u,i,U); // save u0
      } // end if
      if(i>0)
      {
        // ui = ai - proj(u0,ai) - proj(uj,ai) ...
        col_to_vec(A,i,u);
        for(int j=0; j<i; j++)
        {
          col_to_vec(A,i,ai);
          System.out.println("for ai i="+i);
          col_to_vec(U,j,uj);
          System.out.println("for uj j="+j);
          vec_sub(u,proj(uj,ai),u);
        } // end j
        vec_to_col(u,i,U); // save ui
      } // end if
      printv("u",u);
      ulen = len(u);
      System.out.println("ulen="+ulen);
      for(int j=0; j<n; j++)
      {
        e[j] = u[j]/ulen;
        Q[i][j] = e[j];
      } // end j
      printv("e",e);
    } // end i
    System.out.println("U and E in Q built");
    // compute R from inner_products   Q is E 
    for(int i=0; i<n; i++)
    {
      sub = 0.0;
      for(int j=0; j<n; j++) // <ei,aj>
      {
        R[i][j] = 0.0;
        if(j>=i)
        {
          col_to_vec(A,j,a);
          col_to_vec(Q,i,e);
          R[i][j] = inner_product(e,a);
        } // end if
      } // end i
    } // end i
    System.out.println("qr_from_a finished");
  } // end qr_from_a

  static void print(String name, double M[][])
  {
    int row = M.length;
    int col = M[0].length;
    System.out.println(name+"["+row+"]["+col+"]=");
    for(int i=0; i<row; i++)
    {
      System.out.print("[");
      for(int j=0; j<col; j++)
      {
        System.out.print(M[i][j]+" ");
      } // end j
      System.out.println("]");
    } // end i
    System.out.println(" ");
  } // end print

  static void printv(String name, double V[])
  {
    int n = V.length;
    System.out.println(name+"["+n+"]=");
    System.out.print("[");
    for(int i=0; i<n; i++)
    {
      System.out.print(V[i]+" ");
    } // end i
    System.out.println("]");
    System.out.println(" ");
  } // end printv

  static void mat_mul(double A[][], double B[][], double C[][])
  {
    int n = A.length;
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
	C[i][j] = 0.0;
        for(int k=0; k<n; k++)
        {
          C[i][j] = C[i][j] + A[i][k] * B[k][j];
        } // end k
      } // end j
    } // end i
  } // end mat_mul
    
} // test_qrd.java