// decomp3a.java  given A create B1 and B2 ... such that B1 * B2 ... = A
// from mathoverflow.net
// uses inverse.java 
// uses Matrix.java for utility routines. zero based subscripts by row

public class decomp3a
{
  int n = 3;
  double det;
  double A[][] = {{7.0, 5.0, 3.0},  //   0,0  0,1  0,2
       	          {2.0, 9.0, 4.0},  //   1,0  1,1  1,2
		  {6.0, 1.0, 8.0}}; //   2,0  2,1  2,2
  double E21[][] = new double[n][n];   // eye +  -A[1][0]/A[0][0]
  double E31[][] = new double[n][n];   // eye +  -A[2][0]/A[0][0]
  double E31p[][] = new double[n][n];  // E31 * E21 * A  product
  double E32[][] = new double[n][n];   // eye + E31p[2][1]/-E31p[1][1]
  double E32p[][] = new double[n][n];  // E32 * E31 * E21 * A
  double E23[][] = new double[n][n];   // eye + 
  double E13[][] = new double[n][n];   // eye + 
  double E13p[][] = new double[n][n];  // E13 * E23 * E32 * E31 * E21 * A
  double E12[][] = new double[n][n];   // eye + 
  double E12p[][] = new double[n][n];  // E12*E13*E23*E32*E31*E21*A
  double D[][] = new double[n][n];     // just diagonal of E12p

  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 decomp3a()
  {
    System.out.println("decomp3a.java running");

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

    ident(E21);
    E21[1][0] = -A[1][0]/A[0][0];
    det = Matrix.determinant(E21);
    System.out.println("E21=         det="+det);
    prt_mat("E21",E21);

    ident(E31);
    E31[2][0] = -A[2][0]/A[0][0];
    System.out.println("E31=   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);

    ident(E32);
    E32[2][1] = E31p[2][1]/(-E31p[1][1]); 
    det = Matrix.determinant(E32);
    System.out.println("E32=     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);

    ident(E23);
    E23[1][2] = -E32p[1][2]/E32p[2][2]; 
    det = Matrix.determinant(E23);
    System.out.println("E23=     det="+det);
    prt_mat("E23",E23);
    ident(E13);
    E13[0][2] = -E32p[0][2]/E32p[2][2];
    det = Matrix.determinant(E13);
    System.out.println("E13=     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);

    ident(E12);
    E12[0][1] = E13p[0][1]/(-E13p[1][1]);
    det = Matrix.determinant(E12);
    System.out.println("E12=     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);

    ident(D);
    D[0][0] = E12p[0][0];
    D[1][1] = E12p[1][1];
    D[2][2] = E12p[2][2];
    prt_mat("D", D);
    
    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("decomp3a.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 decomp3a();
  }
} // end decomp3a.java