Cholesky decomposition for symmetric and positive definite matrix.
Only the lower triangular part will be used in the decomposition.

@throws IllegalArgumentException if the matrix is not positive definite.
345030
@Override
    public Cholesky cholesky() {
        if (nrows() != ncols()) {
            throw new UnsupportedOperationException("Cholesky decomposition on non-square matrix");
        }

        int n = nrows();

        // Main loop.
        for (int j = 0; j < n; j++) {
            double d = 0.0;
            for (int k = 0; k < j; k++) {
                double s = 0.0;
                for (int i = 0; i < k; i++) {
                    s += get(k, i) * get(j, i);
                }
                s = (get(j, k) - s) / get(k, k);
                set(j, k, s);
                d = d + s * s;
            }
            d = get(j, j) - d;

            if (d < 0.0) {
                throw new IllegalArgumentException("The matrix is not positive definite.");
            }

            set(j, j, Math.sqrt(d));
        }

        return new Cholesky(this);
    }
459264
@Override
	public double calculatePhase1InitialFeasiblePoint(double[] originalNotFeasiblePoint, double tolerance){
		RealMatrix F = this.buildS(originalNotFeasiblePoint).scalarMultiply(-1);
		RealMatrix S = F.scalarMultiply(-1);
		try{
			new CholeskyDecomposition(S);
			//already feasible
			return -1;
		}catch(NonPositiveDefiniteMatrixException ee){
			//it does NOT mean that F is negative, it can be not definite
			EigenDecomposition eFact = new EigenDecomposition(F);
			double[] eValues = eFact.getRealEigenvalues();
			double minEigenValue = eValues[Utils.getMinIndex(DoubleFactory1D.dense.make(eValues))];
			return  -Math.min(minEigenValue * Math.pow(tolerance, -0.5), 0.);
		}
	}
617974
@Override
    public Matrix[] decompose() {

        if (matrix.is(Matrices.SYMMETRIC_MATRIX)) {
            return decomposeSymmetricMatrix(matrix);
        } else if (matrix.rows() == matrix.columns()) {
            return decomposeNonSymmetricMatrix(matrix);
        } else {
            throw new IllegalArgumentException("Can't decompose rectangle matrix");
        }
    }
435553
public void gaussianElimination(int[][] matrix) {
    assert(matrix != null && matrix.length > 0 && matrix[0].length > 0
           && matrix.length < matrix[0].length);
    int height = matrix.length;
    int width = matrix[0].length;
    for (int i = 0; i < height; i++) {
      boolean pivotFound = false;
      // scan the column for a nonzero pivot and swap it to the diagonal
      for (int j = i; j < height; j++) {
        if (matrix[i][j] != 0) {
          int[] tmp = matrix[i];
          matrix[i] = matrix[j];
          matrix[j] = tmp;
          pivotFound = true;
          break;
        }
      }
      if (!pivotFound) {
        continue;
      }
      int pivot = matrix[i][i];
      for (int j = i; j < width; j++) {
        matrix[i][j] = divide(matrix[i][j], pivot);
      }
      for (int j = i + 1; j < height; j++) {
        int lead = matrix[j][i];
        for (int k = i; k < width; k++) {
          matrix[j][k] = add(matrix[j][k], multiply(lead, matrix[i][k]));
        }
      }
    }
    for (int i = height - 1; i >=0; i--) {
      for (int j = 0; j < i; j++) {
        int lead = matrix[j][i];
        for (int k = i; k < width; k++) {
          matrix[j][k] = add(matrix[j][k], multiply(lead, matrix[i][k]));
        }
      }
    }
  }
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;
            }
        }
    }
480943
@Override
    public boolean setA(DMatrixRBlock A) {
        // Extract a lower triangular solution
        if( !decomposer.decompose(A) )
            return false;

        blockLength = A.blockLength;

        return true;
    }
