#NEWLINE# def eig_tr_r( ctx , A ) : #NEWLINE# #INDENT# '\n    This routine calculates the right eigenvectors of an upper right triangular matrix.\n\n    input:\n      A      an upper right triangular matrix\n\n    output:\n      ER     a matrix whose columns form the right eigenvectors of A\n\n    return value: ER\n    ' #NEWLINE# n = A . rows #NEWLINE# ER = ctx . eye ( n ) #NEWLINE# eps = ctx . eps #NEWLINE# unfl = ctx . ldexp ( ctx . one , ( ( - ctx . prec ) * 30 ) ) #NEWLINE# smlnum = ( unfl * ( n / eps ) ) #NEWLINE# simin = ( 1 / ctx . sqrt ( eps ) ) #NEWLINE# rmax = 1 #NEWLINE# for i in xrange ( 1 , n ) : #NEWLINE# #INDENT# s = A [ ( i , i ) ] #NEWLINE# smin = max ( ( eps * abs ( s ) ) , smlnum ) #NEWLINE# for j in xrange ( ( i - 1 ) , ( - 1 ) , ( - 1 ) ) : #NEWLINE# #INDENT# r = 0 #NEWLINE# for k in xrange ( ( j + 1 ) , ( i + 1 ) ) : #NEWLINE# #INDENT# r += ( A [ ( j , k ) ] * ER [ ( k , i ) ] ) #NEWLINE# #UNINDENT# t = ( A [ ( j , j ) ] - s ) #NEWLINE# if ( abs ( t ) < smin ) : #NEWLINE# #INDENT# t = smin #NEWLINE# #UNINDENT# r = ( ( - r ) / t ) #NEWLINE# ER [ ( j , i ) ] = r #NEWLINE# rmax = max ( rmax , abs ( r ) ) #NEWLINE# if ( rmax > simin ) : #NEWLINE# #INDENT# for k in xrange ( j , ( i + 1 ) ) : #NEWLINE# #INDENT# ER [ ( k , i ) ] /= rmax #NEWLINE# #UNINDENT# rmax = 1 #NEWLINE# #UNINDENT# #UNINDENT# if ( rmax != 1 ) : #NEWLINE# #INDENT# for k in xrange ( 0 , ( i + 1 ) ) : #NEWLINE# #INDENT# ER [ ( k , i ) ] /= rmax #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# return ER