import numpy as np

m = 4
n = 3
k = 5

C = np.matrix( np.random.random( (m, n) ) )
print( 'C = ' )
print( C )

Cold = np.matrix( np.zeros( (m,n ) ) )
Cold = np.matrix( np.copy( C ) )   # an alternative way of doing a "hard" copy, in this case of a matrix
    
A = np.matrix( np.random.random( (m, k) ) )
print( 'A = ' )
print( A )

B = np.matrix( np.random.random( (k, n) ) )
print( 'B = ' )
print( B )

def MMmult_lots_of_loops( A, B, C ):

    m, n = np.shape( C )
    m, k = np.shape( A )
    
    # i,j,p
    for i in range( m ):                     
        for j in range( n ):                    
            for p in range( k ):                    
                C[ i,j ] = A[ i,p ] * B[ p, j ] + C[ i,j ]
                
    # i,p,j
#    for i in range( m ):                     
#        for p in range( k ):                    
#            for j in range( n ):                    
#                C[ i,j ] = A[ i,p ] * B[ p, j ] + C[ i,j ]
                
    # j,i,p
#    for j in range( n ):                     
#        for i in range( m ):                    
#            for p in range( k ):                    
#                C[ i,j ] = A[ i,p ] * B[ p, j ] + C[ i,j ]

    # j,p,i
#    for j in range( n ):                     
#        for p in range( k ):                    
#            for i in range( m ):                    
#                C[ i,j ] = A[ i,p ] * B[ p, j ] + C[ i,j ]

    # p,i,j
#    for p in range( k ):                     
#        for i in range( m ):                    
#            for j in range( n ):                    
#                C[ i,j ] = A[ i,p ] * B[ p, j ] + C[ i,j ]

    # p,j,i
#    for p in range( k ):                     
#        for j in range( n ):                    
#            for i in range( m ):                    
#                C[ i,j ] = A[ i,p ] * B[ p, j ] + C[ i,j ]

C = np.matrix( np.copy( Cold ) )             # restore C

MMmult_lots_of_loops( A, B, C )

print( 'C - ( Cold + A * B )' )
print( C - ( Cold + A * B ) )

def MMmult_C_eq_AB( A, B, C ):

    m, n = np.shape( C )
    m, k = np.shape( A )
    
    for i in range( m ):                     
        for j in range( n ):   
            C[ i,j ] = 0.0
            for p in range( k ):                    
                C[ i,j ] = A[ i,p ] * B[ p, j ] + C[ i,j ]

C = np.matrix( np.copy( Cold ) )             # restore C

MMmult_C_eq_AB( A, B, C )

print( 'C - ( A * B )' )
print( C - ( A * B ) )

C = np.matrix( np.zeros( np.shape( C ) ) )          # initialize C = 0 

MMmult_lots_of_loops( A, B, C )

print( 'C - ( A * B )' )
print( C - ( A * B ) )