private void factor() {
        int n = LU.numRows();

        // Internal CRS matrix storage
        int[] colind = LU.getColumnIndices();
        int[] rowptr = LU.getRowPointers();
        double[] data = LU.getData();

        // Find the indices to the diagonal entries
        int[] diagind = findDiagonalIndices(n, colind, rowptr);

        // Go down along the main diagonal
        for (int k = 1; k < n; ++k)
            for (int i = rowptr[k]; i < diagind[k]; ++i) {

                // Get the current diagonal entry
                int index = colind[i];
                double LUii = data[diagind[index]];

                if (LUii == 0)
                    throw new RuntimeException("Zero pivot encountered on row "
                            + (i + 1) + " during ILU process");

                // Elimination factor
                double LUki = (data[i] /= LUii);

                // Traverse the sparse row i, reducing on row k
                for (int j = diagind[index] + 1, l = rowptr[k] + 1; j < rowptr[index + 1]; ++j) {

                    while (l < rowptr[k + 1] && colind[l] < colind[j])
                        l++;

                    if (colind[l] == colind[j])
                        data[l] -= LUki * data[j];
                }
            }

        L = new UnitLowerCompRowMatrix(LU, diagind);
        U = new UpperCompRowMatrix(LU, diagind);
    }