public AbstractMatrix[] decomposePenta() throws IllegalArgumentException {
        final int n = rows();
        if (this.isPentadiagonale() == false)
            throw new IllegalArgumentException("Matrix is not pentadiagonale. Cannot do the pentadiagonal decomposition.");

        BandMatrix L = new BandMatrix(n, new int[]{0, -1, -2});
        BandMatrix R = new BandMatrix(n, new int[]{0, 1, 2});
        double[] m = new double[n];
        double[] l = new double[n];
        double[] k = new double[n];
        double[] p = new double[n];
        double[] r = new double[n];
        // initialization
        m[0] = get(0, 0);
        r[0] = get(0, 1);
        l[1] = get(1, 0) / m[0];
        m[1] = get(1, 1) - l[1] * r[0];

        // p_i s
        for (int i = 0; i < n - 2; i++) {
            p[i] = get(i, i + 2);
        }

        for (int i = 2; i < n; i++) {
            k[i] = get(i, i - 2) / m[i - 2];
            l[i] = (get(i, i - 1) - k[i] * r[i - 2]) / m[i - 1];
            r[i - 1] = get(i - 1, i) - (l[i - 1] * (p[i - 2]));
            m[i] = get(i, i) - k[i] * p[i - 2] - l[i] * r[i - 1];
        }

        // build L
        for (int i = 0; i < n; i++) {
            L.set(i, i, 1);
            if (i >= 1)
                L.set(i, i - 1, l[i]);
            if (i >= 2)
                L.set(i, i - 2, k[i]);
        }

        //build R
        for (int i = 0; i < n; i++) {
            R.set(i, i, m[i]);
            if (i < n - 1)
                R.set(i, i + 1, r[i]);
            if (i < n - 2)
                R.set(i, i + 2, p[i]);
        }

        return new AbstractMatrix[]{L, R};
    }