In this notebook, you will implement the alternative Gauss Jordan algorithm that overwrites $ A $ in one sweep with the identity matrix and $ B $ with the inverse of the original matrix $ A $.
Be sure to make a copy!!!!
import numpy as np
import laff
import flame
L = np.matrix( ' 1, 0, 0. 0;\
-2, 1, 0, 0;\
1,-3, 1, 0;\
2, 3,-1, 1' )
U = np.matrix( ' 2,-1, 3,-2;\
0,-2, 1,-1;\
0, 0, 1, 2;\
0, 0, 0, 3' )
A = L * U
Aold = np.matrix( np.copy( A ) )
B = np.matrix( np.eye( 4 ) )
print( 'A = ' )
print( A )
print( 'B = ' )
print( B )
Here is the algorithm:
Important: if you make a mistake, rerun ALL cells above the cell in which you were working, and then the one where you are working.
Create the routine
GJ_Inverse_alt
with the Spark webpage for the algorithm
import flame
def GJ_Inverse_alt(A, B):
ATL, ATR, \
ABL, ABR = flame.part_2x2(A, \
0, 0, 'TL')
BTL, BTR, \
BBL, BBR = flame.part_2x2(B, \
0, 0, 'TL')
while ATL.shape[0] < A.shape[0]:
A00, a01, A02, \
a10t, alpha11, a12t, \
A20, a21, A22 = flame.repart_2x2_to_3x3(ATL, ATR, \
ABL, ABR, \
1, 1, 'BR')
B00, b01, B02, \
b10t, beta11, b12t, \
B20, b21, B22 = flame.repart_2x2_to_3x3(BTL, BTR, \
BBL, BBR, \
1, 1, 'BR')
#------------------------------------------------------------#
# a01 := a01 / alpha11
# a21 := a21 / alpha11
laff.invscal( alpha11, a01 )
laff.invscal( alpha11, a21 )
# A02 := A02 - a01 * a12t
laff.ger( -1.0, a01, a12t, A02 )
# A22 := A22 - a21 * a12t
laff.ger( -1.0, a21, a12t, A22 )
# B00 := B00 - a01 * b01t
laff.ger( -1.0, a01, b10t, B00 )
# B20 := B20 - a21 * b01t
laff.ger( -1.0, a21, b10t, B20 )
# b01 := - a01 (= - u01 in the discussion)
laff.copy( a01, b01 )
laff.scal( -1.0, b01 )
# b21 := - a21 (= - l21 in the discussion)
laff.copy( a21, b21 )
laff.scal( -1.0, b21 )
# a12t:= a21t / alpha11
laff.invscal( alpha11, a12t )
# b10t:= b10t / alpha11
laff.invscal( alpha11, b10t )
# beta11 := 1.0 / alpha11
laff.invscal( alpha11, beta11 )
# a01 = 0 (zero vector)
# alpha11 = 1
# a21 = 0 (zero vector)
laff.zerov( a01 )
laff.onev( alpha11 )
laff.zerov( a21 )
#------------------------------------------------------------#
ATL, ATR, \
ABL, ABR = flame.cont_with_3x3_to_2x2(A00, a01, A02, \
a10t, alpha11, a12t, \
A20, a21, A22, \
'TL')
BTL, BTR, \
BBL, BBR = flame.cont_with_3x3_to_2x2(B00, b01, B02, \
b10t, beta11, b12t, \
B20, b21, B22, \
'TL')
flame.merge_2x2(ATL, ATR, \
ABL, ABR, A)
flame.merge_2x2(BTL, BTR, \
BBL, BBR, B)
Important: if you make a mistake, rerun ALL cells above the cell in which you were working, and then the one where you are working.
GJ_Inverse_alt( A, B )
print( A )
print( B )
Matrix $ A $ should now be an identity matrix and $ B $ should no longer be an identity matrix.
Check if $ B $ now equals (approximately) the inverse of the original matrix $ A $:
print( Aold * B )