491081
private boolean _decompose() {
        double h[] = QH.data;

        for( int k = 0; k < N-2; k++ ) { // k = column
            u[k*2] = 0;
            u[k*2+1] = 0;
            double max = QrHelperFunctions_ZDRM.extractColumnAndMax(QH,k+1,N,k,u,0);

            if( max > 0 ) {
                // -------- set up the reflector Q_k

                double gamma = QrHelperFunctions_ZDRM.computeTauGammaAndDivide(k+1,N,u,max,tau);
                gammas[k] = gamma;

                // divide u by u_0
                double real_u_0 = u[(k+1)*2]   + tau.real;
                double imag_u_0 = u[(k+1)*2+1] + tau.imaginary;
                QrHelperFunctions_ZDRM.divideElements(k + 2, N, u, 0, real_u_0,imag_u_0 );

                // write the reflector into the lower left column of the matrix
                for (int i = k+2; i < N; i++) {
                    h[(i*N+k)*2]   = u[i*2];
                    h[(i*N+k)*2+1] = u[i*2+1];
                }

                u[(k+1)*2]   = 1;
                u[(k+1)*2+1] = 0;

                // ---------- multiply on the left by Q_k
                QrHelperFunctions_ZDRM.rank1UpdateMultR(QH, u,0, gamma, k + 1, k + 1, N, b);

                // ---------- multiply on the right by Q_k
                QrHelperFunctions_ZDRM.rank1UpdateMultL(QH, u,0, gamma, 0, k + 1, N);

                // since the first element in the householder vector is known to be 1
                // store the full upper hessenberg
                h[((k+1)*N+k)*2]   = -tau.real*max;
                h[((k+1)*N+k)*2+1] = -tau.imaginary*max;

            } else {
                gammas[k] = 0;
            }

        }

        return true;
    }
838408
public Cholesky CholeskyOfAtA() {
        DenseMatrix VD = Matrix.zeros(V.nrows(), V.ncols());
        for (int i = 0; i < V.nrows(); i++) {
            for (int j = 0; j < V.ncols(); j++) {
                VD.set(i, j, V.get(i, j) * s[j]);
            }
        }

        return new Cholesky(VD.aat());
    }
241819
public static void symmLowerToFull( DMatrixSparseCSC A , DMatrixSparseCSC B , @Nullable IGrowArray gw )
    {
        if( A.numCols != A.numRows )
            throw new IllegalArgumentException("Must be a lower triangular square matrix");

        int N = A.numCols;
        int w[] = adjust(gw,N,N);
        B.reshape(N,N,A.nz_length*2);
        B.indicesSorted = false;

        //=== determine the row counts of the full matrix
        for (int col = 0; col < N; col++) {
            int idx0 = A.col_idx[col];
            int idx1 = A.col_idx[col+1];

            // We know the length of the lower part of this column already
            w[col] += idx1 - idx0;

            // add elements to the top of the other columns along row with index 'col'
            for (int i = idx0; i < idx1; i++) {
                int row = A.nz_rows[i];
                if( row > col ) {
                    w[row]++;
                }
            }
        }

        // Update the structure of B
        B.histogramToStructure(w);

        // Zero W again. It's being used to keep track of how many elements have been added to a column already
        Arrays.fill(w,0,N,0);
        // Fill in matrix
        for (int col = 0; col < N; col++) {

            int idx0 = A.col_idx[col];
            int idx1 = A.col_idx[col+1];

            int lengthA = idx1 - idx0;
            int lengthB = B.col_idx[col+1] - B.col_idx[col];

            // Copy the non-zero values from A into B along the columns while taking in account the upper
            // elements already copied
            System.arraycopy(A.nz_values,idx0,B.nz_values,B.col_idx[col]+lengthB-lengthA,lengthA);
            System.arraycopy(A.nz_rows,idx0,B.nz_rows,B.col_idx[col]+lengthB-lengthA,lengthA);

            // Copy this column into the upper portion of B
            for (int i = idx0; i < idx1; i++) {
                int row = A.nz_rows[i];
                if( row > col ) {
                    int indexB = B.col_idx[row] + w[row]++;
                    B.nz_rows[indexB] = col;
                    B.nz_values[indexB] = A.nz_values[i];
                }
            }
        }

    }
