Creates an LU decomposition of the given matrix

@param A
Matrix to decompose. Overwritten with the decomposition
@return The current decomposition
196275
public DenseLU factor(DenseMatrix A) {
        singular = false;

        intW info = new intW(0);
        LAPACK.getInstance().dgetrf(A.numRows(), A.numColumns(), A.getData(),
                Matrices.ld(A.numRows()), piv, info);

        if (info.val > 0)
            singular = true;
        else if (info.val < 0)
            throw new IllegalArgumentException();

        LU.set(A);

        return this;
    }
443386
public static SymmDenseEVD factorize(Matrix A) throws NotConvergedException {
        return new SymmDenseEVD(A.numRows(), true)
                .factor(new UpperSymmDenseMatrix(A));
    }
276256
public static Matrix copy(Matrix m) {
		if (m==null) throw new NullPointerException();
		Matrix copy=createMatrix(m.getRowCount(), m.getColumnCount());
		for (int row=0; row<m.getRowCount(); row++) for (int col=0; col<m.getColumnCount(); col++)
			copy.set(row, col, m.getQuick(row, col));
		return copy;
	}

	/**
	 * Creates a copy of a matrix.
	 * @param m the matrix to copy
	 * @return a new matrix with the same dimensions and contents as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Matrix2 copy(Matrix2 m) {
		return (Matrix2) copy((Matrix)m);
	}

	/**
	 * Creates a copy of a matrix.
	 * @param m the matrix to copy
	 * @return a new matrix with the same dimensions and contents as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Matrix3 copy(Matrix3 m) {
		return (Matrix3) copy((Matrix)m);
	}

	/**
	 * Creates a copy of a matrix.
	 * @param m the matrix to copy
	 * @return a new matrix with the same dimensions and contents as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Vector copy(Vector m) {
		return (Vector) copy((Matrix)m);
	}

	/**
	 * Creates a copy of a matrix.
	 * @param m the matrix to copy
	 * @return a new matrix with the same dimensions and contents as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Vector2 copy(Vector2 m) {
		return (Vector2) copy((Matrix)m);
	}

	/**
	 * Creates a copy of a matrix.
	 * @param m the matrix to copy
	 * @return a new matrix with the same dimensions and contents as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Vector3 copy(Vector3 m) {
		return (Vector3) copy((Matrix)m);
	}

	/**
	 * Creates a copy of a matrix.
	 * @param m the matrix to copy
	 * @return a new matrix with the same dimensions and contents as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Vector4 copy(Vector4 m) {
		return (Vector4) copy((Matrix)m);
	}
	
	/**
	 * Creates a matrix with the same dimensions as the argument.
	 * @param matrix the template matrix
	 * @return a new zero matrix with the same dimensions as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Matrix createLike(Matrix matrix) {
		return createMatrix(matrix.getRowCount(), matrix.getColumnCount());
	}
	
	/**
	 * Creates a matrix with the same dimensions as the argument.
	 * @param matrix the template matrix
	 * @return a new zero matrix with the same dimensions as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Matrix2 createLike(Matrix2 matrix) {
		return (Matrix2) createMatrix(matrix.getRowCount(), matrix.getColumnCount());
	}
	
	/**
	 * Creates a matrix with the same dimensions as the argument.
	 * @param matrix the template matrix
	 * @return a new zero matrix with the same dimensions as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Matrix3 createLike(Matrix3 matrix) {
		return (Matrix3) createMatrix(matrix.getRowCount(), matrix.getColumnCount());
	}
	
	/**
	 * Creates a matrix with the same dimensions as the argument.
	 * @param matrix the template matrix
	 * @return a new zero matrix with the same dimensions as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Vector createLike(Vector matrix) {
		return (Vector) createMatrix(matrix.getRowCount(), matrix.getColumnCount());
	}
	
	/**
	 * Creates a matrix with the same dimensions as the argument.
	 * @param matrix the template matrix
	 * @return a new zero matrix with the same dimensions as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Vector2 createLike(Vector2 matrix) {
		return (Vector2) createMatrix(matrix.getRowCount(), matrix.getColumnCount());
	}
	
	/**
	 * Creates a matrix with the same dimensions as the argument.
	 * @param matrix the template matrix
	 * @return a new zero matrix with the same dimensions as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Vector3 createLike(Vector3 matrix) {
		return (Vector3) createMatrix(matrix.getRowCount(), matrix.getColumnCount());
	}
	
	/**
	 * Creates a matrix with the same dimensions as the argument.
	 * @param matrix the template matrix
	 * @return a new zero matrix with the same dimensions as the argument
	 * @throws NullPointerException if {@code matrix} is {@code null}
	 */
	public static Vector4 createLike(Vector4 matrix) {
		return (Vector4) createMatrix(matrix.getRowCount(), matrix.getColumnCount());
	}
	
	/**
	 * Creates an identity matrix of the specified dimension.
	 * @param size the dimension of the matrix
	 * @return a <code>size</code>×<code>size</code> identity matrix
	 * @throws IllegalArgumentException if the size is not positive
	 */
	public static Matrix identity(int size) {
		Matrix m = MatrixFactory.createMatrix(size, size);
		for (int i = 0; i < size; i++)
			m.set(i, i, 1);
		return m;
	}

	/**
	 * Creates a matrix of the specified size whose elements are all 1.
	 * @param rows the number of rows
	 * @param columns the number of columns
	 * @return a <code>rows</code>×<code>columns</code> matrix whose elements are all 1
	 * @throws IllegalArgumentException if either argument is non-positive
	 */
	public static Matrix ones(int rows, int columns) {
		Matrix m = MatrixFactory.createMatrix(rows, columns);
		for (int i = 0; i < m.getRowCount(); i++)
			for (int j = 0; j < m.getColumnCount(); j++)
				m.set(i, j, 1);
		return m;
	}

	/**
	 * Creates a matrix of the specified size whose elements are pseudorandom, uniformly distributed
	 * values between 0 (inclusive) and 1 (exclusive).
	 * @param rows the number of rows
	 * @param columns the number of columns
	 * @return a <code>rows</code>×<code>columns</code> matrix of random elements
	 * @throws IllegalArgumentException if either argument is non-positive
	 */
	public static Matrix random(int rows, int columns) {
		Matrix m = MatrixFactory.createMatrix(rows, columns);
		for (int i = 0; i < m.getRowCount(); i++)
			for (int j = 0; j < m.getColumnCount(); j++)
				m.set(i, j, Math.random());
		return m;
	}
	
	//
	// Matrix creation from arrays
	//
	
	/**
	 * Creates a matrix and initializes its elements from an array. The parameter should be a {@code double[r][c]}
	 * array with {@code r>0} and {@code c>0}. The matrix element at row {@code row} and column {@code col} will be
	 * {@code elements[row][col]}.
	 * @param elements the array containing the elements of the matrix
	 * @return the new matrix initialized from the array
	 * @throws NullPointerException if the array is {@code null}
	 * @throws IllegalArgumentException if the dimensions of the array are not valid
	 */
	public static Matrix createMatrix(double[][] elements) {
		Matrix m=createMatrix(elements.length, elements[0].length);
		for (int row=0; row<m.getRowCount(); row++) {
			double[] rowData=elements[row];
			if (rowData.length!=m.getColumnCount()) throw new IllegalArgumentException("Invalid array");
			for (int col=0; col<m.getColumnCount(); col++) {
				m.set(row, col, rowData[col]);
			}
		}
		return m;
	}
	
	/**
	 * Creates a matrix and initializes its elements from an array. The parameter should be a {@code double[2][2]}
	 * array. The matrix element at row {@code row} and column {@code col} will be
	 * {@code elements[row][col]}.
	 * @param elements the array containing the elements of the matrix
	 * @return the new matrix initialized from the array
	 * @throws NullPointerException if the array is {@code null}
	 * @throws IllegalArgumentException if the dimensions of the array are not correct
	 */
	public static Matrix2 createMatrix2(double[][] elements) {
		Matrix m=createMatrix(elements);
		if (m.getRowCount()!=2 || m.getColumnCount()!=2)
			throw new IllegalArgumentException("Array describes a matrix of size "+m.getRowCount()+"x"+m.getColumnCount()+" instead of 2x2");
		return (Matrix2)m;
	}
	
	/**
	 * Creates a matrix and initializes its elements from an array. The parameter should be a {@code double[3][3]}
	 * array. The matrix element at row {@code row} and column {@code col} will be
	 * {@code elements[row][col]}.
	 * @param elements the array containing the elements of the matrix
	 * @return the new matrix initialized from the array
	 * @throws NullPointerException if the array is {@code null}
	 * @throws IllegalArgumentException if the dimensions of the array are not correct
	 */
	public static Matrix3 createMatrix3(double[][] elements) {
		Matrix m=createMatrix(elements);
		if (m.getRowCount()!=3 || m.getColumnCount()!=3)
			throw new IllegalArgumentException("Array describes a matrix of size "+m.getRowCount()+"x"+m.getColumnCount()+" instead of 3x3");
		return (Matrix3)m;
	}
	
	/**
	 * Creates a vector and initializes its elements from an array.
	 * @param elements the array containing the elements of the vector
	 * @return the new vector initialized from the array
	 * @throws NullPointerException if the array is {@code null}
	 * @throws IllegalArgumentException if the length of the array is zero
	 */
	public static Vector createVector(double[] elements) {
		Vector v=createVector(elements.length);
		for (int i=0; i<elements.length; i++)
			v.setCoordQuick(i, elements[i]);
		return v;
	}
	
	/**
	 * Creates a matrix backed by the specified data structure.
	 * @param core the data structure for the new matrix
	 * @return a matrix of proper type using the specified matrix core for its implementation
	 * @throws NullPointerException if {@code core} is {@code null}
	 */
	public static Matrix create(MatrixCore core) {
		if (core.getColumnCount()==1) {
			switch (core.getRowCount()) {
				case 2: return new CoredVector2Impl(core);
				case 3: return new CoredVector3Impl(core);
				case 4: return new CoredVector4Impl(core);
				default: return new CoredVectorImpl(core);
			}
		} else {
			if (core.getRowCount()==2 && core.getColumnCount()==2)
				return new CoredMatrix2Impl(core);
			if (core.getRowCount()==3 && core.getColumnCount()==3)
				return new CoredMatrix3Impl(core);
			return new CoredMatrixImpl(core);
		}
	}
}
630169
public static DMatrixRMaj extractColumn(DMatrixRMaj a , int column , DMatrixRMaj out ) {
        if( out == null)
            out = new DMatrixRMaj(a.numRows,1);
        else if( !MatrixFeatures_DDRM.isVector(out) || out.getNumElements() != a.numRows )
            throw new MatrixDimensionException("Output must be a vector of length "+a.numRows);

        int index = column;
        for (int i = 0; i < a.numRows; i++, index += a.numCols ) {
            out.data[i] = a.data[index];
        }
        return out;
    }
