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 )


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)





GJ_Inverse_alt( A, B )

print( A )
print( B )

print( Aold * B )