906542
public static boolean cholU( DMatrix4x4 A ) {

        A.a11 = Math.sqrt(A.a11);
        A.a21 = 0;
        A.a31 = 0;
        A.a41 = 0;
        A.a12 = (A.a12)/A.a11;
        A.a22 = Math.sqrt(A.a22-A.a12*A.a12);
        A.a32 = 0;
        A.a42 = 0;
        A.a13 = (A.a13)/A.a11;
        A.a23 = (A.a23-A.a12*A.a13)/A.a22;
        A.a33 = Math.sqrt(A.a33-A.a13*A.a13-A.a23*A.a23);
        A.a43 = 0;
        A.a14 = (A.a14)/A.a11;
        A.a24 = (A.a24-A.a12*A.a14)/A.a22;
        A.a34 = (A.a34-A.a13*A.a14-A.a23*A.a24)/A.a33;
        A.a44 = Math.sqrt(A.a44-A.a14*A.a14-A.a24*A.a24-A.a34*A.a34);
        return !UtilEjml.isUncountable(A.a44);
    }
737655
private void updateCholesky() {
    // TODO: further improve handling of degenerated cases?
    CholeskyDecomposition chol = new CholeskyDecomposition(covariance);
    if(!chol.isSPD()) {
      // Add a small value to the diagonal, to reduce some rounding problems.
      double s = 0.;
      for(int i = 0; i < covariance.length; i++) {
        s += covariance[i][i];
      }
      s *= SINGULARITY_CHEAT / covariance.length;
      for(int i = 0; i < covariance.length; i++) {
        covariance[i][i] += s;
      }
      chol = new CholeskyDecomposition(covariance);
    }
    if(!chol.isSPD()) {
      LOG.warning("A cluster has degenerated, likely due to lack of variance in a subset of the data or too extreme magnitude differences.\n" + //
          "The algorithm will likely stop without converging, and fail to produce a good fit.");
      chol = this.chol != null ? this.chol : chol; // Prefer previous
    }
    this.chol = chol;
    logNormDet = FastMath.log(weight) - .5 * logNorm - getHalfLogDeterminant(this.chol);
  }
784154
@Override
	public DoubleMatrix1D[] solve() throws Exception {

		final CholeskyFactorization HFact = new CholeskyFactorization(H, new Matrix1NornRescaler());
		boolean isHReducible = true;
		try{
			HFact.factorize();
		}catch(Exception e){
			isHReducible = false;
		}


        final DoubleMatrix1D v;// dim equals cols of A
        final DoubleMatrix1D w;// dim equals rank of A

		if (isHReducible) {
			// Solving KKT system via elimination
			DoubleMatrix1D HInvg = HFact.solve(g);
			if (A != null) {
				DoubleMatrix2D HInvAT = HFact.solve(AT);
				DoubleMatrix2D MenoSLower = ColtUtils.subdiagonalMultiply(A, HInvAT);
				DoubleMatrix1D AHInvg = ALG.mult(A, HInvg);
				
				CholeskyFactorization MSFact = new CholeskyFactorization(MenoSLower, new Matrix1NornRescaler());
				MSFact.factorize();
				if(h == null){
					w = MSFact.solve(ColtUtils.scalarMult(AHInvg, -1));
				}else{
					w = MSFact.solve(ColtUtils.add(h, AHInvg, -1));
				}
				
				v = HInvg.assign(ALG.mult(HInvAT, w), Functions.plus).assign(Mult.mult(-1));
			} else {
				w = null;
				v = HInvg.assign(Mult.mult(-1));
			}
		} else {
			// H is singular
			// Solving the full KKT system
			if(A!=null){
				DoubleMatrix1D[] fullSol =  this.solveAugmentedKKT();
				v = fullSol[0];
				w = fullSol[1];
			}else{
				//@TODO: try with rescaled H
				throw new Exception("KKT solution failed");
			}
		}

		// solution checking
		if (this.checkKKTSolutionAccuracy && !this.checkKKTSolutionAccuracy(v, w)) {
			log.error("KKT solution failed");
			throw new Exception("KKT solution failed");
		}

		DoubleMatrix1D[] ret = new DoubleMatrix1D[2];
		ret[0] = v;
		ret[1] = w;
		return ret;
	}