860019
public void decompose( DMatrixRMaj A ) {

        this.QR = A.copy();

        int N = Math.min(A.numCols,A.numRows);

        gammas = new double[ A.numCols ];

        DMatrixRMaj A_small = new DMatrixRMaj(A.numRows,A.numCols);
        DMatrixRMaj A_mod = new DMatrixRMaj(A.numRows,A.numCols);
        DMatrixRMaj v = new DMatrixRMaj(A.numRows,1);
        DMatrixRMaj Q_k = new DMatrixRMaj(A.numRows,A.numRows);

        for( int i = 0; i < N; i++ ) {
            // reshape temporary variables
            A_small.reshape(QR.numRows-i,QR.numCols-i);
            v.reshape(A_small.numRows,1);
            Q_k.reshape(v.getNumElements(),v.getNumElements());

            // use extract matrix to get the column that is to be zeroed
            CommonOps_DDRM.extract(QR,i,QR.numRows,i,i+1,v,0,0);

            double max = CommonOps_DDRM.elementMaxAbs(v);

            if( max > 0 && v.getNumElements() > 1 ) {
                // normalize to reduce overflow issues
                CommonOps_DDRM.divide(v,max);

                // compute the magnitude of the vector
                double tau = NormOps_DDRM.normF(v);

                if( v.get(0) < 0 )
                    tau *= -1.0;

                double u_0 = v.get(0) + tau;
                double gamma = u_0/tau;

                CommonOps_DDRM.divide(v,u_0);
                v.set(0,1.0);

                // extract the submatrix of A which is being operated on
                CommonOps_DDRM.extract(QR,i,QR.numRows,i,QR.numCols,A_small,0,0);

                // A = (I - &gamma;*u*u<sup>T</sup>)A
                CommonOps_DDRM.setIdentity(Q_k);
                CommonOps_DDRM.multAddTransB(-gamma,v,v,Q_k);
                CommonOps_DDRM.mult(Q_k,A_small,A_mod);

                // save the results
                CommonOps_DDRM.insert(A_mod, QR, i,i);
                CommonOps_DDRM.insert(v, QR, i,i);
                QR.unsafe_set(i,i,-tau*max);

                // save gamma for recomputing Q later on
                gammas[i] = gamma;
            }
        }
    }
