// decomp2.java  given A create B1 and B2 ... such that B1 * B2 ... = A
// from mathoverflow.net
// uses inverse.java 
// uses Matrix.java for utility routines.

public class decomp2
{
  int n = 3;
  double det;
  double A[][] = {{2.0, 9.0, 4.0},
       	          {7.0, 5.0, 3.0},
		  {6.0, 1.0, 8.0}};
  double M21[][] = {{0.0,      0.0, 0.0},  // below (1,1)(1,1)
		    {-7.0/2.0, 0.0, 0.0},  // -A[2][1]/A[1][1]
		    {0.0,      0.0, 0.0}};
  double M31[][] = {{0.0,  0.0, 0.0},
		    {0.0,  0.0, 0.0},
		    {-3.0, 0.0, 0.0}}; // -A[1][3]/A[1][1]
  double E21[][] = new double[n][n];   // eye + M21
  double E31[][] = new double[n][n];   // eye + M31 
  double eye[][] = new double[n][n];
  double E31p[][] = new double[n][n];  // E31 * E21 * A  product
    
  double M32[][] = {{0.0, 0.0,        0.0},  // below (2,2)(2,2)
                    {0.0, 0.0,        0.0},
		    {0.0, -26.0/26.5, 0.0}}; // E31p[1][3]/-E31p[1][2]
  double E32[][] = new double[n][n];   // eye + M32
  double E32p[][] = new double[n][n];  // E32 * E31 * E21 * A

  double M23[][] = {{0.0, 0.0, 0.0},  // above  (3,3)(3,3)
		    {0.0, 0.0, 11.0/6.79245}, // -E32p[3][2]/E32p[3][3]
		    {0.0, 0.0, 0.0}};
  double M13[][] = {{0.0,  0.0, -4.0/6.79245}, // -E32p[3][1]/E32p[3][3] 
		    {0.0,  0.0, 0.0},
		    {0.0,  0.0, 0.0}};
  double E23[][] = new double[n][n];   // eye + M23
  double E13[][] = new double[n][n];   // eye + M13
  double E13p[][] = new double[n][n];  // E13 * E23 * E32 * E31 * E21 * A

  double M12[][] = {{0.0, 9.0/26.5, 0.0},  // above (2,2)(2,2)
		    {0.0, 0.0, 0.0},       // E23p[2][1]/-E23p[2][2]
		    {0.0, 0.0, 0.0}};
  double E12[][] = new double[n][n];   // eye + M12
  double E12p[][] = new double[n][n];  // E12*E13*E23*E32*E31*E21*A

  double D[][] = {{2.0, 0.0,   0.0}, // just diagonal of E23p
		  {0.0, -26.5, 0.0},
		  {0.0, 0.0,   6.79245}};// diagonal

  double E21i[][] = new double[n][n];   // inverse(E21)
  double E31i[][] = new double[n][n];   // inverse(E31)
  double E32i[][] = new double[n][n];   // inverse(E32)
  double E23i[][] = new double[n][n];   // inverse(E23)
  double E13i[][] = new double[n][n];   // inverse(E13)
  double E12i[][] = new double[n][n];   // inverse(E12)
  double T1[][] = new double[n][n];     // temp for product
  double T2[][] = new double[n][n];     // temp for product
  double T3[][] = new double[n][n];     // temp for product
  double T4[][] = new double[n][n];     // temp for product
  double T5[][] = new double[n][n];     // temp for product
  double Ack[][] = new double[n][n];    // check product ...
  double AI[][] = new double[n][n];      // inverse
    
