Numeric Cholesky factorization LL=PAP'.

@param A
column-compressed matrix, only upper triangular part is used
@param S
symbolic Cholesky analysis, pinv is optional
@return numeric Cholesky factorization, null on error
164237
public static Scsn cs_chol(Scs A, Scss S) {
        float d, lki, Lx[], x[], Cx[];
        int top, i, p, k, n, Li[], Lp[], cp[], pinv[], s[], c[], parent[], Cp[], Ci[];
        Scs L, C;
        Scsn N;
        if (!Scs_util.CS_CSC(A) || S == null || S.cp == null || S.parent == null)
            return (null);
        n = A.n;
        N = new Scsn(); /* allocate result */
        c = new int[2 * n]; /* get int workspace */
        x = new float[n]; /* get float workspace */
        cp = S.cp;
        pinv = S.pinv;
        parent = S.parent;
        C = pinv != null ? Scs_symperm.cs_symperm(A, pinv, true) : A;
        s = c;
        int s_offset = n;
        Cp = C.p;
        Ci = C.i;
        Cx = C.x;
        N.L = L = Scs_util.cs_spalloc(n, n, cp[n], true, false); /* allocate result */
        Lp = L.p;
        Li = L.i;
        Lx = L.x;
        for (k = 0; k < n; k++)
            Lp[k] = c[k] = cp[k];
        for (k = 0; k < n; k++) /* compute L(k,:) for L*L' = C */
        {
            /* --- Nonzero pattern of L(k,:) ------------------------------------ */
            top = Scs_ereach.cs_ereach(C, k, parent, s, s_offset, c); /* find pattern of L(k,:) */
            x[k] = 0; /* x (0:k) is now zero */
            for (p = Cp[k]; p < Cp[k + 1]; p++) /* x = full(triu(C(:,k))) */
            {
                if (Ci[p] <= k)
                    x[Ci[p]] = Cx[p];
            }
            d = x[k]; /* d = C(k,k) */
            x[k] = 0; /* clear x for k+1st iteration */
            /* --- Triangular solve --------------------------------------------- */
            for (; top < n; top++) /* solve L(0:k-1,0:k-1) * x = C(:,k) */
            {
                i = s[s_offset + top]; /* s [top..n-1] is pattern of L(k,:) */
                lki = x[i] / Lx[Lp[i]]; /* L(k,i) = x (i) / L(i,i) */
                x[i] = 0; /* clear x for k+1st iteration */
                for (p = Lp[i] + 1; p < c[i]; p++) {
                    x[Li[p]] -= Lx[p] * lki;
                }
                d -= lki * lki; /* d = d - L(k,i)*L(k,i) */
                p = c[i]++;
                Li[p] = k; /* store L(k,i) in column i */
                Lx[p] = lki;
            }
            /* --- Compute L(k,k) ----------------------------------------------- */
            if (d <= 0)
                return null; /* not pos def */
            p = c[k]++;
            Li[p] = k; /* store L(k,k) = sqrt (d) in column k */
            Lx[p] = (float)Math.sqrt(d);
        }
        Lp[n] = cp[n]; /* finalize L */
        return N;
    }
799143
public static boolean isRowsLinearIndependent( DMatrixRMaj A )
    {
        // LU decomposition
        LUDecomposition<DMatrixRMaj> lu = DecompositionFactory_DDRM.lu(A.numRows,A.numCols);
        if( lu.inputModified() )
            A = A.copy();

        if( !lu.decompose(A))
            throw new RuntimeException("Decompositon failed?");

        // if they are linearly independent it should not be singular
        return !lu.isSingular();
    }
929254
public static FloatMatrix SVDValues(FloatMatrix A) {
        int m = A.rows;
        int n = A.columns;
        FloatMatrix S = new FloatMatrix(min(m, n));

        int info = NativeBlas.sgesvd('N', 'N', m, n, A.dup().data, 0, m, S.data, 0, null, 0, 1, null, 0, 1);

        if (info > 0) {
          throw new LapackConvergenceException("GESVD", info + " superdiagonals of an intermediate bidiagonal form failed to converge.");
        }

        return S;
    }
818862
public SymmDenseEVD factor(UpperSymmDenseMatrix A)
            throws NotConvergedException {
        if (uplo != UpLo.Upper)
            throw new IllegalArgumentException(
                    "Eigenvalue computer configured for upper-symmetrical matrices");

        return factor(A, A.getData());
    }
102747
public static DMatrixSparseCSC extractRows(DMatrixSparseCSC A , int row0 , int row1 , DMatrixSparseCSC out ) {

        if( out == null )
            out = new DMatrixSparseCSC(1,1,1);

        out.reshape(row1-row0,A.numCols,A.nz_length);
//        out.col_idx[0] = 0;
//        out.nz_length = 0;

        for (int col = 0; col < A.numCols; col++) {
            int idx0 = A.col_idx[col];
            int idx1 = A.col_idx[col+1];

            for (int i = idx0; i < idx1; i++) {
                int row = A.nz_rows[i];
                if( row >= row0 && row < row1 ) {
                    out.nz_values[out.nz_length] = A.nz_values[i];
                    out.nz_rows[out.nz_length++] = row-row0;
                }
            }
            out.col_idx[col+1] = out.nz_length;
        }

        return out;
    }