343852
public Matrix getL () {
		Matrix X = new Matrix(m,n);
		double[][] L = X.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				if (i > j) {
					L[i][j] = LU[i][j];
				} else if (i == j) {
					L[i][j] = 1.0;
				} else {
					L[i][j] = 0.0;
				}
			}
		}
		return X;
	}
1157354
public static DenseMatrix matrix(double[] A) {
        if (nlmatrixZeros != null) {
            try {
                return (DenseMatrix) nlmatrixArray.newInstance(A);
            } catch (Exception e) {
                logger.error("Failed to call NLMatrix(double[]): {}", e);
            }
        }

        return new JMatrix(A);
    }
123512
public static Matrix diag(Vec a)
    {
        DenseMatrix A = new DenseMatrix(a.length(), a.length());
        for(Iterator<IndexValue> iter = a.getNonZeroIterator(); iter.hasNext();)
        {
            IndexValue iv = iter.next();
            A.set(iv.getIndex(), iv.getIndex(), iv.getValue());
        }
            
        return A;
    }
958542
@Override
    public Matrix[] decompose() {

        Matrix l = matrix.blankOfShape(matrix.rows(), matrix.rows());

        for (int j = 0; j < l.rows(); j++) {

            double d = 0.0;

            for (int k = 0; k < j; k++) {

                double s = 0.0;

                for (int i = 0; i < k; i++) {
                    s += l.get(k, i) * l.get(j, i);
                }

                s = (matrix.get(j, k) - s) / l.get(k, k);

                l.set(j, k, s);

                d = d + s * s;
            }

            d = matrix.get(j, j) - d;

            l.set(j, j, Math.sqrt(Math.max(d, 0.0)));

            for (int k = j + 1; k < l.rows(); k++) {
                l.set(j, k, 0.0);
            }
        }

        return new Matrix[] { l };
    }