384351
protected void decompose() throws MatrixException
    {
        // Return if the decomposition is valid.
        if (LU != null) return;

        // Create a new LU matrix and permutation vector.
        // LU is initially just a copy of the values of this system.
        LU = new SquareMatrix(this.copyValues2D());
        permutation = new int[nRows];

        double scales[] = new double[nRows];

        // Loop to initialize the permutation vector and scales.
        for (int r = 0; r < nRows; ++r) {
            permutation[r] = r;     // initially no row exchanges

            // Find the largest row element.
            double largestRowElmt = 0;
            for (int c = 0; c < nRows; ++c) {
                double elmt = Math.abs(LU.at(r, c));
                if (largestRowElmt < elmt) largestRowElmt = elmt;
            }

            // Set the scaling factor for row equilibration.
            if (largestRowElmt != 0) {
                scales[r] = 1/largestRowElmt;
            }
            else {
                throw new MatrixException(MatrixException.ZERO_ROW);
            }
        }

        // Do forward elimination with scaled partial row pivoting.
        forwardElimination(scales);

        // Check bottom right element of the permuted matrix.
        if (LU.at(permutation[nRows - 1], nRows - 1) == 0) {
            throw new MatrixException(MatrixException.SINGULAR);
        }
    }
776373
@Override
	public void computeHessian(DMatrixRMaj jacLeft , DMatrixRMaj jacRight) {
		A.reshape(jacLeft.numCols,jacLeft.numCols);
		B.reshape(jacLeft.numCols,jacRight.numCols);
		D.reshape(jacRight.numCols,jacRight.numCols);

		// take advantage of the inner product's symmetry when possible to reduce
		// the number of calculations
		MatrixMultProduct_DDRM.inner_reorder_lower(jacLeft,A);
		CommonOps_DDRM.symmLowerToFull(A);
		CommonOps_DDRM.multTransA(jacLeft,jacRight,B);
		MatrixMultProduct_DDRM.inner_reorder_lower(jacRight,D);
		CommonOps_DDRM.symmLowerToFull(D);
	}
367879
public void multT(DoubleMatrix bMat, DoubleMatrix cMat, boolean transposeA, boolean transposeB) {
        if (!transposeA) {            
            throw new IllegalArgumentException("Multiplication where this matrix NOT transposed is not implemented.");
        }
        
        SparseRowDoubleMatrix.checkMultDimensions(this.viewTranspose(), transposeB ? bMat.viewTranspose() : bMat, cMat);

        int numRowsC = cMat.getNumRows();
        int numColsC = cMat.getNumColumns();
        
        if (bMat instanceof DenseDoubleMatrix) {
            DenseDoubleMatrix denseBMat = (DenseDoubleMatrix) bMat;
            for (int row = 0; row < numRowsC; row++) {
                for (int col = 0; col < numColsC; col++) {
                    if (transposeB) {
                        cMat.set(row, col, this.cols[row].dot(denseBMat.matrix[col]));
                    } else {
                        cMat.set(row, col, this.cols[row].dot(denseBMat.matrix, col));
                    }
                }
            }
        } else if (bMat instanceof SparseColDoubleMatrix && !transposeB) {
            SparseColDoubleMatrix sparseBMat = (SparseColDoubleMatrix) bMat;
            for (int row = 0; row < numRowsC; row++) {
                for (int col = 0; col < numColsC; col++) {
                    cMat.set(row, col, this.cols[row].dot(sparseBMat.cols[col]));
                }
            }
        } else if (bMat instanceof SparseRowDoubleMatrix && transposeB) {
            SparseRowDoubleMatrix sparseBMat = (SparseRowDoubleMatrix) bMat;
            for (int row = 0; row < numRowsC; row++) {
                for (int col = 0; col < numColsC; col++) {
                    cMat.set(row, col, this.cols[row].dot(sparseBMat.rows[col]));
                }
            }
        } else {
            throw new IllegalArgumentException("unhandled type: " + bMat.getClass().getCanonicalName());
        }
    }
