///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Faddeev.java                                       //
//                                                                       //
// © Tao Pang 2006                                                       //
//                                                                       //
// Last modified: August 1, 2007                                         //
//                                                                       //
// (1) This Java program is part of the book, "An Introduction to        //
//     Computational Physics, 2nd Edition," written by Tao Pang and      //
//     published by Cambridge University Press on January 19, 2006.      //
//                                                                       //
// (2) No warranties, express or implied, are made for this program.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

// An example of using the Faddeev-Leverrier method to
// obtain the inversion of a matrix.

import java.lang.*;
public class Faddeev {
  public static void main(String argv[]) {
    double a[][]= {{ 1,  3,  2},
                   {-2,  3, -1},
                   {-3, -1,  2}};
    int n = a.length;
    double c[] = new double[n];
    double d[][] = new double[n][n];
    double s[][][] = (double[][][]) fl(a, c);
    for (int i=0; i<n; ++i)
      for (int j=0; j<n; ++j)
        d[i][j] = -s[n-1][i][j]/c[0];
    for (int i=0; i<n; ++i)
      for (int j=0; j<n; ++j)
        System.out.println(d[i][j]);
  }

// Method to complete the Faddeev-Leverrier recursion.
 
  public static double[][][] fl(double a[][],
    double c[]) {
    int n = c.length;
    double s[][][] = new double[n][n][n]; 
    for (int i=0; i<n; ++i) s[0][i][i] = 1;
    for (int k=1; k<n; ++k) {
      s[k] = mm(a,s[k-1]);
      c[n-k] = -tr(mm(a,s[k-1]))/k;
      for (int i=0; i<n; ++i)
      s[k][i][i] += c[n-k];
    }
      c[0] = -tr(mm(a,s[n-1]))/n;
    return s;
  }
       
// Method to calculate the trace of a matrix.
 
  public static double tr(double a[][]) { 
    int n = a.length;
    double sum = 0;
    for (int i=0; i<n; ++i) sum += a[i][i];
    return sum;
  }

// Method to evaluate the product of two matrices.

  public static double[][] mm(double a[][], double b[][]) { 
    int n = a.length;
    double c[][] = new double[n][n];
    for (int i=0; i<n; ++i)
      for (int j=0; j<n; ++j) 
        for (int k=0; k<n; ++k) 
          c[i][j] += a[i][k]*b[k][j];
    return c;
  }
}