43520
public static DoubleMatrix cholesky(DoubleMatrix A) {
      DoubleMatrix result = A.dup();
      int info = NativeBlas.dpotrf('U', A.rows, result.data, 0, A.rows);
      if (info < 0) {
          throw new LapackArgumentException("DPOTRF", -info);
      } else if (info > 0) {
          throw new LapackPositivityException("DPOTRF", "Minor " + info + " was negative. Matrix must be positive definite.");
      }
      clearLower(result);
      return result;
  }
480990
function lu(A) {
    var size = A.length;
    //var L=jStat.diagonal(jStat.ones(1,size)[0]);
    var L = jStat.identity(size);
    var R = jStat.zeros(A.length, A[0].length);
    var parts;
    jStat.arange(size).forEach(function(t) {
      R[0][t] = A[0][t];
    });
    jStat.arange(1, size).forEach(function(l) {
      jStat.arange(l).forEach(function(i) {
        parts = jStat.arange(i).map(function(jj) {
          return L[l][jj] * R[jj][i];
        });
        L[l][i] = (A[l][i] - jStat.sum(parts)) / R[i][i];
      });
      jStat.arange(l, size).forEach(function(j) {
        parts = jStat.arange(l).map(function(jj) {
          return L[l][jj] * R[jj][j];
        });
        R[l][j] = A[i][j] - jStat.sum(parts);
      });
    });
    return [L, R];
  }
