What Numerical Stability conditioning is this algorithm using for Gaussian-elimination matrix inversion?

Question

import numpy as np

def MATINV(N, A):
    '''
    A = [B | D]
    returns
    A = [I | B^-1 * D]

    if A[:,:N] \= eye(N)
        then deter(B) == 0

    N - number of rows (N <= number of columns)

        BTRY    - largest value in a row
        BNEXT   - second largest value in a row
        BMAX    - ratio of the two largest magnitude values in a row (BNEXT / BTRY)
        ID      - keeps track of available pivot rows (array)

    pivot row    - row with the smallest BMAX value (excluding already used pivot rows)
    pivot column - max magnitude index in pivot row (excluding already used pivot columns)

    '''
    Determ = 0.0
    ID = np.zeros(N)

    for NN in range(N):
#########################################################################################
# Don't quite understand why this is done 
#########################################################################################
        BMAX = 1.1

        # print(A[:,:N])

        for I in range(N):                  # I - row index
            if (ID[I] == 0):
                BNEXT = 0.0
                BTRY = 0.0

                for J in range(N):          # J - column index
                    if (ID[J] == 0):
                        if not(abs(A[I,J]) <= BNEXT):       # if (abs(A[I,J]) > BNEXT)):
                            BNEXT = abs(A[I,J])
                            if not(BNEXT <= BTRY):          # if (BNEXT > BTRY):
                                BNEXT = BTRY
                                BTRY = abs(A[I,J])
                                JC = J
                    # print('({},{})'.format(I,J), JC, 'BMAX', BMAX, 'BNEXT', BNEXT, 'BTRY', BTRY)

                if not(BNEXT >= BMAX*BTRY):                 # if (BNEXT < BMAX*BTRY):
                    BMAX = BNEXT/BTRY                       # if (BNEXT / BTRY < BMAX)
                    IROW = I
                    JCOL = JC
#########################################################################################


        if (ID[JC] != 0):
            DETERM = 0.0
            return

        # print('ID', ID)
        # print("IROW", "JCOL:", IROW, JCOL)
        # print('BMAX', BMAX, 'BNEXT', BNEXT, 'BTRY', BTRY)

        ID[JCOL] = 1


        # else:
        # swap rows(JCOL, IROW)
        if JCOL != IROW:
            A_temp     = A[JCOL,:].copy()
            A[JCOL,:]  = A[IROW,:]
            A[IROW,:]  = A_temp

        # make the leading value 1
        A[JCOL, :] /= A[JCOL, JCOL]

        # subtract pivot row from all rows BELOW / ABOVE pivot row
        for ii in range(N):
            if ii == JCOL:
                continue
            f = A[ii,JCOL]
            A[ii,:] -= f * A[JCOL,:]

    return A

I'm working to convert some fortran code to python and one of the subroutines is MATINV. I understand the code at a high level:

• input: matrix A that is composed of two matrices B and D. Matrix B is N by N and matrix D has N rows, but can have more columns.
• returns B^-1 * D using Gaussian-elimination. (Matrix A is returned as A = [ I | B^-1 * D], where I is the identity matrix and the relevant portion is A[:,N:])

The part I do not quite understand is marked with hash symbols (#). I understand what the code is doing:

• it is looking for the matrix row with the smallest ratio of BNEXT / BTRY, where BNEXT is the value with the second largest magnitude, and BTRY is the value with the largest magnitude (BMAX = BNEXT / BTRY)
• The row with with min(BMAX) is selected as the pivot row

I believe that this conditioning of the pivot rows is done for numerical stability, but the specific numerical theory behind it is not clear to me. The current structure of the code seems a bit nebulous to me, and I think understanding the fundamental idea guiding this algorithm might help me write clearer code. Any insights will be helpful.

Answer 1

Seems to be tailored to a particular class of B matrices. More than preconditioning, the code block seems to be looking for instances when the determinant of B can be zero and terminate gracefully.