This notebook helps you implement the operation $ C := U R $ where $ C, U, R \in \mathbb{R}^{n \times n} $, and $ U $ and $ R $ are upper triangular. $ U $ and $ R $ are stored in the upper triangular part of numpy matrices U
and R
. The upper triangular part of matrix C
is to be overwritten with the resulting upper triangular matrix.
First, we create some matrices.
import numpy as np
n = 5
C = np.matrix( np.random.random( (n, n) ) )
print( 'C = ' )
print( C )
Cold = np.matrix( np.zeros( (n,n ) ) )
Cold = np.matrix( np.copy( C ) ) # an alternative way of doing a "hard" copy, in this case of a matrix
U = np.matrix( np.random.random( (n, n) ) )
print( 'U = ' )
print( U )
R = np.matrix( np.random.random( (n, n) ) )
print( 'R = ' )
print( R )
Trtrmm_uu_unb_var1( U, R, C )
This routine computes $ C := U R + C $. The "_uu_" means that $ U $ and $ R $ are upper triangular matrices (which means $ C $ is too). However, the lower triangular parts of numpy matrices U
, R
, and C
are not to be "touched".
The specific laff function you will want to use is some subset of
laff.gemv( trans, alpha, A, x, beta, y )
which computes $ y := \alpha A x + \beta y $ or $ y := \alpha A^T x + \beta y $ depending on whether trans = 'No transpose'
or trans = 'Transpose'
laff.ger( alpha, x, y, A )
which computes the rank-1 update (adds a multiple of an outer product to a matrix)
$ A := \alpha x y^T + A $. laff.axpy( alpha, x, y )
laff.dots( x, y, alpha )
Hint: If you multiply with $ U_{00} $, you will want to use np.triu( U00 )
to make sure you don't compute with the nonzeroes below the diagonal.
Use the Spark webpage to generate a code skeleton. (Make sure you adjust the name of the routine.)
import laff as laff
import flame
def Trtrmm_uu_unb_var1(U, R, C):
UTL, UTR, \
UBL, UBR = flame.part_2x2(U, \
0, 0, 'TL')
RTL, RTR, \
RBL, RBR = flame.part_2x2(R, \
0, 0, 'TL')
CTL, CTR, \
CBL, CBR = flame.part_2x2(C, \
0, 0, 'TL')
while UTL.shape[0] < U.shape[0]:
U00, u01, U02, \
u10t, upsilon11, u12t, \
U20, u21, U22 = flame.repart_2x2_to_3x3(UTL, UTR, \
UBL, UBR, \
1, 1, 'BR')
R00, r01, R02, \
r10t, rho11, r12t, \
R20, r21, R22 = flame.repart_2x2_to_3x3(RTL, RTR, \
RBL, RBR, \
1, 1, 'BR')
C00, c01, C02, \
c10t, gamma11, c12t, \
C20, c21, C22 = flame.repart_2x2_to_3x3(CTL, CTR, \
CBL, CBR, \
1, 1, 'BR')
#------------------------------------------------------------#
# imporant: You do not need to recompute C00 = U00 * R00!!!!
laff.gemv( 'No transpose', 1.0, np.triu( U00 ), r01, 0.0, c01 )
laff.axpy( rho11, u01, c01 )
laff.dot( rho11, upsilon11, gamma11 )
#------------------------------------------------------------#
UTL, UTR, \
UBL, UBR = flame.cont_with_3x3_to_2x2(U00, u01, U02, \
u10t, upsilon11, u12t, \
U20, u21, U22, \
'TL')
RTL, RTR, \
RBL, RBR = flame.cont_with_3x3_to_2x2(R00, r01, R02, \
r10t, rho11, r12t, \
R20, r21, R22, \
'TL')
CTL, CTR, \
CBL, CBR = flame.cont_with_3x3_to_2x2(C00, c01, C02, \
c10t, gamma11, c12t, \
C20, c21, C22, \
'TL')
flame.merge_2x2(CTL, CTR, \
CBL, CBR, C)
C = np.matrix( np.copy( Cold ) ) # restore C
Trtrmm_uu_unb_var1( U, R, C )
# compute it using numpy *. This is a little complex, since we want to make sure we
# don't change the lower triangular part of C
# Cref = np.tril( Cold, -1 ) + np.triu( U ) * np.triu( R )
Cref = np.triu( U ) * np.triu( R ) + np.tril( Cold, -1 )
print( 'C - Cref' )
print( C - Cref )
In theory, the result matrix should be (approximately) zero.
Copy and paste the code into PictureFLAME , a webpage where you can watch your routine in action. Just cut and paste into the box.
Disclaimer: we implemented a VERY simple interpreter. If you do something wrong, we cannot guarantee the results. But if you do it right, you are in for a treat.
If you want to reset the problem, just click in the box into which you pasted the code and hit "next" again.