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

        double[] LUi = new double[n];

        // Find the indices to the diagonal entries
        for (int k = 0; k < n; ++k) {
            SparseVector row = LU.getRow(k);
            diagind[k] = findDiagonalIndex(row, k);
            if (diagind[k] < 0)
                throw new RuntimeException("Missing diagonal entry on row "
                        + (k + 1));
        }

        for (int i = 1; i < n; ++i) {

            // Get row i
            SparseVector rowi = LU.getRow(i);

            // Drop tolerance on current row
            double taui = rowi.norm(Vector.Norm.Two) * tau;

            // Store in dense format
            scatter(rowi, LUi);

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

                // Get row k
                SparseVector rowk = LU.getRow(k);
                int[] rowIndex = rowk.getIndex();
                int rowUsed = rowk.getUsed();
                double[] rowData = rowk.getData();

                if (rowData[diagind[k]] == 0)
                    throw new RuntimeException("Zero diagonal entry on row "
                            + (k + 1) + " during ILU process");

                double LUik = LUi[k] / rowData[diagind[k]];

                // Check for small elimination entry
                if (Math.abs(LUik) <= taui)
                    continue;

                // Traverse the sparse row k, reducing row i
                for (int j = diagind[k] + 1; j < rowUsed; ++j)
                    LUi[rowIndex[j]] -= LUik * rowData[j];

                // The above has overwritten LUik, so remedy that
                LUi[k] = LUik;
            }

            // Store back into the LU matrix, dropping as needed
            gather(LUi, rowi, taui, i);

            // Update diagonal index on row i if it is outdated
            if (rowi.getIndex()[diagind[i]] != i) {
                diagind[i] = findDiagonalIndex(rowi, i);
                if (diagind[i] < 0)
                    throw new RuntimeException("Missing diagonal entry on row "
                            + (i + 1) + " during ILU process");
            }
        }

        L = new UnitLowerFlexCompRowMatrix(LU, diagind);
        U = new UpperFlexCompRowMatrix(LU, diagind);
    }