content: package Jama;
  
  import java.io.BufferedReader;
  import java.io.PrintWriter;
  import java.io.StreamTokenizer;
  import java.text.DecimalFormat;
  import java.text.DecimalFormatSymbols;
  import java.text.NumberFormat;
  import java.util.Locale;
  
  import Jama.util.Maths;
  
  /**
   * Jama = Java Matrix class.
   * <p/>
   * The Java Matrix Class provides the fundamental operations of numerical
   * linear algebra.  Various constructors create Matrices from two dimensional
   * arrays of double precision floating point numbers.  Various "gets" and
   * "sets" provide access to submatrices and matrix elements.  Several methods
   * implement basic matrix arithmetic, including matrix addition and
   * multiplication, matrix norms, and element-by-element array operations.
   * Methods for reading and printing matrices are also included.  All the
   * operations in this version of the Matrix Class involve real matrices.
   * Complex matrices may be handled in a future version.
   * <p/>
   * Five fundamental matrix decompositions, which consist of pairs or triples
   * of matrices, permutation vectors, and the like, produce results in five
   * decomposition classes.  These decompositions are accessed by the Matrix
   * class to compute solutions of simultaneous linear equations, determinants,
   * inverses and other matrix functions.  The five decompositions are:
   * <P><UL>
   * <LI>Cholesky Decomposition of symmetric, positive definite matrices.
   * <LI>LU Decomposition of rectangular matrices.
   * <LI>QR Decomposition of rectangular matrices.
   * <LI>Singular Value Decomposition of rectangular matrices.
   * <LI>Eigenvalue Decomposition of both symmetric and nonsymmetric square matrices.
   * </UL>
   * <DL>
   * <DT><B>Example of use:</B></DT>
   * <p/>
   * <DD>Solve a linear system A x = b and compute the residual norm, ||b - A x||.
   * <P><PRE>
   * double[][] vals = {{1.,2.,3},{4.,5.,6.},{7.,8.,10.}};
   * Matrix A = new Matrix(vals);
   * Matrix b = Matrix.random(3,1);
   * Matrix x = A.solve(b);
   * Matrix r = A.times(x).minus(b);
   * double rnorm = r.normInf();
   * </PRE></DD>
   * </DL>
   *
   * @author The MathWorks, Inc. and the National Institute of Standards and Technology.
   * @version 5 August 1998
   */
  
  public class Matrix implements Cloneable, java.io.Serializable {
  
  /* ------------------------
     Class variables
   * ------------------------ */
  
      /**
       * Array for internal storage of elements.
       *
       * @serial internal array storage.
       */
      private double[][] A;
  
      /**
       * Row and column dimensions.
       *
       * @serial row dimension.
       * @serial column dimension.
       */
      private int m, n;
  
  /* ------------------------
     Constructors
   * ------------------------ */
  
      /**
       * Construct an m-by-n matrix of zeros.
       *
       * @param m Number of rows.
       * @param n Number of colums.
       */
  
      public Matrix(int m, int n) {
          this.m = m;
          this.n = n;
          A = new double[m][n];
      }
  
      /**
       * Construct an m-by-n constant matrix.
       *
       * @param m Number of rows.
       * @param n Number of colums.
       * @param s Fill the matrix with this scalar value.
       */
  
      public Matrix(int m, int n, double s) {
          this.m = m;
          this.n = n;
          A = new double[m][n];
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  A[i][j] = s;
              }
          }
      }
  
      /**
       * Construct a matrix from a 2-D array.
       *
       * @param A Two-dimensional array of doubles.
       * @throws IllegalArgumentException All rows must have the same length
       * @see #constructWithCopy
       */
  
      public Matrix(double[][] A) {
          m = A.length;
          n = A[0].length;
          for (int i = 0; i < m; i++) {
              if (A[i].length != n) {
                  throw new IllegalArgumentException("All rows must have the same length.");
              }
          }
          this.A = A;
      }
  
      /**
       * Construct a matrix quickly without checking arguments.
       *
       * @param A Two-dimensional array of doubles.
       * @param m Number of rows.
       * @param n Number of colums.
       */
  
      public Matrix(double[][] A, int m, int n) {
          this.A = A;
          this.m = m;
          this.n = n;
      }
  
      /**
       * Construct a matrix from a one-dimensional packed array
       *
       * @param vals One-dimensional array of doubles, packed by columns (ala Fortran).
       * @param m    Number of rows.
       * @throws IllegalArgumentException Array length must be a multiple of m.
       */
  
      public Matrix(double vals[], int m) {
          this.m = m;
          n = (m != 0 ? vals.length / m : 0);
          if (m * n != vals.length) {
              throw new IllegalArgumentException("Array length must be a multiple of m.");
          }
          A = new double[m][n];
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  A[i][j] = vals[i + j * m];
              }
          }
      }
  
  /* ------------------------
     Public Methods
   * ------------------------ */
  
      /**
       * Construct a matrix from a copy of a 2-D array.
       *
       * @param A Two-dimensional array of doubles.
       * @throws IllegalArgumentException All rows must have the same length
       */
  
      public static Matrix constructWithCopy(double[][] A) {
          int m = A.length;
          int n = A[0].length;
          Matrix X = new Matrix(m, n);
          double[][] C = X.getArray();
          for (int i = 0; i < m; i++) {
              if (A[i].length != n) {
                  throw new IllegalArgumentException
                          ("All rows must have the same length.");
              }
              for (int j = 0; j < n; j++) {
                  C[i][j] = A[i][j];
              }
          }
          return X;
      }
  
      /**
       * Generate matrix with random elements
       *
       * @param m Number of rows.
       * @param n Number of colums.
       * @return An m-by-n matrix with uniformly distributed random elements.
       */
  
      public static Matrix random(int m, int n) {
          Matrix A = new Matrix(m, n);
          double[][] X = A.getArray();
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  X[i][j] = Math.random();
              }
          }
          return A;
      }
  
      /**
       * Generate identity matrix
       *
       * @param m Number of rows.
       * @param n Number of colums.
       * @return An m-by-n matrix with ones on the diagonal and zeros elsewhere.
       */
  
      public static Matrix identity(int m, int n) {
          Matrix A = new Matrix(m, n);
          double[][] X = A.getArray();
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  X[i][j] = (i == j ? 1.0 : 0.0);
              }
          }
          return A;
      }
  
      /**
       * Read a matrix from a stream.  The format is the same the print method,
       * so printed matrices can be read back in (provided they were printed using
       * US Locale).  Elements are separated by
       * whitespace, all the elements for each row appear on a single line,
       * the last row is followed by a blank line.
       *
       * @param input the input stream.
       */
  
      public static Matrix read(BufferedReader input) throws java.io.IOException {
          StreamTokenizer tokenizer = new StreamTokenizer(input);
  
          // Although StreamTokenizer will parse numbers, it doesn't recognize
          // scientific notation (E or D); however, Double.valueOf does.
          // The strategy here is to disable StreamTokenizer's number parsing.
          // We'll only get whitespace delimited words, EOL's and EOF's.
          // These words should all be numbers, for Double.valueOf to parse.
  
          tokenizer.resetSyntax();
          tokenizer.wordChars(0, 255);
          tokenizer.whitespaceChars(0, ' ');
          tokenizer.eolIsSignificant(true);
          java.util.Vector v = new java.util.Vector();
  
          // Ignore initial empty lines
          while (tokenizer.nextToken() == StreamTokenizer.TT_EOL) ;
          if (tokenizer.ttype == StreamTokenizer.TT_EOF)
              throw new java.io.IOException("Unexpected EOF on matrix read.");
          do {
              v.addElement(Double.valueOf(tokenizer.sval)); // Read & store 1st row.
          } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
  
          int n = v.size();  // Now we've got the number of columns!
          double row[] = new double[n];
          for (int j = 0; j < n; j++)  // extract the elements of the 1st row.
              row[j] = ((Double) v.elementAt(j)).doubleValue();
          v.removeAllElements();
          v.addElement(row);  // Start storing rows instead of columns.
          while (tokenizer.nextToken() == StreamTokenizer.TT_WORD) {
              // While non-empty lines
              v.addElement(row = new double[n]);
              int j = 0;
              do {
                  if (j >= n) throw new java.io.IOException
                          ("Row " + v.size() + " is too long.");
                  row[j++] = Double.valueOf(tokenizer.sval).doubleValue();
              } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
              if (j < n) throw new java.io.IOException
                      ("Row " + v.size() + " is too short.");
          }
          int m = v.size();  // Now we've got the number of rows.
          double[][] A = new double[m][];
          v.copyInto(A);  // copy the rows out of the vector
          return new Matrix(A);
      }
  
      /**
       * Make a deep copy of a matrix
       */
  
      public Matrix copy() {
          Matrix X = new Matrix(m, n);
          double[][] C = X.getArray();
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  C[i][j] = A[i][j];
              }
          }
          return X;
      }
  
      /**
       * Clone the Matrix object.
       */
  
      @Override
      public Object clone() {
          return this.copy();
      }
  
      /**
       * Access the internal two-dimensional array.
       *
       * @return Pointer to the two-dimensional array of matrix elements.
       */
  
      public double[][] getArray() {
          return A;
      }
  
      /**
       * Copy the internal two-dimensional array.
       *
       * @return Two-dimensional array copy of matrix elements.
       */
  
      public double[][] getArrayCopy() {
          double[][] C = new double[m][n];
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  C[i][j] = A[i][j];
              }
          }
          return C;
      }
  
      /**
       * Make a one-dimensional column packed copy of the internal array.
       *
       * @return Matrix elements packed in a one-dimensional array by columns.
       */
  
      public double[] getColumnPackedCopy() {
          double[] vals = new double[m * n];
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  vals[i + j * m] = A[i][j];
              }
          }
          return vals;
      }
  
      /**
       * Make a one-dimensional row packed copy of the internal array.
       *
       * @return Matrix elements packed in a one-dimensional array by rows.
       */
  
      public double[] getRowPackedCopy() {
          double[] vals = new double[m * n];
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  vals[i * n + j] = A[i][j];
              }
          }
          return vals;
      }
  
      /**
       * Get row dimension.
       *
       * @return m, the number of rows.
       */
  
      public int getRowDimension() {
          return m;
      }
  
      /**
       * Get column dimension.
       *
       * @return n, the number of columns.
       */
  
      public int getColumnDimension() {
          return n;
      }
  
      /**
       * Get a single element.
       *
       * @param i Row index.
       * @param j Column index.
       * @return A(i, j)
       * @throws ArrayIndexOutOfBoundsException
       */
  
      public double get(int i, int j) {
          return A[i][j];
      }
  
      /**
       * Get a submatrix.
       *
       * @param i0 Initial row index
       * @param i1 Final row index
       * @param j0 Initial column index
       * @param j1 Final column index
       * @return A(i0:i1, j0:j1)
       * @throws ArrayIndexOutOfBoundsException Submatrix indices
       */
  
      public Matrix getMatrix(int i0, int i1, int j0, int j1) {
          Matrix X = new Matrix(i1 - i0 + 1, j1 - j0 + 1);
          double[][] B = X.getArray();
          try {
              for (int i = i0; i <= i1; i++) {
                  for (int j = j0; j <= j1; j++) {
                      B[i - i0][j - j0] = A[i][j];
                  }
              }
          } catch (ArrayIndexOutOfBoundsException e) {
              throw new ArrayIndexOutOfBoundsException("Submatrix indices");
          }
          return X;
      }
  
      /**
       * Get a submatrix.
       *
       * @param r Array of row indices.
       * @param c Array of column indices.
       * @return A(r(:), c(:))
       * @throws ArrayIndexOutOfBoundsException Submatrix indices
       */
  
      public Matrix getMatrix(int[] r, int[] c) {
          Matrix X = new Matrix(r.length, c.length);
          double[][] B = X.getArray();
          try {
              for (int i = 0; i < r.length; i++) {
                  for (int j = 0; j < c.length; j++) {
                      B[i][j] = A[r[i]][c[j]];
                  }
              }
          } catch (ArrayIndexOutOfBoundsException e) {
              throw new ArrayIndexOutOfBoundsException("Submatrix indices");
          }
          return X;
      }
  
      /**
       * Get a submatrix.
       *
       * @param i0 Initial row index
       * @param i1 Final row index
       * @param c  Array of column indices.
       * @return A(i0:i1, c(:))
       * @throws ArrayIndexOutOfBoundsException Submatrix indices
       */
  
      public Matrix getMatrix(int i0, int i1, int[] c) {
          Matrix X = new Matrix(i1 - i0 + 1, c.length);
          double[][] B = X.getArray();
          try {
              for (int i = i0; i <= i1; i++) {
                  for (int j = 0; j < c.length; j++) {
                      B[i - i0][j] = A[i][c[j]];
                  }
              }
          } catch (ArrayIndexOutOfBoundsException e) {
              throw new ArrayIndexOutOfBoundsException("Submatrix indices");
          }
          return X;
      }
  
      /**
       * Get a submatrix.
       *
       * @param r  Array of row indices.
       * @param i0 Initial column index
       * @param i1 Final column index
       * @return A(r(:), j0:j1)
       * @throws ArrayIndexOutOfBoundsException Submatrix indices
       */
  
      public Matrix getMatrix(int[] r, int j0, int j1) {
          Matrix X = new Matrix(r.length, j1 - j0 + 1);
          double[][] B = X.getArray();
          try {
              for (int i = 0; i < r.length; i++) {
                  for (int j = j0; j <= j1; j++) {
                      B[i][j - j0] = A[r[i]][j];
                  }
              }
          } catch (ArrayIndexOutOfBoundsException e) {
              throw new ArrayIndexOutOfBoundsException("Submatrix indices");
          }
          return X;
      }
  
      /**
       * Set a single element.
       *
       * @param i Row index.
       * @param j Column index.
       * @param s A(i,j).
       * @throws ArrayIndexOutOfBoundsException
       */
  
      public void set(int i, int j, double s) {
          A[i][j] = s;
      }
  
      /**
       * Set a submatrix.
       *
       * @param i0 Initial row index
       * @param i1 Final row index
       * @param j0 Initial column index
       * @param j1 Final column index
       * @param X  A(i0:i1,j0:j1)
       * @throws ArrayIndexOutOfBoundsException Submatrix indices
       */
  
      public void setMatrix(int i0, int i1, int j0, int j1, Matrix X) {
          try {
              for (int i = i0; i <= i1; i++) {
                  for (int j = j0; j <= j1; j++) {
                      A[i][j] = X.get(i - i0, j - j0);
                  }
              }
          } catch (ArrayIndexOutOfBoundsException e) {
              throw new ArrayIndexOutOfBoundsException("Submatrix indices");
          }
      }
  
      /**
       * Set a submatrix.
       *
       * @param r Array of row indices.
       * @param c Array of column indices.
       * @param X A(r(:),c(:))
       * @throws ArrayIndexOutOfBoundsException Submatrix indices
       */
  
      public void setMatrix(int[] r, int[] c, Matrix X) {
          try {
              for (int i = 0; i < r.length; i++) {
                  for (int j = 0; j < c.length; j++) {
                      A[r[i]][c[j]] = X.get(i, j);
                  }
              }
          } catch (ArrayIndexOutOfBoundsException e) {
              throw new ArrayIndexOutOfBoundsException("Submatrix indices");
          }
      }
  
      /**
       * Set a submatrix.
       *
       * @param r  Array of row indices.
       * @param j0 Initial column index
       * @param j1 Final column index
       * @param X  A(r(:),j0:j1)
       * @throws ArrayIndexOutOfBoundsException Submatrix indices
       */
  
      public void setMatrix(int[] r, int j0, int j1, Matrix X) {
          try {
              for (int i = 0; i < r.length; i++) {
                  for (int j = j0; j <= j1; j++) {
                      A[r[i]][j] = X.get(i, j - j0);
                  }
              }
          } catch (ArrayIndexOutOfBoundsException e) {
              throw new ArrayIndexOutOfBoundsException("Submatrix indices");
          }
      }
  
      /**
       * Set a submatrix.
       *
       * @param i0 Initial row index
       * @param i1 Final row index
       * @param c  Array of column indices.
       * @param X  A(i0:i1,c(:))
       * @throws ArrayIndexOutOfBoundsException Submatrix indices
       */
  
      public void setMatrix(int i0, int i1, int[] c, Matrix X) {
          try {
              for (int i = i0; i <= i1; i++) {
                  for (int j = 0; j < c.length; j++) {
                      A[i][c[j]] = X.get(i - i0, j);
                  }
              }
          } catch (ArrayIndexOutOfBoundsException e) {
              throw new ArrayIndexOutOfBoundsException("Submatrix indices");
          }
      }
  
      /**
       * Matrix transpose.
       *
       * @return A'
       */
  
      public Matrix transpose() {
          Matrix X = new Matrix(n, m);
          double[][] C = X.getArray();
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  C[j][i] = A[i][j];
              }
          }
          return X;
      }
  
      /**
       * One norm
       *
       * @return maximum column sum.
       */
  
      public double norm1() {
          double f = 0;
          for (int j = 0; j < n; j++) {
              double s = 0;
              for (int i = 0; i < m; i++) {
                  s += Math.abs(A[i][j]);
              }
              f = Math.max(f, s);
          }
          return f;
      }
  
      /**
       * Two norm
       *
       * @return maximum singular value.
       */
  
      public double norm2() {
          return (new SingularValueDecomposition(this).norm2());
      }
  
      /**
       * Infinity norm
       *
       * @return maximum row sum.
       */
  
      public double normInf() {
          double f = 0;
          for (int i = 0; i < m; i++) {
              double s = 0;
              for (int j = 0; j < n; j++) {
                  s += Math.abs(A[i][j]);
              }
              f = Math.max(f, s);
          }
          return f;
      }
  
      /**
       * Frobenius norm
       *
       * @return sqrt of sum of squares of all elements.
       */
  
      public double normF() {
          double f = 0;
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  f = Maths.hypot(f, A[i][j]);
              }
          }
          return f;
      }
  
      /**
       * Unary minus
       *
       * @return -A
       */
  
      public Matrix uminus() {
          Matrix X = new Matrix(m, n);
          double[][] C = X.getArray();
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  C[i][j] = -A[i][j];
              }
          }
          return X;
      }
  
      /**
       * C = A + B
       *
       * @param B another matrix
       * @return A + B
       */
  
      public Matrix plus(Matrix B) {
          checkMatrixDimensions(B);
          Matrix X = new Matrix(m, n);
          double[][] C = X.getArray();
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  C[i][j] = A[i][j] + B.A[i][j];
              }
          }
          return X;
      }
  
      /**
       * A = A + B
       *
       * @param B another matrix
       * @return A + B
       */
  
      public Matrix plusEquals(Matrix B) {
          checkMatrixDimensions(B);
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  A[i][j] = A[i][j] + B.A[i][j];
              }
          }
          return this;
      }
  
      /**
       * C = A - B
       *
       * @param B another matrix
       * @return A - B
       */
  
      public Matrix minus(Matrix B) {
          checkMatrixDimensions(B);
          Matrix X = new Matrix(m, n);
          double[][] C = X.getArray();
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  C[i][j] = A[i][j] - B.A[i][j];
              }
          }
          return X;
      }
  
      /**
       * A = A - B
       *
       * @param B another matrix
       * @return A - B
       */
  
      public Matrix minusEquals(Matrix B) {
          checkMatrixDimensions(B);
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  A[i][j] = A[i][j] - B.A[i][j];
              }
          }
          return this;
      }
  
      /**
       * Element-by-element multiplication, C = A.*B
       *
       * @param B another matrix
       * @return A.*B
       */
  
      public Matrix arrayTimes(Matrix B) {
          checkMatrixDimensions(B);
          Matrix X = new Matrix(m, n);
          double[][] C = X.getArray();
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  C[i][j] = A[i][j] * B.A[i][j];
              }
          }
          return X;
      }
  
      /**
       * Element-by-element multiplication in place, A = A.*B
       *
       * @param B another matrix
       * @return A.*B
       */
  
      public Matrix arrayTimesEquals(Matrix B) {
          checkMatrixDimensions(B);
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  A[i][j] = A[i][j] * B.A[i][j];
              }
          }
          return this;
      }
  
      /**
       * Element-by-element right division, C = A./B
       *
       * @param B another matrix
       * @return A./B
       */
  
      public Matrix arrayRightDivide(Matrix B) {
          checkMatrixDimensions(B);
          Matrix X = new Matrix(m, n);
          double[][] C = X.getArray();
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  C[i][j] = A[i][j] / B.A[i][j];
              }
          }
          return X;
      }
  
      /**
       * Element-by-element right division in place, A = A./B
       *
       * @param B another matrix
       * @return A./B
       */
  
      public Matrix arrayRightDivideEquals(Matrix B) {
          checkMatrixDimensions(B);
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  A[i][j] = A[i][j] / B.A[i][j];
              }
          }
          return this;
      }
  
      /**
       * Element-by-element left division, C = A.\B
       *
       * @param B another matrix
       * @return A.\B
       */
  
      public Matrix arrayLeftDivide(Matrix B) {
          checkMatrixDimensions(B);
          Matrix X = new Matrix(m, n);
          double[][] C = X.getArray();
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  C[i][j] = B.A[i][j] / A[i][j];
              }
          }
          return X;
      }
  
      /**
       * Element-by-element left division in place, A = A.\B
       *
       * @param B another matrix
       * @return A.\B
       */
  
      public Matrix arrayLeftDivideEquals(Matrix B) {
          checkMatrixDimensions(B);
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  A[i][j] = B.A[i][j] / A[i][j];
              }
          }
          return this;
      }
  
      /**
       * Multiply a matrix by a scalar, C = s*A
       *
       * @param s scalar
       * @return s*A
       */
  
      public Matrix times(double s) {
          Matrix X = new Matrix(m, n);
          double[][] C = X.getArray();
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  C[i][j] = s * A[i][j];
              }
          }
          return X;
      }
  
      /**
       * Multiply a matrix by a scalar in place, A = s*A
       *
       * @param s scalar
       * @return replace A by s*A
       */
  
      public Matrix timesEquals(double s) {
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  A[i][j] = s * A[i][j];
              }
          }
          return this;
      }
  
      /**
       * Linear algebraic matrix multiplication, A * B
       *
       * @param B another matrix
       * @return Matrix product, A * B
       * @throws IllegalArgumentException Matrix inner dimensions must agree.
       */
  
      public Matrix times(Matrix B) {
          if (B.m != n) {
              throw new IllegalArgumentException("Matrix inner dimensions must agree.");
          }
          Matrix X = new Matrix(m, B.n);
          double[][] C = X.getArray();
          double[] Bcolj = new double[n];
          for (int j = 0; j < B.n; j++) {
              for (int k = 0; k < n; k++) {
                  Bcolj[k] = B.A[k][j];
              }
              for (int i = 0; i < m; i++) {
                  double[] Arowi = A[i];
                  double s = 0;
                  for (int k = 0; k < n; k++) {
                      s += Arowi[k] * Bcolj[k];
                  }
                  C[i][j] = s;
              }
          }
          return X;
      }
  
      /**
       * LU Decomposition
       *
       * @return LUDecomposition
       * @see LUDecomposition
       */
  
      public LUDecomposition lu() {
          return new LUDecomposition(this);
      }
  
      /**
       * QR Decomposition
       *
       * @return QRDecomposition
       * @see QRDecomposition
       */
  
      public QRDecomposition qr() {
          return new QRDecomposition(this);
      }
  
      /**
       * Cholesky Decomposition
       *
       * @return CholeskyDecomposition
       * @see CholeskyDecomposition
       */
  
      public CholeskyDecomposition chol() {
          return new CholeskyDecomposition(this);
      }
  
      /**
       * Singular Value Decomposition
       *
       * @return SingularValueDecomposition
       * @see SingularValueDecomposition
       */
  
      public SingularValueDecomposition svd() {
          return new SingularValueDecomposition(this);
      }
  
      /**
       * Eigenvalue Decomposition
       *
       * @return EigenvalueDecomposition
       * @see EigenvalueDecomposition
       */
  
      public EigenvalueDecomposition eig() {
          return new EigenvalueDecomposition(this);
      }
  
      /**
       * Solve A*X = B
       *
       * @param B right hand side
       * @return solution if A is square, least squares solution otherwise
       */
  
      public Matrix solve(Matrix B) {
          return (m == n ? (new LUDecomposition(this)).solve(B) :
                  (new QRDecomposition(this)).solve(B));
      }
  
      /**
       * Solve X*A = B, which is also A'*X' = B'
       *
       * @param B right hand side
       * @return solution if A is square, least squares solution otherwise.
       */
  
      public Matrix solveTranspose(Matrix B) {
          return transpose().solve(B.transpose());
      }
  
      /**
       * Matrix inverse or pseudoinverse
       *
       * @return inverse(A) if A is square, pseudoinverse otherwise.
       */
  
      public Matrix inverse() {
          return solve(identity(m, m));
      }
  
      /**
       * Matrix determinant
       *
       * @return determinant
       */
  
      public double det() {
          return new LUDecomposition(this).det();
      }
  
      /**
       * Matrix rank
       *
       * @return effective numerical rank, obtained from SVD.
       */
  
      public int rank() {
          return new SingularValueDecomposition(this).rank();
      }
  
      /**
       * Matrix condition (2 norm)
       *
       * @return ratio of largest to smallest singular value.
       */
  
      public double cond() {
          return new SingularValueDecomposition(this).cond();
      }
  
      /**
       * Matrix trace.
       *
       * @return sum of the diagonal elements.
       */
  
      public double trace() {
          double t = 0;
          for (int i = 0; i < Math.min(m, n); i++) {
              t += A[i][i];
          }
          return t;
      }
  
      /**
       * Print the matrix to stdout.   Line the elements up in columns
       * with a Fortran-like 'Fw.d' style format.
       *
       * @param w Column width.
       * @param d Number of digits after the decimal.
       */
  
      public void print(int w, int d) {
          print(new PrintWriter(System.out, true), w, d);
      }
  
      /**
       * Print the matrix to the output stream.   Line the elements up in
       * columns with a Fortran-like 'Fw.d' style format.
       *
       * @param output Output stream.
       * @param w      Column width.
       * @param d      Number of digits after the decimal.
       */
  
      public void print(PrintWriter output, int w, int d) {
          DecimalFormat format = new DecimalFormat();
          format.setDecimalFormatSymbols(new DecimalFormatSymbols(Locale.US));
          format.setMinimumIntegerDigits(1);
          format.setMaximumFractionDigits(d);
          format.setMinimumFractionDigits(d);
          format.setGroupingUsed(false);
          print(output, format, w + 2);
      }
  
      // DecimalFormat is a little disappointing coming from Fortran or C's printf.
      // Since it doesn't pad on the left, the elements will come out different
      // widths.  Consequently, we'll pass the desired column width in as an
      // argument and do the extra padding ourselves.
  
      /**
       * Print the matrix to stdout.  Line the elements up in columns.
       * Use the format object, and right justify within columns of width
       * characters.
       * Note that is the matrix is to be read back in, you probably will want
       * to use a NumberFormat that is set to US Locale.
       *
       * @param format A  Formatting object for individual elements.
       * @param width  Field width for each column.
       * @see java.text.DecimalFormat#setDecimalFormatSymbols
       */
  
      public void print(NumberFormat format, int width) {
          print(new PrintWriter(System.out, true), format, width);
      }
  
      /**
       * Print the matrix to the output stream.  Line the elements up in columns.
       * Use the format object, and right justify within columns of width
       * characters.
       * Note that is the matrix is to be read back in, you probably will want
       * to use a NumberFormat that is set to US Locale.
       *
       * @param output the output stream.
       * @param format A formatting object to format the matrix elements
       * @param width  Column width.
       * @see java.text.DecimalFormat#setDecimalFormatSymbols
       */
  
      public void print(PrintWriter output, NumberFormat format, int width) {
          output.println();  // start on new line.
          for (int i = 0; i < m; i++) {
              for (int j = 0; j < n; j++) {
                  String s = format.format(A[i][j]); // format the number
                  int padding = Math.max(1, width - s.length()); // At _least_ 1 space
                  for (int k = 0; k < padding; k++)
                      output.print(' ');
                  output.print(s);
              }
              output.println();
          }
          output.println();   // end with blank line.
      }
  
  
  /* ------------------------
     Private Methods
   * ------------------------ */
  
      /**
       * Check if size(A) == size(B)
       **/
  
      private void checkMatrixDimensions(Matrix B) {
          if (B.m != m || B.n != n) {
              throw new IllegalArgumentException("Matrix dimensions must agree.");
          }
      }
  
  }
  