Matrix QRsolve(Matrix B, Matrix X, Transpose trans) {
        int nrhs = B.numColumns();

        // Allocate temporary solution matrix
        DenseMatrix Xtmp = new DenseMatrix(Math.max(numRows, numColumns), nrhs);
        int M = trans == Transpose.NoTranspose ? numRows : numColumns;
        for (int j = 0; j < nrhs; ++j)
            for (int i = 0; i < M; ++i)
                Xtmp.set(i, j, B.get(i, j));
        double[] newData = data.clone();

        // Query optimal workspace
        double[] work = new double[1];
        intW info = new intW(0);
        LAPACK.getInstance().dgels(trans.netlib(), numRows, numColumns, nrhs,
                newData, Matrices.ld(numRows), Xtmp.getData(), Matrices.ld(numRows, numColumns),
                work, -1, info);

        // Allocate workspace
        int lwork = -1;
        if (info.val != 0)
            lwork = Math.max(1, Math.min(numRows, numColumns)
                    + Math.max(Math.min(numRows, numColumns), nrhs));
        else
            lwork = Math.max((int) work[0], 1);
        work = new double[lwork];

        // Compute the factorization
        info.val = 0;
        LAPACK.getInstance().dgels(trans.netlib(), numRows, numColumns, nrhs,
                newData, Matrices.ld(numRows), Xtmp.getData(), Matrices.ld(numRows, numColumns),
                work, lwork, info);

        if (info.val < 0)
            throw new IllegalArgumentException();

        // Extract the solution
        int N = trans == Transpose.NoTranspose ? numColumns : numRows;
        for (int j = 0; j < nrhs; ++j)
            for (int i = 0; i < N; ++i)
                X.set(i, j, Xtmp.get(i, j));
        return X;
    }