#NEWLINE# def svd_r_raw( ctx , A , V = False , calc_u = False ) : #NEWLINE# #INDENT# "\n    This routine computes the singular value decomposition of a matrix A.\n    Given A, two orthogonal matrices U and V are calculated such that\n\n                    A = U S V\n\n    where S is a suitable shaped matrix whose off-diagonal elements are zero.\n    The diagonal elements of S are the singular values of A, i.e. the\n    squareroots of the eigenvalues of A' A or A A'. Here ' denotes the transpose.\n    Householder bidiagonalization and a variant of the QR algorithm is used.\n\n    overview of the matrices :\n\n      A : m*n       A gets replaced by U\n      U : m*n       U replaces A. If n>m then only the first m*m block of U is\n                    non-zero. column-orthogonal: U' U = B\n                    here B is a n*n matrix whose first min(m,n) diagonal\n                    elements are 1 and all other elements are zero.\n      S : n*n       diagonal matrix, only the diagonal elements are stored in\n                    the array S. only the first min(m,n) diagonal elements are non-zero.\n      V : n*n       orthogonal: V V' = V' V = 1\n\n    parameters:\n      A        (input/output) On input, A contains a real matrix of shape m*n.\n               On output, if calc_u is true A contains the column-orthogonal\n               matrix U; otherwise A is simply used as workspace and thus destroyed.\n\n      V        (input/output) if false, the matrix V is not calculated. otherwise\n               V must be a matrix of shape n*n.\n\n      calc_u   (input) If true, the matrix U is calculated and replaces A.\n               if false, U is not calculated and A is simply destroyed\n\n    return value:\n      S        an array of length n containing the singular values of A sorted by\n               decreasing magnitude. only the first min(m,n) elements are non-zero.\n\n    This routine is a python translation of the fortran routine svd.f in the\n    software library EISPACK (see netlib.org) which itself is based on the\n    algol procedure svd described in:\n      - num. math. 14, 403-420(1970) by golub and reinsch.\n      - wilkinson/reinsch: handbook for auto. comp., vol ii-linear algebra, 134-151(1971).\n\n    " #NEWLINE# ( m , n ) = ( A . rows , A . cols ) #NEWLINE# S = ctx . zeros ( n , 1 ) #NEWLINE# work = ctx . zeros ( n , 1 ) #NEWLINE# g = scale = anorm = 0 #NEWLINE# maxits = ( 3 * ctx . dps ) #NEWLINE# for i in xrange ( n ) : #NEWLINE# #INDENT# work [ i ] = ( scale * g ) #NEWLINE# g = s = scale = 0 #NEWLINE# if ( i < m ) : #NEWLINE# #INDENT# for k in xrange ( i , m ) : #NEWLINE# #INDENT# scale += ctx . fabs ( A [ ( k , i ) ] ) #NEWLINE# #UNINDENT# if ( scale != 0 ) : #NEWLINE# #INDENT# for k in xrange ( i , m ) : #NEWLINE# #INDENT# A [ ( k , i ) ] /= scale #NEWLINE# s += ( A [ ( k , i ) ] * A [ ( k , i ) ] ) #NEWLINE# #UNINDENT# f = A [ ( i , i ) ] #NEWLINE# g = ( - ctx . sqrt ( s ) ) #NEWLINE# if ( f < 0 ) : #NEWLINE# #INDENT# g = ( - g ) #NEWLINE# #UNINDENT# h = ( ( f * g ) - s ) #NEWLINE# A [ ( i , i ) ] = ( f - g ) #NEWLINE# for j in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# s = 0 #NEWLINE# for k in xrange ( i , m ) : #NEWLINE# #INDENT# s += ( A [ ( k , i ) ] * A [ ( k , j ) ] ) #NEWLINE# #UNINDENT# f = ( s / h ) #NEWLINE# for k in xrange ( i , m ) : #NEWLINE# #INDENT# A [ ( k , j ) ] += ( f * A [ ( k , i ) ] ) #NEWLINE# #UNINDENT# #UNINDENT# for k in xrange ( i , m ) : #NEWLINE# #INDENT# A [ ( k , i ) ] *= scale #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# S [ i ] = ( scale * g ) #NEWLINE# g = s = scale = 0 #NEWLINE# if ( ( i < m ) and ( i != ( n - 1 ) ) ) : #NEWLINE# #INDENT# for k in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# scale += ctx . fabs ( A [ ( i , k ) ] ) #NEWLINE# #UNINDENT# if scale : #NEWLINE# #INDENT# for k in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# A [ ( i , k ) ] /= scale #NEWLINE# s += ( A [ ( i , k ) ] * A [ ( i , k ) ] ) #NEWLINE# #UNINDENT# f = A [ ( i , ( i + 1 ) ) ] #NEWLINE# g = ( - ctx . sqrt ( s ) ) #NEWLINE# if ( f < 0 ) : #NEWLINE# #INDENT# g = ( - g ) #NEWLINE# #UNINDENT# h = ( ( f * g ) - s ) #NEWLINE# A [ ( i , ( i + 1 ) ) ] = ( f - g ) #NEWLINE# for k in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# work [ k ] = ( A [ ( i , k ) ] / h ) #NEWLINE# #UNINDENT# for j in xrange ( ( i + 1 ) , m ) : #NEWLINE# #INDENT# s = 0 #NEWLINE# for k in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# s += ( A [ ( j , k ) ] * A [ ( i , k ) ] ) #NEWLINE# #UNINDENT# for k in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# A [ ( j , k ) ] += ( s * work [ k ] ) #NEWLINE# #UNINDENT# #UNINDENT# for k in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# A [ ( i , k ) ] *= scale #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# anorm = max ( anorm , ( ctx . fabs ( S [ i ] ) + ctx . fabs ( work [ i ] ) ) ) #NEWLINE# #UNINDENT# if ( not isinstance ( V , bool ) ) : #NEWLINE# #INDENT# for i in xrange ( ( n - 2 ) , ( - 1 ) , ( - 1 ) ) : #NEWLINE# #INDENT# V [ ( ( i + 1 ) , ( i + 1 ) ) ] = 1 #NEWLINE# if ( work [ ( i + 1 ) ] != 0 ) : #NEWLINE# #INDENT# for j in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# V [ ( i , j ) ] = ( ( A [ ( i , j ) ] / A [ ( i , ( i + 1 ) ) ] ) / work [ ( i + 1 ) ] ) #NEWLINE# #UNINDENT# for j in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# s = 0 #NEWLINE# for k in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# s += ( A [ ( i , k ) ] * V [ ( j , k ) ] ) #NEWLINE# #UNINDENT# for k in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# V [ ( j , k ) ] += ( s * V [ ( i , k ) ] ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# for j in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# V [ ( j , i ) ] = V [ ( i , j ) ] = 0 #NEWLINE# #UNINDENT# #UNINDENT# V [ ( 0 , 0 ) ] = 1 #NEWLINE# #UNINDENT# if ( m < n ) : #NEWLINE# #INDENT# minnm = m #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# minnm = n #NEWLINE# #UNINDENT# if calc_u : #NEWLINE# #INDENT# for i in xrange ( ( minnm - 1 ) , ( - 1 ) , ( - 1 ) ) : #NEWLINE# #INDENT# g = S [ i ] #NEWLINE# for j in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# A [ ( i , j ) ] = 0 #NEWLINE# #UNINDENT# if ( g != 0 ) : #NEWLINE# #INDENT# g = ( 1 / g ) #NEWLINE# for j in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# s = 0 #NEWLINE# for k in xrange ( ( i + 1 ) , m ) : #NEWLINE# #INDENT# s += ( A [ ( k , i ) ] * A [ ( k , j ) ] ) #NEWLINE# #UNINDENT# f = ( ( s / A [ ( i , i ) ] ) * g ) #NEWLINE# for k in xrange ( i , m ) : #NEWLINE# #INDENT# A [ ( k , j ) ] += ( f * A [ ( k , i ) ] ) #NEWLINE# #UNINDENT# #UNINDENT# for j in xrange ( i , m ) : #NEWLINE# #INDENT# A [ ( j , i ) ] *= g #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# for j in xrange ( i , m ) : #NEWLINE# #INDENT# A [ ( j , i ) ] = 0 #NEWLINE# #UNINDENT# #UNINDENT# A [ ( i , i ) ] += 1 #NEWLINE# #UNINDENT# #UNINDENT# for k in xrange ( ( n - 1 ) , ( - 1 ) , ( - 1 ) ) : #NEWLINE# #INDENT# its = 0 #NEWLINE# while 1 : #NEWLINE# #INDENT# its += 1 #NEWLINE# flag = True #NEWLINE# for l in xrange ( k , ( - 1 ) , ( - 1 ) ) : #NEWLINE# #INDENT# nm = ( l - 1 ) #NEWLINE# if ( ( ctx . fabs ( work [ l ] ) + anorm ) == anorm ) : #NEWLINE# #INDENT# flag = False #NEWLINE# break #NEWLINE# #UNINDENT# if ( ( ctx . fabs ( S [ nm ] ) + anorm ) == anorm ) : #NEWLINE# #INDENT# break #NEWLINE# #UNINDENT# #UNINDENT# if flag : #NEWLINE# #INDENT# c = 0 #NEWLINE# s = 1 #NEWLINE# for i in xrange ( l , ( k + 1 ) ) : #NEWLINE# #INDENT# f = ( s * work [ i ] ) #NEWLINE# work [ i ] *= c #NEWLINE# if ( ( ctx . fabs ( f ) + anorm ) == anorm ) : #NEWLINE# #INDENT# break #NEWLINE# #UNINDENT# g = S [ i ] #NEWLINE# h = ctx . hypot ( f , g ) #NEWLINE# S [ i ] = h #NEWLINE# h = ( 1 / h ) #NEWLINE# c = ( g * h ) #NEWLINE# s = ( ( - f ) * h ) #NEWLINE# if calc_u : #NEWLINE# #INDENT# for j in xrange ( m ) : #NEWLINE# #INDENT# y = A [ ( scale , nm ) ] #NEWLINE# z = A [ ( j , i ) ] #NEWLINE# A [ ( j , nm ) ] = ( ( y * c ) + ( z * s ) ) #NEWLINE# A [ ( j , i ) ] = ( ( z * c ) - ( y * s ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# z = S [ k ] #NEWLINE# if ( l == k ) : #NEWLINE# #INDENT# if ( z < 0 ) : #NEWLINE# #INDENT# S [ k ] = ( - z ) #NEWLINE# if ( not isinstance ( V , bool ) ) : #NEWLINE# #INDENT# for j in xrange ( n ) : #NEWLINE# #INDENT# V [ ( k , j ) ] = ( - V [ ( k , j ) ] ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# break #NEWLINE# #UNINDENT# if ( its >= maxits ) : #NEWLINE# #INDENT# raise RuntimeError ( ( 'svd: no convergence to an eigenvalue after %d iterations' % its ) ) #NEWLINE# #UNINDENT# x = S [ l ] #NEWLINE# nm = ( k - 1 ) #NEWLINE# y = S [ nm ] #NEWLINE# g = work [ nm ] #NEWLINE# h = work [ k ] #NEWLINE# f = ( ( ( ( y - z ) * ( y + z ) ) + ( ( g - h ) * ( g + h ) ) ) / ( ( 2 * h ) * y ) ) #NEWLINE# g = ctx . hypot ( f , 1 ) #NEWLINE# if ( f >= 0 ) : #NEWLINE# #INDENT# f = ( ( ( ( x - z ) * ( x + z ) ) + ( h * ( ( y / ( f + g ) ) - h ) ) ) / x ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# f = ( ( ( ( x - z ) * ( x + z ) ) + ( h * ( ( y / ( f - g ) ) - h ) ) ) / x ) #NEWLINE# #UNINDENT# c = s = 1 #NEWLINE# for j in xrange ( l , ( nm + 1 ) ) : #NEWLINE# #INDENT# g = work [ ( j + 1 ) ] #NEWLINE# y = S [ ( j + 1 ) ] #NEWLINE# h = ( s * g ) #NEWLINE# g = ( c * g ) #NEWLINE# z = ctx . hypot ( f , h ) #NEWLINE# work [ j ] = z #NEWLINE# c = ( f / z ) #NEWLINE# s = ( h / z ) #NEWLINE# f = ( ( x * c ) + ( g * s ) ) #NEWLINE# g = ( ( g * c ) - ( x * s ) ) #NEWLINE# h = ( y * s ) #NEWLINE# y *= c #NEWLINE# if ( not isinstance ( V , bool ) ) : #NEWLINE# #INDENT# for jj in xrange ( n ) : #NEWLINE# #INDENT# x = V [ ( j , jj ) ] #NEWLINE# z = V [ ( ( j + 1 ) , jj ) ] #NEWLINE# V [ ( j , jj ) ] = ( ( x * c ) + ( z * s ) ) #NEWLINE# V [ ( ( j + 1 ) , jj ) ] = ( ( z * c ) - ( x * s ) ) #NEWLINE# #UNINDENT# #UNINDENT# z = ctx . hypot ( f , h ) #NEWLINE# S [ j ] = z #NEWLINE# if ( z != 0 ) : #NEWLINE# #INDENT# z = ( 1 / z ) #NEWLINE# c = ( f * z ) #NEWLINE# s = ( h * z ) #NEWLINE# #UNINDENT# f = ( ( c * g ) + ( s * y ) ) #NEWLINE# x = ( ( c * y ) - ( s * g ) ) #NEWLINE# if calc_u : #NEWLINE# #INDENT# for jj in xrange ( m ) : #NEWLINE# #INDENT# y = A [ ( jj , j ) ] #NEWLINE# z = A [ ( jj , ( j + 1 ) ) ] #NEWLINE# A [ ( jj , j ) ] = ( ( y * c ) + ( z * s ) ) #NEWLINE# A [ ( jj , ( j + 1 ) ) ] = ( ( z * c ) - ( y * s ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# work [ l ] = 0 #NEWLINE# work [ k ] = f #NEWLINE# S [ k ] = x #NEWLINE# #UNINDENT# #UNINDENT# for i in xrange ( n ) : #NEWLINE# #INDENT# imax = i #NEWLINE# s = ctx . fabs ( S [ i ] ) #NEWLINE# for j in xrange ( ( i + 1 ) , n ) : #NEWLINE# #INDENT# c = ctx . fabs ( S [ j ] ) #NEWLINE# if ( c > s ) : #NEWLINE# #INDENT# s = c #NEWLINE# imax = j #NEWLINE# #UNINDENT# #UNINDENT# if ( imax != i ) : #NEWLINE# #INDENT# z = S [ i ] #NEWLINE# S [ i ] = S [ imax ] #NEWLINE# S [ imax ] = z #NEWLINE# if calc_u : #NEWLINE# #INDENT# for j in xrange ( m ) : #NEWLINE# #INDENT# z = A [ ( j , i ) ] #NEWLINE# A [ ( j , i ) ] = A [ ( j , imax ) ] #NEWLINE# A [ ( j , imax ) ] = z #NEWLINE# #UNINDENT# #UNINDENT# if ( not isinstance ( V , bool ) ) : #NEWLINE# #INDENT# for j in xrange ( n ) : #NEWLINE# #INDENT# z = V [ ( i , j ) ] #NEWLINE# V [ ( i , j ) ] = V [ ( imax , j ) ] #NEWLINE# V [ ( imax , j ) ] = z #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# return S