Improve and expand applicability of RREF algorithms

[dsvandet]

List of QEC requirements being added
List the names and id's of the QEC requirements being addressed. If there is no requirement then propose a new requirement first.

Describe the solution you'd proposing

The rref implemented so far, in linear.utils.py is not parallel and is not sophisticated. Will work just fine for relatively small matrices but as matrices grow, as in LDPC code check matrices, we should have a much better implementation. Specifically we should implement a general rref method that calls the appropriate methods depending on the characteristics of the matrix. I would start by implementing the algorithms in the following paper: arXiv:1111.6549 from 2011 for matrices over GF(2). For matrices over GF(q) and alike the approach below can be easily modified as a start.

The current rref implementation is as follows:

def _rref_complete(matrix):
    """Compute the Row Reduced Echelon Form for a boolean matrix and associated 
    results.

    Args:
        matrix (2d matrix): numpy.ndarry of a boolean matrix

    Returns:
        heads (list): A 0-1 list with value of position k being 1 if the kth column is a pivot
        rref_mat (numpy.ndarray) : Row Reduced Echelon Form of input matrix
        transform_mat : Transform used to tranform input matrix into RREF form
        rank :  rank of input matrix
    """
    matrix = augment_mat(matrix, "right")

    nrows = matrix.shape[0]
    ncols = matrix.shape[1]
    hncols = ncols-nrows

    co_row = []
    null_space = []
    heads = hncols*[-1]
    
    # First clear entries below pivot
    for i in range(nrows):
        row = matrix[i]

        for j in range(hncols):
            if heads[j] != -1:
                if row[j] == True:
                    row = row ^ co_row[heads[j]]

        # Find the first non-zero value
        j = np.nonzero(row)[0][0]
        if j < hncols:
            co_row.append(row)
            heads[j] = len(co_row)-1

    # Clear entries above pivots
    for j in range(hncols-1,-1,-1):
        if heads[j] != -1:
            diff = [item for item in range(heads[j]) if item not in heads[j+1:]]
            for i in diff:
                row = co_row[i]
                if row[j] == True:
                    co_row[i] = co_row[i] ^ co_row[heads[j]]
                    
    non_zero_heads = [i for i in heads if i != -1]
    co_row = [co_row[i] for i in non_zero_heads]
    co_row = np.array(co_row, dtype=bool)
    rank = sum(non_zero_heads)
    
    rref_mat = co_row[:,:hncols]
    transform_mat = co_row[:,hncols:]
    heads = [int(bool(x+1)) for x in heads]
    
    return heads, rref_mat, transform_mat, rank