/////////////////////////////////////////////////////////////////////////////////////////////
//
//   Basis.java
//   Determines bases for the image and kernal of a matrix
//
//   Save this file as Basis.java on any system with the javac compiler, available
//   from Sun Microsystems: http://java.sun.com/javase/downloads/index.jsp
//   If you're on a unix system also save the file "Makefile" in the same directory, then
//        Compile by typing "make" or "gmake" at the command line (leave out the quotes).
//        Run the executable by doing "Basis infile outfile", where "infile" is
//        a properly formated (described below) input file.
//   If you're on any other system, then
//        Compile by typing "javac Basis.java" at the command line (leave out quotes).
//        Run the executable by doing "java Basis infile outfile, where "infile" is
//        a properly formated (described below) input file.
//
//   Note: if "outfile" already exists in your working directory, then the above commands will
//         overwrite it.
//
//   Input file format.
//   The first line of an input file consists of two integers n and m, separated by a space,
//   giving the number of rows and columns, respectively.  The next n lines of the input file
//   specify the n rows of the matrix.  Each line is a space separated list of m real numbers.
//   Note that the number of spaces separating the numbers in each row is irrelevant.
//
//   Example 1  (Save the following lines as a text file called in1.)
//   3 5
//    2  0  4  2  9
//    3  4  0  0 -7
//    0  1  0  5  0
//
//   Example 2  (Save as text file in2.)
//   4 7
//    3.1  4.2  7.9  0.0  0.0  1.1  8.2
//   -6.7  1.2  3.4  5.6  7.8  0.0  0.0
//    9.8  7.6  5.4  3.2  2.1  2.3  4.5
//    6.7  8.8 -9.1  2.5  3.6  4.7 -5.5
//
//   Output file format.
//   The output will consist of three matrices.  The input matrix (nxm) will be reprinted,
//   followed by a matrix with n rows and r columns, where r is the rank of the input matrix.
//   This second matrix will consist of a subset of the columns of the input matrix which
//   form a basis for it's image.  The third matrix will have m rows and m-r columns, which
//   form a basis for the kernal of the input matrix.  The format of the matrix entries can
//   be adjusted by changing the vaules of the constants PRECISION and COLUMN_WIDTH below
//   (lines 60 and 61 of this file).  PRECISION specifies the number of decimal digits to
//   the right of the decimal point (default 3), and COLUMN_WIDTH gives the minimum number
//   of characters to be printed in each column (default 10).  If the number printed exceeds
//   this width, the field will expand to accomodate the extra characters.  There is a
//   further space included after each column.  After changing these constants, you must
//   recompile the program.
//
/////////////////////////////////////////////////////////////////////////////////////////////

import java.io.*;
import java.util.StringTokenizer;

class Basis{

   static final int PRECISION = 3;
   static final int COLUMN_WIDTH = 10;

   // swap()
   // swap row i with row k
   // pre: A[i][q]==A[k][q]==0 for 1<=q<j
   static void swap(double[][] A, int i, int k, int j){
      int m = A[0].length - 1;
      double temp;
      for(int q=j; q<=m; q++){
         temp = A[i][q];
         A[i][q] = A[k][q];
         A[k][q] = temp;
      }
   }

   // divide()
   // divide row i by A[i][j]
   // pre: A[i][j]!=0, A[i][q]==0 for 1<=q<j
   // post: A[i][j]==1;
   static void divide(double[][] A, int i, int j){
      int m = A[0].length - 1;
      for(int q=j+1; q<=m; q++) A[i][q] /= A[i][j];
      A[i][j] = 1;
   }

   // eliminate()
   // subtract an appropriate multiple of row i from every other row
   // pre: A[i][j]==1, A[i][q]==0 for 1<=q<j
   // post: A[p][j]==0 for p!=i
   static void eliminate(double[][] A, int i, int j){
      int n = A.length - 1;
      int m = A[0].length - 1;
      for(int p=1; p<=n; p++){
         if( p!=i && A[p][j]!=0 ){
            for(int q=j+1; q<=m; q++){
               A[p][q] -= A[p][j]*A[i][q];
            }
            A[p][j] = 0;
         }
      }
   }

