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.)
# insert code here
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, you should get a matrix of (approximately) all zeroes.
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.