1221438
public static QR factorize(Matrix A) {
        return new QR(A.numRows(), A.numColumns()).factor(new DenseMatrix(A));
    }
928518
@Nonnull
  @ReturnsMutableCopy
  public Matrix solve (@Nonnull final Matrix aMatrix)
  {
    if (aMatrix.getRowDimension () != m_nRows)
      throw new IllegalArgumentException ("Matrix row dimensions must agree.");
    if (!isNonsingular ())
      throw new IllegalStateException ("Matrix is singular.");

    // Copy right hand side with pivoting
    final int nCols = aMatrix.getColumnDimension ();
    final Matrix aNewMatrix = aMatrix.getMatrix (m_aPivot, 0, nCols - 1);
    final double [] [] aNewArray = aNewMatrix.internalGetArray ();

    // Solve L*Y = B(piv,:)
    for (int k = 0; k < m_nCols; k++)
    {
      final double [] aNewk = aNewArray[k];
      for (int i = k + 1; i < m_nCols; i++)
      {
        final double [] aLUi = m_aLU[i];
        final double [] aNewi = aNewArray[i];
        for (int j = 0; j < nCols; j++)
          aNewi[j] -= aNewk[j] * aLUi[k];
      }
    }
    // Solve U*X = Y;
    for (int k = m_nCols - 1; k >= 0; k--)
    {
      final double [] aLUk = m_aLU[k];
      final double [] aNewk = aNewArray[k];
      for (int j = 0; j < nCols; j++)
        aNewk[j] /= aLUk[k];
      for (int i = 0; i < k; i++)
      {
        final double [] aLUi = m_aLU[i];
        final double [] aNewi = aNewArray[i];
        for (int j = 0; j < nCols; j++)
          aNewi[j] -= aNewk[j] * aLUi[k];
      }
    }
    return aNewMatrix;
  }
246255
public static DoubleMatrix[] symmetricEigenvectors(DoubleMatrix A) {
        A.assertSquare();
        DoubleMatrix eigenvalues = new DoubleMatrix(A.rows);
        DoubleMatrix eigenvectors = A.dup();
        int isuppz[] = new int[2 * A.rows];
        SimpleBlas.syevr('V', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, eigenvectors, isuppz);
        return new DoubleMatrix[]{eigenvectors, DoubleMatrix.diag(eigenvalues)};
    }