   // printMatrix()
   // print the present state of Matrix A to file out
   static void printMatrix(PrintWriter out, double[][] A){
      int n = A.length - 1;
      int m = A[0].length - 1;
      StringBuffer format = new StringBuffer("%");
      format.append(COLUMN_WIDTH).append(".").append(PRECISION).append("f ");
      for(int i=1; i<=n; i++){
         for(int j=1; j<=m; j++) out.printf(format.toString(), A[i][j]);
         out.println();
      }
      out.println();
      out.println();
   }

   // main()
   // read input file, initialize matrix, perform Gauss-Jordan Elimination,
   // and write resulting matrices to output file
   public static void main(String[] args) throws IOException {
      int n, m, r, p, i, j, k, q, s;
      String line;
      StringTokenizer st;

      // check command line arguments, open input and output files
      if( args.length!=2 ){
         System.out.println("Usage: Basis infile outfile");
         System.exit(1);
      }
      BufferedReader in = new BufferedReader(new FileReader(args[0]));
      PrintWriter out = new PrintWriter(new FileWriter(args[1]));

      // read first line of input file
      line = in.readLine();
      st = new StringTokenizer(line);
      n = Integer.parseInt(st.nextToken());
      m = Integer.parseInt(st.nextToken());

      // declare A and B to be of size (n+1)x(m+1) and do not use index 0
      double[][] A = new double[n+1][m+1];
      double[][] B = new double[n+1][m+1];

      // read next n lines of input file and initialize arrays A and B
      for(i=1; i<=n; i++){
         line = in.readLine();
         st = new StringTokenizer(line);
         for(j=1; j<=m; j++){
            B[i][j] = A[i][j] = Double.parseDouble(st.nextToken());
         }
      }

      // close input file
      in.close();

      // perform Gauss-Jordan Elimination algorithm on matrix B, mark columns which
      // contain a leading 1, determine the rank and nullity, save matrix A for later use
      i = 1; j = 1;
      while( i<=n && j<=m ){
         B[0][j] = 0;  // mark col j as not containing a leading 1
         k = i; while( k<=n && B[k][j]==0 ) k++;
         if( k<=n ){
            if( k!=i ) swap(B, i, k, j);
            if( B[i][j]!=1 ) divide(B, i, j);
            eliminate(B, i, j);
            B[0][j] = 1; // mark col j as containing a leading 1
            i++;
         }
         j++;
      }
      r = i-1;  // rank(A)
      p = m-r;  // nullity(A)

      // determine an mxp matrix C whose columns form a basis of ker(A)
      double[][] C = new double[m+1][p+1];
      j = 1; k = 1;
      for(i=1; i<=m; i++){
         if( B[0][i]==1 ){
            for(q=1; q<j; q++) C[i][q] = 0;
            s = i;
            for(q=j; q<=p; q++){
               s++;
               while(B[0][s]==1) s++;
               C[i][q] = -B[k][s];
            }
            k++;
         }else{
            for(q=1; q<j; q++) C[i][q] = 0;
            C[i][j] = 1;
            for(q=j+1; q<=p; q++) C[i][q] = 0;
            j++;
         }
      }

      // determine a nxr matrix D whose columns form a basis of im(A)
      double[][] D = new double[n+1][r+1];
      k = 1;
      for(j=1; j<=m; j++){
         if( B[0][j]==1 ){
            for(i=1; i<=n; i++) D[i][k] = A[i][j];
            k++;
         }
      }

      // print A, C and D
      out.println("Input matrix A:");
      printMatrix(out, A);
      out.println("Basis for im(A):");
      printMatrix(out, D);
      out.println("Basis for ker(A):");
      printMatrix(out, C);


      // close output file
      out.close();
   }
}