  public decomp2()
  {
    System.out.println("decomp2.java running");


    det = Matrix.determinant(A);
    System.out.println("A 3 x 3, det="+det);
    prt_mat("A",A);

    det = Matrix.determinant(M21);
    System.out.println("M21 below (1,1)(1,1) det="+det);
    prt_mat("M21",M21);

    det = Matrix.determinant(M31);
    System.out.println("M31                  det="+det);
    prt_mat("M31",M31);


    ident(eye);
    det = Matrix.determinant(eye);
    System.out.println("eye   det="+det);
    prt_mat("eye",eye);

    mat_add(eye, M21, E21);
    det = Matrix.determinant(E21);
    System.out.println("E21=eye+M21   det="+det);
    prt_mat("E21",E21);

    mat_add(eye, M31, E31);
    det = Matrix.determinant(E31);
    System.out.println("E31=eye+M31   det="+det);
    prt_mat("E31",E31);
    mat_mul(E31, E21, T1);
    mat_mul(T1, A, E31p);
    det = Matrix.determinant(E31p);
    System.out.println("E31p=E31 * E21 * A   det="+det);
    prt_mat("E31p",E31p);

    mat_add(eye, M32, E32);
    det = Matrix.determinant(E32);
    System.out.println("E32=eye+M32   det="+det);
    prt_mat("E32",E32);
    mat_mul(E32, E31, T1);
    mat_mul(T1, E21, T2);
    mat_mul(T2, A, E32p);
    det = Matrix.determinant(E32p);
    System.out.println("E32p=E32*E31*E21*A     det="+det);
    prt_mat("E32p",E32p);

    mat_add(eye, M23, E23);
    det = Matrix.determinant(E23);
    System.out.println("E23=eye+M23   det="+det);
    prt_mat("E23",E23);
    mat_add(eye, M13, E13);
    det = Matrix.determinant(E13);
    System.out.println("E13=eye+M13   det="+det);
    prt_mat("E13",E13);
    mat_mul(E13, E23, T1);
    mat_mul(T1, E32, T2);
    mat_mul(T2, E31, T3);
    mat_mul(T3, E21, T4);
    mat_mul(T4, A, E13p);
    det = Matrix.determinant(E13p);
    System.out.println("E13p=E13*E23*E32*E31*E21*A     det="+det);
    prt_mat("E13p",E13p);

    mat_add(eye, M12, E12);
    det = Matrix.determinant(E12);
    System.out.println("E12=eye+M12   det="+det);
    prt_mat("E12",E12);
    mat_mul(E12, E13, T1);
    mat_mul(T1, E23, T2);
    mat_mul(T2, E32, T3);
    mat_mul(T3, E31, T4);
    mat_mul(T4, E21, T5);
    mat_mul(T5, A, E12p);
    det = Matrix.determinant(E12p);
    System.out.println("E12p=E12*E13*E23*E32*E31*E21*A     det="+det);
    prt_mat("E12p",E12p);

    System.out.println("A is product of sparse matrix, any can be combined into 2 matrix");
    System.out.println("A=E21i*E31i*E32i*E23i*E13i*E12i*D");
    new inverse(E21, E21i);
    prt_mat("E21i",E21i);
    new inverse(E31, E31i);
    prt_mat("E31i",E31i);
    new inverse(E32, E32i);
    prt_mat("E32i",E32i);
    new inverse(E23, E23i);
    prt_mat("E23i",E23i);
    new inverse(E13, E13i);
    prt_mat("E13i",E13i);
    new inverse(E12, E12i);
    prt_mat("E12i",E12i);
    prt_mat("D",D);
    mat_mul(E21i, E31i, T1);
    mat_mul(T1, E32i, T2);
    mat_mul(T2, E23i, T3);
    mat_mul(T3, E13i, T4);
    mat_mul(T4, E12i, T5);
    mat_mul(T5, D, Ack);
    det = Matrix.determinant(Ack);
    System.out.println("Ack = A   det="+det);
    prt_mat("Ack",Ack);
    new inverse(Ack, AI);
    det = Matrix.determinant(AI);
    System.out.println("AI = inverse(A)   det="+det);
    prt_mat("AI",AI);

    System.out.println("decomp2.java finished");
  }

  void mat_mul(double M1[][], double M2[][], double MP[][])
  {
    int n = MP.length;
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        MP[i][j] = 0.0;
	for(int k=0; k<n; k++)
	{
	  MP[i][j] = MP[i][j] + M1[i][k]*M2[k][j];
	}
      }
    }		      
  }

  public void mat_add(double B[][], double C[][], double A[][])
  {
    int n = A.length;	
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        A[i][j] = B[i][j] + C[i][j];
      }
    }
  } // end mat_add

  public void prt_mat(String name, double A[][])
  {
    int n = A.length;	
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        System.out.println(name+"["+i+"]["+j+"]="+A[i][j]);
      }
    }
    System.out.println(" ");
  } // end prt_mat

  void ident(double A[][])
  {
    int n = A.length;	
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
	if(i==j)
	{
	  A[i][j] = 1.0;
        }
      else
      {
	  A[i][j] = 0.0;
      }
      }
    }      
  } // end ident
    
  public static void main (String[] args)
  {
    new decomp2();
  }
} // end decomp2.java