public static Dcsn cs_lu(Dcs A, Dcss S, double tol) {
        Dcs L, U;
        Dcsn N;
        double pivot, Lx[], Ux[], x[], a, t;
        int Lp[], Li[], Up[], Ui[], pinv[], xi[], q[], n, ipiv, k, top, p, i, col, lnz, unz;
        if (!Dcs_util.CS_CSC(A) || S == null)
            return (null); /* check inputs */
        n = A.n;
        q = S.q;
        lnz = S.lnz;
        unz = S.unz;
        x = new double[n]; /* get double workspace */
        xi = new int[2 * n]; /* get int workspace */
        N = new Dcsn(); /* allocate result */
        N.L = L = Dcs_util.cs_spalloc(n, n, lnz, true, false); /* allocate result L */
        N.U = U = Dcs_util.cs_spalloc(n, n, unz, true, false); /* allocate result U */
        N.pinv = pinv = new int[n]; /* allocate result pinv */
        Lp = L.p;
        Up = U.p;
        for (i = 0; i < n; i++)
            x[i] = 0; /* clear workspace */
        for (i = 0; i < n; i++)
            pinv[i] = -1; /* no rows pivotal yet */
        for (k = 0; k <= n; k++)
            Lp[k] = 0; /* no cols of L yet */
        lnz = unz = 0;
        for (k = 0; k < n; k++) /* compute L(:,k) and U(:,k) */
        {
            /* --- Triangular solve --------------------------------------------- */
            Lp[k] = lnz; /* L(:,k) starts here */
            Up[k] = unz; /* U(:,k) starts here */
            if (lnz + n > L.nzmax) {
                Dcs_util.cs_sprealloc(L, 2 * L.nzmax + n);
            }
            if (unz + n > U.nzmax) {
                Dcs_util.cs_sprealloc(U, 2 * U.nzmax + n);
            }
            Li = L.i;
            Lx = L.x;
            Ui = U.i;
            Ux = U.x;
            col = q != null ? (q[k]) : k;
            top = Dcs_spsolve.cs_spsolve(L, A, col, xi, x, pinv, true); /* x = L\A(:,col) */
            /* --- Find pivot --------------------------------------------------- */
            ipiv = -1;
            a = -1;
            for (p = top; p < n; p++) {
                i = xi[p]; /* x(i) is nonzero */
                if (pinv[i] < 0) /* row i is not yet pivotal */
                {
                    if ((t = Math.abs(x[i])) > a) {
                        a = t; /* largest pivot candidate so far */
                        ipiv = i;
                    }
                } else /* x(i) is the entry U(pinv[i],k) */
                {
                    Ui[unz] = pinv[i];
                    Ux[unz++] = x[i];
                }
            }
            if (ipiv == -1 || a <= 0)
                return (null);
            if (pinv[col] < 0 && Math.abs(x[col]) >= a * tol)
                ipiv = col;
            /* --- Divide by pivot ---------------------------------------------- */
            pivot = x[ipiv]; /* the chosen pivot */
            Ui[unz] = k; /* last entry in U(:,k) is U(k,k) */
            Ux[unz++] = pivot;
            pinv[ipiv] = k; /* ipiv is the kth pivot row */
            Li[lnz] = ipiv; /* first entry in L(:,k) is L(k,k) = 1 */
            Lx[lnz++] = 1;
            for (p = top; p < n; p++) /* L(k+1:n,k) = x / pivot */
            {
                i = xi[p];
                if (pinv[i] < 0) /* x(i) is an entry in L(:,k) */
                {
                    Li[lnz] = i; /* save unpermuted row in L */
                    Lx[lnz++] = x[i] / pivot; /* scale pivot column */
                }
                x[i] = 0; /* x [0..n-1] = 0 for next k */
            }
        }
        /* --- Finalize L and U ------------------------------------------------- */
        Lp[n] = lnz;
        Up[n] = unz;
        Li = L.i; /* fix row indices of L for final pinv */
        for (p = 0; p < lnz; p++)
            Li[p] = pinv[Li[p]];
        Dcs_util.cs_sprealloc(L, 0); /* remove extra space from L and U */
        Dcs_util.cs_sprealloc(U, 0);
        return N;
    }