1247897
public static void symmLowerToFull( DMatrixRMaj A )
    {
        if( A.numRows != A.numCols )
            throw new MatrixDimensionException("Must be a square matrix");

        final int cols = A.numCols;

        for (int row = 0; row < A.numRows; row++) {
            for (int col = row+1; col < cols; col++) {
                A.data[row*cols+col] = A.data[col*cols+row];
            }
        }
    }
948801
@Override
    protected void sscal(long N, double a, INDArray X, int incx) {
        cblas_sscal((int) N, (float) a, (FloatPointer) X.data().addressPointer(), incx);
    }
364301
public static Scss cs_sqr(int order, Scs A, boolean qr) {
        int n, k, post[];
        Scss S;
        boolean ok = true;
        if (!Scs_util.CS_CSC(A))
            return (null); /* check inputs */
        n = A.n;
        S = new Scss(); /* allocate result S */
        S.q = Scs_amd.cs_amd(order, A); /* fill-reducing ordering */
        if (order > 0 && S.q == null)
            return (null);
        if (qr) /* QR symbolic analysis */
        {
            Scs C = order > 0 ? Scs_permute.cs_permute(A, null, S.q, false) : A;
            S.parent = Scs_etree.cs_etree(C, true); /* etree of C'*C, where C=A(:,q) */
            post = Scs_post.cs_post(S.parent, n);
            S.cp = Scs_counts.cs_counts(C, S.parent, post, true); /* col counts chol(C'*C) */
            ok = C != null && S.parent != null && S.cp != null && cs_vcount(C, S);
            if (ok)
                for (S.unz = 0, k = 0; k < n; k++)
                    S.unz += S.cp[k];
            ok = ok && S.lnz >= 0 && S.unz >= 0; /* int overflow guard */
        } else {
            S.unz = 4 * (A.p[n]) + n; /* for LU factorization only, */
            S.lnz = S.unz; /* guess nnz(L) and nnz(U) */
        }
        return (ok ? S : null); /* return result S */
    }
750145
static Integer denominator(final String s) {
        if ("".equals(s)) {
            return null;
        }
        String[] tokens = s.split("/");
        return (tokens.length < 2) ? null : emptyToNullInteger(tokens[1]);
    }
189945
def _cho_factor(A, lower=True, check_finite=True):
    """"""

    return cp.linalg.cholesky(A), True
314748
public static BandCholesky factorize(LowerSPDBandMatrix A) {
        return new BandCholesky(A.numRows(), A.kl, false).factor(A);
    }
323784
public static LinearSolverDense<DMatrixRMaj> symmPosDef(int matrixWidth ) {
        if(matrixWidth < EjmlParameters.SWITCH_BLOCK64_CHOLESKY )  {
            CholeskyDecompositionCommon_DDRM decomp = new CholeskyDecompositionInner_DDRM(true);
            return new LinearSolverChol_DDRM(decomp);
        } else {
            if( EjmlParameters.MEMORY == EjmlParameters.MemoryUsage.FASTER )
                return new LinearSolverChol_DDRB();
            else {
                CholeskyDecompositionCommon_DDRM decomp = new CholeskyDecompositionInner_DDRM(true);
                return new LinearSolverChol_DDRM(decomp);
            }
        }
    }
726928
public static boolean cholL( DMatrix6x6 A ) {

        A.a11 = Math.sqrt(A.a11);
        A.a12 = 0;
        A.a13 = 0;
        A.a14 = 0;
        A.a15 = 0;
        A.a16 = 0;
        A.a21 = (A.a21)/A.a11;
        A.a22 = Math.sqrt(A.a22-A.a21*A.a21);
        A.a23 = 0;
        A.a24 = 0;
        A.a25 = 0;
        A.a26 = 0;
        A.a31 = (A.a31)/A.a11;
        A.a32 = (A.a32-A.a31*A.a21)/A.a22;
        A.a33 = Math.sqrt(A.a33-A.a31*A.a31-A.a32*A.a32);
        A.a34 = 0;
        A.a35 = 0;
        A.a36 = 0;
        A.a41 = (A.a41)/A.a11;
        A.a42 = (A.a42-A.a41*A.a21)/A.a22;
        A.a43 = (A.a43-A.a41*A.a31-A.a42*A.a32)/A.a33;
        A.a44 = Math.sqrt(A.a44-A.a41*A.a41-A.a42*A.a42-A.a43*A.a43);
        A.a45 = 0;
        A.a46 = 0;
        A.a51 = (A.a51)/A.a11;
        A.a52 = (A.a52-A.a51*A.a21)/A.a22;
        A.a53 = (A.a53-A.a51*A.a31-A.a52*A.a32)/A.a33;
        A.a54 = (A.a54-A.a51*A.a41-A.a52*A.a42-A.a53*A.a43)/A.a44;
        A.a55 = Math.sqrt(A.a55-A.a51*A.a51-A.a52*A.a52-A.a53*A.a53-A.a54*A.a54);
        A.a56 = 0;
        A.a61 = (A.a61)/A.a11;
        A.a62 = (A.a62-A.a61*A.a21)/A.a22;
        A.a63 = (A.a63-A.a61*A.a31-A.a62*A.a32)/A.a33;
        A.a64 = (A.a64-A.a61*A.a41-A.a62*A.a42-A.a63*A.a43)/A.a44;
        A.a65 = (A.a65-A.a61*A.a51-A.a62*A.a52-A.a63*A.a53-A.a64*A.a54)/A.a55;
        A.a66 = Math.sqrt(A.a66-A.a61*A.a61-A.a62*A.a62-A.a63*A.a63-A.a64*A.a64-A.a65*A.a65);
        return !UtilEjml.isUncountable(A.a66);
    }