1262849
public void householderSymmetric( int row , double gamma )
    {
        int startU = row*N;

        // compute v = -gamma*A*u
        for( int i = row+1; i < N; i++ ) {
            double totalReal = 0;
            double totalImag = 0;

            // the lower triangle is not written to so it needs to traverse upwards
            // to get the information.  Reduces the number of matrix writes need
            // improving large matrix performance
            for( int j = row+1; j < i; j++ ) {
                double realA = QT.data[(j*N+i)*2];
                double imagA = -QT.data[(j*N+i)*2+1];

                double realU = QT.data[(startU+j)*2];
                double imagU = QT.data[(startU+j)*2+1];

                totalReal += realA*realU - imagA*imagU;
                totalImag += realA*imagU + imagA*realU;
            }
            for( int j = i; j < N; j++ ) {
                double realA = QT.data[(i*N+j)*2];
                double imagA = QT.data[(i*N+j)*2+1];

                double realU = QT.data[(startU+j)*2];
                double imagU = QT.data[(startU+j)*2+1];

                totalReal += realA*realU - imagA*imagU;
                totalImag += realA*imagU + imagA*realU;
            }
            w[i*2]   = -gamma*totalReal;
            w[i*2+1] = -gamma*totalImag;
        }

        // alpha = -0.5*gamma*u^T*v
        double realAplha = 0;
        double imageAlpha = 0;

        for( int i = row+1; i < N; i++ ) {
            double realU = QT.data[(startU+i)*2];
            double imagU = -QT.data[(startU+i)*2+1];

            double realV = w[i*2];
            double imagV = w[i*2+1];

            realAplha += realU*realV - imagU*imagV;
            imageAlpha += realU*imagV + imagU*realV;
        }
        realAplha *= -0.5*gamma;
        imageAlpha *= -0.5*gamma;

        // w = v + alpha*u
        for( int i = row+1; i < N; i++ ) {
            double realU = QT.data[(startU+i)*2];
            double imagU = QT.data[(startU+i)*2+1];

            w[i*2]   += realAplha*realU - imageAlpha*imagU;
            w[i*2+1] += realAplha*imagU + imageAlpha*realU;
        }

        // A = A + w*u^T + u*w^T
        for( int i = row+1; i < N; i++ ) {

            double realWW = w[i*2];
            double imagWW = w[i*2+1];

            double realUU = QT.data[(startU+i)*2];
            double imagUU = QT.data[(startU+i)*2+1];

            int indA = (i*N+i)*2;
            for( int j = i; j < N; j++ ) {
                // only write to the upper portion of the matrix
                // this reduces the number of cache misses

                double realU = QT.data[(startU+j)*2];
                double imagU = -QT.data[(startU+j)*2+1];

                double realW = w[j*2];
                double imagW = -w[j*2+1];

                QT.data[indA++] += realWW*realU - imagWW*imagU + realW*realUU - imagW*imagUU;
                QT.data[indA++] += realWW*imagU + imagWW*realU + realW*imagUU + imagW*realUU;
            }
        }
    }
581972
public void solve(DenseMatrix B) {
        int m = lu.nrows();
        int n = lu.ncols();
        int nrhs = B.ncols();

        if (B.nrows() != m)
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but B is %d x %d", lu.nrows(), lu.ncols(), B.nrows(), B.ncols()));

        if (isSingular()) {
            throw new RuntimeException("Matrix is singular.");
        }

        DenseMatrix X = Matrix.zeros(B.nrows(), B.ncols());

        // Copy right hand side with pivoting
        for (int j = 0; j < nrhs; j++) {
            for (int i = 0; i < m; i++) {
                X.set(i, j, B.get(piv[i], j));
            }
        }

        // Solve L*Y = B(piv,:)
        for (int k = 0; k < n; k++) {
            for (int i = k + 1; i < n; i++) {
                for (int j = 0; j < nrhs; j++) {
                    X.sub(i, j, X.get(k, j) * lu.get(i, k));
                }
            }
        }

        // Solve U*X = Y;
        for (int k = n - 1; k >= 0; k--) {
            for (int j = 0; j < nrhs; j++) {
                X.div(k, j, lu.get(k, k));
            }

            for (int i = 0; i < k; i++) {
                for (int j = 0; j < nrhs; j++) {
                    X.sub(i, j, X.get(k, j) * lu.get(i, k));
                }
            }
        }

        // Copy the result back to B.
        for (int j = 0; j < nrhs; j++) {
            for (int i = 0; i < m; i++) {
                B.set(i, j, X.get(i, j));
            }
        }
    }