public static Dcsn cs_lu ( Dcs A , Dcss S , double tol ) { Dcs L , U ; Dcsn N ; double pivot , Lx [ ] , Ux [ ] , x [ ] , a , t ; int Lp [ ] , Li [ ] , Up [ ] , Ui [ ] , pinv [ ] , xi [ ] , q [ ] , n , ipiv , k , top , p , i , col , lnz , unz ; if ( ! Dcs_util . CS_CSC ( A ) || S == null ) return ( null ) ; n = A . n ; q = S . q ; lnz = S . lnz ; unz = S . unz ; x = new double [ n ] ; xi = new int [ 2 * n ] ; N = new Dcsn ( ) ; N . L = L = Dcs_util . cs_spalloc ( n , n , lnz , true , false ) ; N . U = U = Dcs_util . cs_spalloc ( n , n , unz , true , false ) ; N . pinv = pinv = new int [ n ] ; Lp = L . p ; Up = U . p ; for ( i = 0 ; i < n ; i ++ ) x [ i ] = 0 ; for ( i = 0 ; i < n ; i ++ ) pinv [ i ] = - 1 ; for ( k = 0 ; k <= n ; k ++ ) Lp [ k ] = 0 ; lnz = unz = 0 ; for ( k = 0 ; k < n ; k ++ ) { Lp [ k ] = lnz ; Up [ k ] = unz ; if ( lnz + n > L . nzmax ) { Dcs_util . cs_sprealloc ( L , 2 * L . nzmax + n ) ; } if ( unz + n > U . nzmax ) { Dcs_util . cs_sprealloc ( U , 2 * U . nzmax + n ) ; } Li = L . i ; Lx = L . x ; Ui = U . i ; Ux = U . x ; col = q != null ? ( q [ k ] ) : k ; top = Dcs_spsolve . cs_spsolve ( L , A , col , xi , x , pinv , true ) ; ipiv = - 1 ; a = - 1 ; for ( p = top ; p < n ; p ++ ) { i = xi [ p ] ; if ( pinv [ i ] < 0 ) { if ( ( t = Math . abs ( x [ i ] ) ) > a ) { a = t ; ipiv = i ; } } else { Ui [ unz ] = pinv [ i ] ; Ux [ unz ++ ] = x [ i ] ; } } if ( ipiv == - 1 || a <= 0 ) return ( null ) ; if ( pinv [ col ] < 0 && Math . abs ( x [ col ] ) >= a * tol ) ipiv = col ; pivot = x [ ipiv ] ; Ui [ unz ] = k ; Ux [ unz ++ ] = pivot ; pinv [ ipiv ] = k ; Li [ lnz ] = ipiv ; Lx [ lnz ++ ] = 1 ; for ( p = top ; p < n ; p ++ ) { i = xi [ p ] ; if ( pinv [ i ] < 0 ) { Li [ lnz ] = i ; Lx [ lnz ++ ] = x [ i ] / pivot ; } x [ i ] = 0 ; } } Lp [ n ] = lnz ; Up [ n ] = unz ; Li = L . i ; for ( p = 0 ; p < lnz ; p ++ ) Li [ p ] = pinv [ Li [ p ] ] ; Dcs_util . cs_sprealloc ( L , 0 ) ; Dcs_util . cs_sprealloc ( U , 0 ) ; return N ; }
Sparse LU factorization of a square matrix, PAQ = LU.

@param A
column-compressed matrix
@param S
symbolic LU analysis
@param tol
partial pivoting threshold (1 for partial pivoting)
@return numeric LU factorization, null on error