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(:))   * @exception  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");    }  }    /**   * Returns true if the matrix is symmetric.   *   * @return boolean true if matrix is symmetric.   * @author FracPete, taken from old weka.core.Matrix class   */  public boolean isSymmetric() {    int nr = A.length, nc = A[0].length;    if (nr != nc)      return false;    for (int i = 0; i < nc; i++) {      for (int j = 0; j < i; j++) {        if (A[i][j] != A[j][i])          return false;      }    }    return true;  }  /**   * returns whether the matrix is a square matrix or not.   *   * @return true if the matrix is a square matrix   * @author FracPete   */  public boolean isSquare() {    return (getRowDimension() == getColumnDimension());  }  /**    * 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   * @exception  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));  }  /**   * returns the square root of the matrix, i.e., X from the equation   * X*X = A.<br/>   * Steps in the Calculation (see <a href="http://www.mathworks.com/access/helpdesk/help/techdoc/ref/sqrtm.html" target="blank"><code>sqrtm</code></a> in Matlab):<br/>   * <ol>   *   <li>perform eigenvalue decomposition<br/>[V,D]=eig(A)</li>   *   <li>take the square root of all elements in D (only the ones with    *       positive sign are considered for further computation)<br/>   *       S=sqrt(D)</li>   *   <li>calculate the root<br/>   *       X=V*S/V, which can be also written as X=(V'\(V*S)')'</li>   * </ol>   * <p/>   * <b>Note:</b> since this method uses other high-level methods, it generates   * several instances of matrices. This can be problematic with large   * matrices.   * <p/>   * Examples:   * <ol>   *   <li>   *   <pre>   *  X =   *   5   -4    1    0    0   *  -4    6   -4    1    0   *   1   -4    6   -4    1   *   0    1   -4    6   -4   *   0    0    1   -4    5   *    *  sqrt(X) =   *   2   -1   -0   -0   -0    *  -1    2   -1    0   -0    *   0   -1    2   -1    0    *  -0    0   -1    2   -1    *  -0   -0   -0   -1    2    *     *  Matrix m = new Matrix(new double[][]{{5,-4,1,0,0},{-4,6,-4,1,0},{1,-4,6,-4,1},{0,1,-4,6,-4},{0,0,1,-4,5}});   *   </pre>   *   </li>   *   <li>   *   <pre>   *  X =   *   7   10   *  15   22   *     *  sqrt(X) =   *  1.5667    1.7408   *  2.6112    4.1779   *    *  Matrix m = new Matrix(new double[][]{{7, 10},{15, 22}});   *   </pre>   *   </li>   * </ol>   *   * @return    sqrt(A)   * @author    FracPete   */  public Matrix sqrt() {    EigenvalueDecomposition   evd;    Matrix                    s;    Matrix                    v;    Matrix                    d;    Matrix                    result;    Matrix                    a;    Matrix                    b;    int                       i;    int                       n;    result = null;        // eigenvalue decomp.    // [V, D] = eig(A) with A = this    evd = this.eig();    v   = evd.getV();    d   = evd.getD();    // S = sqrt of cells of D    s = new Matrix(d.getRowDimension(), d.getColumnDimension());    for (i = 0; i < s.getRowDimension(); i++)