1088348
@Override
    public Matrix[] decompose() {

        // AHTIUNG: this code derived from Jama

        Matrix a = matrix.copy();

        int n = Math.min(a.rows(), a.columns());

        Matrix u = matrix.blankOfShape(a.rows(), n);
        Matrix s = matrix.blankOfShape(a.columns(), a.columns());
        Matrix v = matrix.blankOfShape(a.columns(), a.columns());

        Vector e = DenseVector.zero(a.columns());
        Vector work = DenseVector.zero(a.rows());

        int nct = Math.min(a.rows() - 1, a.columns());
        int nrt = Math.max(0, Math.min(a.columns() - 2, a.rows()));

        for (int k = 0; k < Math.max(nct, nrt); k++) {

            if (k < nct) {

                s.set(k, k, 0.0);

                for (int i = k; i < a.rows(); i++) {
                    s.set(k, k, hypot(s.get(k, k), a.get(i, k)));
                }

                if (Math.abs(s.get(k, k)) > Matrices.EPS) {

                    if (a.get(k, k) < 0.0) {
                        s.updateAt(k, k, Matrices.INV_FUNCTION);
                    }

                    double skk = s.get(k, k);
                    for (int i = k; i < a.rows(); i++) {
                        a.updateAt(i, k, Matrices.asDivFunction(skk));
                    }

                    a.updateAt(k, k, Matrices.INC_FUNCTION);
                }

                s.updateAt(k, k, Matrices.INV_FUNCTION);
            }

            for (int j = k + 1; j < a.columns(); j++) {

                if ((k < nct) & (Math.abs(s.get(k, k)) > Matrices.EPS)) {

                    double t = 0;

                    for (int i = k; i < a.rows(); i++) {
                        t += a.get(i, k) * a.get(i, j);
                    }

                    t = -t / a.get(k, k);

                    for (int i = k; i < a.rows(); i++) {
                        a.updateAt(i, j, Matrices.asPlusFunction(t * a.get(i, k)));
                    }
                }

                e.set(j, a.get(k, j));
            }

            if (k < nct) {

                for (int i = k; i < a.rows(); i++) {
                    u.set(i, k, a.get(i, k));
                }

            }

            if (k < nrt) {

                e.set(k, 0);

                for (int i = k + 1; i < a.columns(); i++) {
                    e.set(k, hypot(e.get(k), e.get(i)));
                }

                if (Math.abs(e.get(k)) > Matrices.EPS) {

                    if (e.get(k + 1) < 0.0) {

                        e.updateAt(k, Vectors.INV_FUNCTION);
                    }

                    double ek = e.get(k);
                    for (int i = k + 1; i < a.columns(); i++) {
                        e.updateAt(i, Vectors.asDivFunction(ek));
                    }

                    e.updateAt(k + 1, Vectors.INC_FUNCTION);
                }

                e.updateAt(k, Vectors.INV_FUNCTION);

                if ((k + 1 < a.rows()) && (Math.abs(e.get(k)) > Matrices.EPS)) {

                    for (int i = k + 1; i < a.rows(); i++) {
                        work.set(i, 0.0);
                    }

                    for (int j = k + 1; j < a.columns(); j++) {
                        for (int i = k + 1; i < a.rows(); i++) {
                            work.updateAt(i, Vectors.asPlusFunction(e.get(j) *
                                a.get(i, j)));
                        }
                    }

                    for (int j = k + 1; j < a.columns(); j++) {

                        double t = -e.get(j) / e.get(k + 1);

                        for (int i = k + 1; i < a.rows(); i++) {
                            a.updateAt(i, j, Matrices.asPlusFunction(t *
                                work.get(i)));
                        }
                    }
                }

                for (int i = k + 1; i < a.columns(); i++) {
                    v.set(i, k, e.get(i));
                }
            }
        }

        int p = Math.min(a.columns(), a.rows() + 1);

        if (nct < a.columns()) {
            s.set(nct, nct, a.get(nct, nct));
        }

        if (a.rows() < p) {
            s.set(p - 1, p - 1, 0.0);
        }

        if (nrt + 1 < p) {
            e.set(nrt, a.get(nrt, p - 1));
        }

        e.set(p - 1, 0.0);

        for (int j = nct; j < n; j++) {

            for (int i = 0; i < a.rows(); i++) {
                u.set(i, j, 0.0);
            }

            u.set(j, j, 1.0);
        }

        for (int k = nct - 1; k >= 0; k--) {

            if (Math.abs(s.get(k, k)) > Matrices.EPS) {

                for (int j = k + 1; j < n; j++) {

                    double t = 0;
                    for (int i = k; i < a.rows(); i++) {
                        t += u.get(i, k) * u.get(i, j);
                    }

                    t = -t / u.get(k, k);

                    for (int i = k; i < a.rows(); i++) {
                        u.updateAt(i, j, Matrices.asPlusFunction(t * u.get(i, k)));
                    }
                }

                for (int i = k; i < a.rows(); i++) {
                    u.updateAt(i, k, Matrices.INV_FUNCTION);
                }

                u.updateAt(k, k, Matrices.INC_FUNCTION);

                for (int i = 0; i < k - 1; i++) {
                    u.set(i, k, 0.0);
                }

            } else {

                for (int i = 0; i < a.rows(); i++) {
                    u.set(i, k, 0.0);
                }

                u.set(k, k, 1.0);
            }
        }

        for (int k = n - 1; k >= 0; k--) {

            if ((k < nrt) & (Math.abs(e.get(k)) > Matrices.EPS)) {

                for (int j = k + 1; j < n; j++) {

                    double t = 0;

                    for (int i = k + 1; i < a.columns(); i++) {
                        t += v.get(i, k) * v.get(i, j);
                    }

                    t = -t / v.get(k + 1, k);

                    for (int i = k + 1; i < a.columns(); i++) {
                        v.updateAt(i, j, Matrices.asPlusFunction(t * v.get(i, k)));
                    }
                }
            }

            for (int i = 0; i < a.columns(); i++) {
                v.set(i, k, 0.0);
            }

            v.set(k, k, 1.0);
        }

        int pp = p - 1;
        int iter = 0;
        double eps = Math.pow(2.0, -52.0);
        double tiny = Math.pow(2.0, -966.0);

        while (p > 0) {

            int k;
            int kase;

            for (k = p - 2; k >= -1; k--) {
                if (k == -1)
                    break;

                if (Math.abs(e.get(k)) <= tiny + eps * 
                        (Math.abs(s.get(k, k)) + Math.abs(s.get(k + 1, k + 1)))) {

                    e.set(k, 0.0);
                    break;
                }
            }

            if (k == p - 2) {

                kase = 4;

            } else {

                int ks;

                for (ks = p - 1; ks >= k; ks--) {

                    if (ks == k)
                        break;

                    double t = (ks != p ? Math.abs(e.get(ks)) : 0.)
                            + (ks != k + 1 ? Math.abs(e.get(ks - 1)) 
                            			   : 0.);

                    if (Math.abs(s.get(ks, ks)) <= tiny + eps * t) {
                        s.set(ks, ks, 0.0);
                        break;
                    }
                }

                if (ks == k) {
                    kase = 3;
                } else if (ks == p - 1) {
                    kase = 1;
                } else {
                    kase = 2;
                    k = ks;
                }
            }

            k++;

            switch (kase) {

            case 1: {
                double f = e.get(p - 2);
                e.set(p - 2, 0.0);

                for (int j = p - 2; j >= k; j--) {

                    double t = hypot(s.get(j, j), f);
                    double cs = s.get(j, j) / t;
                    double sn = f / t;

                    s.set(j, j, t);

                    if (j != k) {
                        f = -sn * e.get(j - 1);
                        e.set(j - 1, cs * e.get(j - 1));
                    }

                    for (int i = 0; i < a.columns(); i++) {
                        t = cs * v.get(i, j) + sn 
                        		* v.get(i, p - 1);
                        v.set(i, p - 1,
                                -sn * v.get(i, j) 
                                + cs * v.get(i, p - 1));
                        v.set(i, j, t);
                    }
                }
            }
                break;

            case 2: {
                double f = e.get(k - 1);
                e.set(k - 1, 0.0);

                for (int j = k; j < p; j++) {

                    double t = hypot(s.get(j, j), f);
                    double cs = s.get(j, j) / t;
                    double sn = f / t;

                    s.set(j, j, t);
                    f = -sn * e.get(j);
                    e.set(j, cs * e.get(j));

                    for (int i = 0; i < a.rows(); i++) {
                        t = cs * u.get(i, j) 
                        		+ sn * u.get(i, k - 1);
                        u.set(i, k - 1,
                                -sn * u.get(i, j) 
                                + cs * u.get(i, k - 1));
                        u.set(i, j, t);
                    }
                }
            }
                break;

            case 3: {

                double scale = Math.max(
                        Math.max(Math.max(
                                Math.max(Math.abs(s.get(p - 1, p - 1)),
                                        Math.abs(s.get(p - 2, p - 2))),
                                    Math.abs(e.get(p - 2))), 
                                Math.abs(s.get(k, k))),
                        Math.abs(e.get(k)));

                double sp = s.get(p - 1, p - 1) / scale;
                double spm1 = s.get(p - 2, p - 2) / scale;
                double epm1 = e.get(p - 2) / scale;
                double sk = s.get(k, k) / scale;
                double ek = e.get(k) / scale;
                double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
                double c = (sp * epm1) * (sp * epm1);
                double shift = 0.0;

                if ((b != 0.0) | (c != 0.0)) {
                    shift = Math.sqrt(b * b + c);
                    if (b < 0.0) {
                        shift = -shift;
                    }
                    shift = c / (b + shift);
                }

                double f = (sk + sp) * (sk - sp) + shift;
                double g = sk * ek;

                for (int j = k; j < p - 1; j++) {
                    double t = hypot(f, g);
                    double cs = f / t;
                    double sn = g / t;

                    if (j != k) {
                        e.set(j - 1, t);
                    }

                    f = cs * s.get(j, j) + sn * e.get(j);
                    e.set(j, 
                    		cs * e.get(j) - sn * s.get(j, j));
                    g = sn * s.get(j + 1, j + 1);
                    s.set(j + 1, j + 1, 
                    		cs * s.get(j + 1, j + 1));

                    for (int i = 0; i < a.columns(); i++) {
                        t = cs * v.get(i, j) 
                        		+ sn * v.get(i, j + 1);
                        v.set(i, j + 1,
                                -sn * v.get(i, j) 
                                + cs * v.get(i, j + 1));
                        v.set(i, j, t);
                    }

                    t = hypot(f, g);
                    cs = f / t;
                    sn = g / t;
                    s.set(j, j, t);
                    f = cs * e.get(j) + sn * s.get(j + 1, j + 1);
                    s.set(j + 1, j + 1, 
                            -sn * e.get(j) 
                            + cs * s.get(j + 1, j + 1));
                    g = sn * e.get(j + 1);
                    e.updateAt(j + 1, Vectors.asMulFunction(cs));

                    if (j < a.rows() - 1) {
                        for (int i = 0; i < a.rows(); i++) {
                            t = cs * u.get(i, j) 
                                    + sn * u.get(i, j + 1);
                            u.set(i, j + 1,
                                    -sn * u.get(i, j) 
                                    + cs * u.get(i, j + 1));
                            u.set(i, j, t);
                        }
                    }
                }

                e.set(p - 2, f);
                iter = iter + 1;
            }
                break;

            case 4: {

                if (s.get(k, k) <= 0.0) {
                    s.set(k, k, s.get(k, k) < 0.0 ? -s.get(k, k) : 0.0);
                    for (int i = 0; i <= pp; i++) {
                        v.updateAt(i, k, Matrices.INV_FUNCTION);
                    }
                }

                while (k < pp) {

                    if (s.get(k, k) >= s.get(k + 1, k + 1)) {
                        break;
                    }

                    double t = s.get(k, k);
                    s.set(k, k, s.get(k + 1, k + 1));
                    s.set(k + 1, k + 1, t);

                    if (k < a.columns() - 1) {
                        for (int i = 0; i < a.columns(); i++) {
                            t = v.get(i, k + 1);
                            v.set(i, k + 1, v.get(i, k));
                            v.set(i, k, t);
                        }
                    }

                    if (k < a.rows() - 1) {
                        for (int i = 0; i < a.rows(); i++) {
                            t = u.get(i, k + 1);
                            u.set(i, k + 1, u.get(i, k));
                            u.set(i, k, t);
                        }
                    }

                    k++;
                }

                iter = 0;
                p--;
            }
                break;
            }
        }

        // TODO:
        // We need to change the logic of this code in order do not use resize
        // It is also not a good idea to use SVN for [m < n] matrices
        // We need to change it in further releases
        return new Matrix[] { u, s.copyOfShape(n, a.columns()), v };
    }
439845
public static SimpleMatrix convertSimple( DMatrixRBlock A ) {
        DMatrixRMaj B = ConvertDMatrixStruct.convert(A, null);

        return SimpleMatrix.wrap(B);
    }
964415
@Override
    public boolean setA(DMatrixRMaj A) {
        blockA.reshape(A.numRows,A.numCols,false);
        MatrixOps_DDRB.convert(A,blockA);

        return alg.setA(blockA);
    }