math4610

Routine Name: RowEchelonForm

Author: David Merkley

Language: Python

Description/Purpose: Computes the row echelon form of a matrix

Input:

Required:

A - a matrix

Optional:

b - a solution vector

mutation - whether to modify A

leaveUpperTriangularImplicit - True means zero out lower triangle

scalePartialPivoting - Whether or not to apply scale partial pivoting

Output: Square diagonally dominant matrix

Implementation/Code:

def rowEchelonForm(A, b=[], mutation=False, leaveUpperTriangularImplicit=False, scalePartialPivoting=False):
    m = len(A)
    n = len(A[0])

    if (mutation):
        A2 = A
    else:
        A2 = [[A[i][j] for j in range(n)] for i in range(m)]

    if (scalePartialPivoting):
        s = [max(row) for row in A2]
        indexMap = [x for x in range(m)]
        for k in range(m - 1):
            maxk = k
            maxScale = A2[indexMap[k]][k] / s[indexMap[k]]
            for i in range(k + 1, m):
                iScale = A2[indexMap[i]][k] / s[indexMap[i]]
                if (iScale > maxScale):
                    maxk = i
                    maxScale = iScale
            indexMap[k], indexMap[maxk] = indexMap[maxk], indexMap[k]
    else:
        indexMap = range(m)

    if (len(b) != 0):
        for i in range(m):
            A2[i].append(b[i])
        n += 1

    for k in range(m - 1):
        krow = indexMap[k]
        for i in range(k + 1, m):
            irow = indexMap[i]
            factor = A2[irow][k] / A2[krow][k]
            for j in range(k + 1, n):
                A2[irow][j] -= factor * A2[krow][j]

    if (not leaveUpperTriangularImplicit):
        for i in range(m):
            imap = indexMap[i]
            for j in range(i):
                A2[imap][j] = 0
    return [A2[indexMap[i]] for i in range(m)]