512370
public static double elementMinAbs( DMatrixSparseCSC A ) {
        if( A.nz_length == 0)
            return 0;

        double min = A.isFull() ? Math.abs(A.nz_values[0]) : 0;
        for(int i = 0; i < A.nz_length; i++ ) {
            double val = Math.abs(A.nz_values[i]);
            if( val < min ) {
                min = val;
            }
        }

        return min;
    }
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];
                }
            }
        }

    }
1079806
def cholesky(X):
    

    assert isinstance(X, cspmatrix) and X.is_factor is False, "X must be a cspmatrix"

    n = X.symb.n
    snpost = X.symb.snpost
    snptr = X.symb.snptr
    chptr = X.symb.chptr
    chidx = X.symb.chidx

    relptr = X.symb.relptr
    relidx = X.symb.relidx
    blkptr = X.symb.blkptr
    blkval = X.blkval

    stack = []

    for k in snpost:

        nn = snptr[k+1]-snptr[k]       # |Nk|
        na = relptr[k+1]-relptr[k]     # |Ak|
        nj = na + nn                   

        # build frontal matrix
        F = matrix(0.0, (nj, nj))
        lapack.lacpy(blkval, F, offsetA = blkptr[k], m = nj, n = nn, ldA = nj, uplo = 'L')

        # add update matrices from children to frontal matrix
        for i in range(chptr[k+1]-1,chptr[k]-1,-1):
            Ui = stack.pop()
            frontal_add_update(F, Ui, relidx, relptr, chidx[i])

        # factor L_{Nk,Nk}
        lapack.potrf(F, n = nn, ldA = nj)

        # if supernode k is not a root node, compute and push update matrix onto stack
        if na > 0:   
            # compute L_{Ak,Nk} := A_{Ak,Nk}*inv(L_{Nk,Nk}')
            blas.trsm(F, F, m = na, n = nn, ldA = nj, 
                      ldB = nj, offsetB = nn, transA = 'T', side = 'R')

            # compute Uk = Uk - L_{Ak,Nk}*inv(D_{Nk,Nk})*L_{Ak,Nk}'
            if nn == 1:
                blas.syr(F, F, n = na, offsetx = nn, \
                         offsetA = nn*nj+nn, ldA = nj, alpha = -1.0)
            else:
                blas.syrk(F, F, k = nn, n = na, offsetA = nn, ldA = nj,
                          offsetC = nn*nj+nn, ldC = nj, alpha = -1.0, beta = 1.0)

            # compute L_{Ak,Nk} := L_{Ak,Nk}*inv(L_{Nk,Nk})
            blas.trsm(F, F, m = na, n = nn,\
                      ldA = nj, ldB = nj, offsetB = nn, side = 'R')

            # add Uk to stack
            Uk = matrix(0.0,(na,na))
            lapack.lacpy(F, Uk, m = na, n = na, uplo = 'L', offsetA = nn*nj+nn, ldA = nj)
            stack.append(Uk)

        # copy the leading Nk columns of frontal matrix to blkval
        lapack.lacpy(F, blkval, uplo = "L", offsetB = blkptr[k], m = nj, n = nn, ldB = nj)        

    X.is_factor = True

    return
984780
@Override
    public void sgesvd(byte jobu, byte jobvt, int M, int N, INDArray A, INDArray S, INDArray U, INDArray VT,
                    INDArray INFO) {
        INDArray superb = Nd4j.create(DataType.FLOAT, M < N ? M : N ) ;
        int status = LAPACKE_sgesvd(getColumnOrder(A), jobu, jobvt, M, N, 
                        (FloatPointer)A.data().addressPointer(), getLda(A),
                        (FloatPointer)S.data().addressPointer(), 
                        U == null ? null : (FloatPointer)U.data().addressPointer(), U == null ? 1 : getLda(U),
                        VT == null ? null : (FloatPointer)VT.data().addressPointer(), VT == null ? 1 : getLda(VT), 
                        (FloatPointer)superb.data().addressPointer() 
                        );
        if( status != 0 ) {
            throw new BlasException( "Failed to execute sgesvd", status ) ;
        }
    }