0
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.

Nick Brady
  • 107
  • 8

1 Answers1

0

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.

Raghunathan Ramakrishnan
  • 141
  • 4
  • Yes, it seems to be looking to exit gracefully when `DETERM = 0`. But even when the determinant is not zero, it still chooses the pivot rows in a way that is not intuitive to me. I would just go from the top to the bottom and swap rows only when there is a zero at the pivot point. But the process of this code looks for the minimum value of `BMAX`. Mostly I am curious about why it does this: what are the advantages and the reasoning behind it. What special types of `B` matrices do you think it is tailored to? In my use cases, `B` is less than 20x20, but can be densely filled. – Nick Brady Feb 22 '22 at 11:26
  • I found the original reference where this method was first introduced. The algorithm was proposed in Appendix C of the book "Electrochemical Systems" by John Newman and Karen Thomas-Alyea. They have some descriptions of the problem with equations. – Raghunathan Ramakrishnan Feb 22 '22 at 11:46
  • That is correct, but the only description given in that text is "The subroutine `MATINV` is used to solve the linear equations C.19, C.23, C.24, and C.27, which arise at each value of j . If, at any value of j, the determinant of the matrix of these equations is found to be zero, this fact is reported in the output." There is no further information about `MATINV` - it does not even indicate that the method of inversion is Gaussian-elimination. I have also looked at the reference sources in Appendix C, but they are similarly sparse in their descriptions of `MATINV`. – Nick Brady Feb 22 '22 at 12:56