(ctx, A, V = False, calc_u = False)
return S
This routine computes the singular value decomposition of a matrix A.
    Given A, two orthogonal matrices U and V are calculated such that

                    A = U S V

    where S is a suitable shaped matrix whose off-diagonal elements are zero.
    The diagonal elements of S are the singular values of A, i.e. the
    squareroots of the eigenvalues of A' A or A A'. Here ' denotes the transpose.
    Householder bidiagonalization and a variant of the QR algorithm is used.

    overview of the matrices :

      A : m*n       A gets replaced by U
      U : m*n       U replaces A. If n>m then only the first m*m block of U is
                    non-zero. column-orthogonal: U' U = B
                    here B is a n*n matrix whose first min(m,n) diagonal
                    elements are 1 and all other elements are zero.
      S : n*n       diagonal matrix, only the diagonal elements are stored in
                    the array S. only the first min(m,n) diagonal elements are non-zero.
      V : n*n       orthogonal: V V' = V' V = 1

    parameters:
      A        (input/output) On input, A contains a real matrix of shape m*n.
               On output, if calc_u is true A contains the column-orthogonal
               matrix U; otherwise A is simply used as workspace and thus destroyed.

      V        (input/output) if false, the matrix V is not calculated. otherwise
               V must be a matrix of shape n*n.

      calc_u   (input) If true, the matrix U is calculated and replaces A.
               if false, U is not calculated and A is simply destroyed

    return value:
      S        an array of length n containing the singular values of A sorted by
               decreasing magnitude. only the first min(m,n) elements are non-zero.

    This routine is a python translation of the fortran routine svd.f in the
    software library EISPACK (see netlib.org) which itself is based on the
    algol procedure svd described in:
      - num. math. 14, 403-420(1970) by golub and reinsch.
      - wilkinson/reinsch: handbook for auto. comp., vol ii-linear algebra, 134-151(1971).
This routine computes the singular value decomposition of a matrix A.
    Given A, two orthogonal matrices U and V are calculated such that
This routine computes the singular value decomposition of a matrix A . Given A two orthogonal matrices U and V are calculated such that
def svd_r_raw(ctx, A, V = False, calc_u = False):
    """
    This routine computes the singular value decomposition of a matrix A.
    Given A, two orthogonal matrices U and V are calculated such that

                    A = U S V

    where S is a suitable shaped matrix whose off-diagonal elements are zero.
    The diagonal elements of S are the singular values of A, i.e. the
    squareroots of the eigenvalues of A' A or A A'. Here ' denotes the transpose.
    Householder bidiagonalization and a variant of the QR algorithm is used.

    overview of the matrices :

      A : m*n       A gets replaced by U
      U : m*n       U replaces A. If n>m then only the first m*m block of U is
                    non-zero. column-orthogonal: U' U = B
                    here B is a n*n matrix whose first min(m,n) diagonal
                    elements are 1 and all other elements are zero.
      S : n*n       diagonal matrix, only the diagonal elements are stored in
                    the array S. only the first min(m,n) diagonal elements are non-zero.
      V : n*n       orthogonal: V V' = V' V = 1

    parameters:
      A        (input/output) On input, A contains a real matrix of shape m*n.
               On output, if calc_u is true A contains the column-orthogonal
               matrix U; otherwise A is simply used as workspace and thus destroyed.

      V        (input/output) if false, the matrix V is not calculated. otherwise
               V must be a matrix of shape n*n.

      calc_u   (input) If true, the matrix U is calculated and replaces A.
               if false, U is not calculated and A is simply destroyed

    return value:
      S        an array of length n containing the singular values of A sorted by
               decreasing magnitude. only the first min(m,n) elements are non-zero.

    This routine is a python translation of the fortran routine svd.f in the
    software library EISPACK (see netlib.org) which itself is based on the
    algol procedure svd described in:
      - num. math. 14, 403-420(1970) by golub and reinsch.
      - wilkinson/reinsch: handbook for auto. comp., vol ii-linear algebra, 134-151(1971).

    """

    m, n = A.rows, A.cols

    S = ctx.zeros(n, 1)

    # work is a temporary array of size n
    work = ctx.zeros(n, 1)

    g = scale = anorm = 0
    maxits = 3 * ctx.dps

    for i in xrange(n):     # householder reduction to bidiagonal form
        work[i] = scale*g
        g = s = scale = 0
        if i < m:
            for k in xrange(i, m):
                scale += ctx.fabs(A[k,i])
            if scale != 0:
                for k in xrange(i, m):
                    A[k,i] /= scale
                    s += A[k,i] * A[k,i]
                f = A[i,i]
                g = -ctx.sqrt(s)
                if f < 0:
                    g = -g
                h = f * g - s
                A[i,i] = f - g
                for j in xrange(i+1, n):
                    s = 0
                    for k in xrange(i, m):
                        s += A[k,i] * A[k,j]
                    f = s / h
                    for k in xrange(i, m):
                        A[k,j] += f * A[k,i]
                for k in xrange(i,m):
                    A[k,i] *= scale

        S[i] = scale * g
        g = s = scale = 0

        if i < m and i != n - 1:
            for k in xrange(i+1, n):
                scale += ctx.fabs(A[i,k])
            if scale:
                for k in xrange(i+1, n):
                    A[i,k] /= scale
                    s += A[i,k] * A[i,k]
                f = A[i,i+1]
                g = -ctx.sqrt(s)
                if f < 0:
                    g = -g
                h = f * g - s
                A[i,i+1] = f - g

                for k in xrange(i+1, n):
                    work[k] = A[i,k] / h

                for j in xrange(i+1, m):
                    s = 0
                    for k in xrange(i+1, n):
                        s += A[j,k] * A[i,k]
                    for k in xrange(i+1, n):
                        A[j,k] += s * work[k]

                for k in xrange(i+1, n):
                    A[i,k] *= scale

        anorm = max(anorm, ctx.fabs(S[i]) + ctx.fabs(work[i]))

    if not isinstance(V, bool):
        for i in xrange(n-2, -1, -1):     # accumulation of right hand transformations
            V[i+1,i+1] = 1

            if work[i+1] != 0:
                for j in xrange(i+1, n):
                    V[i,j] = (A[i,j] / A[i,i+1]) / work[i+1]
                for j in xrange(i+1, n):
                    s = 0
                    for k in xrange(i+1, n):
                        s += A[i,k] * V[j,k]
                    for k in xrange(i+1, n):
                        V[j,k] += s * V[i,k]

            for j in xrange(i+1, n):
                V[j,i] = V[i,j] = 0

        V[0,0] = 1

    if m<n : minnm = m
    else   : minnm = n

    if calc_u:
        for i in xrange(minnm-1, -1, -1): # accumulation of left hand transformations
            g = S[i]
            for j in xrange(i+1, n):
                A[i,j] = 0
            if g != 0:
                g = 1 / g
                for j in xrange(i+1, n):
                    s = 0
                    for k in xrange(i+1, m):
                        s += A[k,i] * A[k,j]
                    f = (s / A[i,i]) * g
                    for k in xrange(i, m):
                        A[k,j] += f * A[k,i]
                for j in xrange(i, m):
                    A[j,i] *= g
            else:
                for j in xrange(i, m):
                    A[j,i] = 0
            A[i,i] += 1

    for k in xrange(n - 1, -1, -1):
        # diagonalization of the bidiagonal form:
        #   loop over singular values, and over allowed itations

        its = 0
        while 1:
            its += 1
            flag = True

            for l in xrange(k, -1, -1):
                nm = l-1

                if ctx.fabs(work[l]) + anorm == anorm:
                    flag = False
                    break

                if ctx.fabs(S[nm]) + anorm == anorm:
                    break

            if flag:
                c = 0
                s = 1
                for i in xrange(l, k + 1):
                    f = s * work[i]
                    work[i] *= c
                    if ctx.fabs(f) + anorm == anorm:
                        break
                    g = S[i]
                    h = ctx.hypot(f, g)
                    S[i] = h
                    h = 1 / h
                    c = g * h
                    s = - f * h

                    if calc_u:
                        for j in xrange(m):
                            y = A[j,nm]
                            z = A[j,i]
                            A[j,nm] = y * c + z * s
                            A[j,i]  = z * c - y * s

            z = S[k]

            if l == k:               # convergence
                if z < 0:            # singular value is made nonnegative
                    S[k] = -z
                    if not isinstance(V, bool):
                        for j in xrange(n):
                            V[k,j] = -V[k,j]
                break

            if its >= maxits:
                raise RuntimeError("svd: no convergence to an eigenvalue after %d iterations" % its)

            x = S[l]         # shift from bottom 2 by 2 minor
            nm = k-1
            y = S[nm]
            g = work[nm]
            h = work[k]
            f = ((y - z) * (y + z) + (g - h) * (g + h))/(2 * h * y)
            g = ctx.hypot(f, 1)
            if f >= 0: f = ((x - z) * (x + z) + h * ((y / (f + g)) - h)) / x
            else:      f = ((x - z) * (x + z) + h * ((y / (f - g)) - h)) / x

            c = s = 1         # next qt transformation

            for j in xrange(l, nm + 1):
                g = work[j+1]
                y = S[j+1]
                h = s * g
                g = c * g
                z = ctx.hypot(f, h)
                work[j] = z
                c = f / z
                s = h / z
                f = x * c + g * s
                g = g * c - x * s
                h = y * s
                y *= c
                if not isinstance(V, bool):
                    for jj in xrange(n):
                        x = V[j  ,jj]
                        z = V[j+1,jj]
                        V[j    ,jj]= x * c + z * s
                        V[j+1  ,jj]= z * c - x * s
                z = ctx.hypot(f, h)
                S[j] = z
                if z != 0:            # rotation can be arbitray if z=0
                    z = 1 / z
                    c = f * z
                    s = h * z
                f = c * g + s * y
                x = c * y - s * g

                if calc_u:
                    for jj in xrange(m):
                        y = A[jj,j  ]
                        z = A[jj,j+1]
                        A[jj,j    ] = y * c + z * s
                        A[jj,j+1  ] = z * c - y * s

            work[l] = 0
            work[k] = f
            S[k] = x

    ##########################

    # Sort singular values into decreasing order (bubble-sort)

    for i in xrange(n):
        imax = i
        s = ctx.fabs(S[i])         # s is the current maximal element

        for j in xrange(i + 1, n):
            c = ctx.fabs(S[j])
            if c > s:
                s = c
                imax = j

        if imax != i:
            # swap singular values

            z = S[i]
            S[i] = S[imax]
            S[imax] = z

            if calc_u:
                for j in xrange(m):
                    z = A[j,i]
                    A[j,i] = A[j,imax]
                    A[j,imax] = z

            if not isinstance(V, bool):
                for j in xrange(n):
                    z = V[i,j]
                    V[i,j] = V[imax,j]
                    V[imax,j] = z

    return S
def svd_r_raw ( ctx , A , V = False , calc_u = False ) : m , n = A . rows , A . cols S = ctx . zeros ( n , 1 ) # work is a temporary array of size n work = ctx . zeros ( n , 1 ) g = scale = anorm = 0 maxits = 3 * ctx . dps for i in xrange ( n ) : # householder reduction to bidiagonal form work [ i ] = scale * g g = s = scale = 0 if i < m : for k in xrange ( i , m ) : scale += ctx . fabs ( A [ k , i ] ) if scale != 0 : for k in xrange ( i , m ) : A [ k , i ] /= scale s += A [ k , i ] * A [ k , i ] f = A [ i , i ] g = - ctx . sqrt ( s ) if f < 0 : g = - g h = f * g - s A [ i , i ] = f - g for j in xrange ( i + 1 , n ) : s = 0 for k in xrange ( i , m ) : s += A [ k , i ] * A [ k , j ] f = s / h for k in xrange ( i , m ) : A [ k , j ] += f * A [ k , i ] for k in xrange ( i , m ) : A [ k , i ] *= scale S [ i ] = scale * g g = s = scale = 0 if i < m and i != n - 1 : for k in xrange ( i + 1 , n ) : scale += ctx . fabs ( A [ i , k ] ) if scale : for k in xrange ( i + 1 , n ) : A [ i , k ] /= scale s += A [ i , k ] * A [ i , k ] f = A [ i , i + 1 ] g = - ctx . sqrt ( s ) if f < 0 : g = - g h = f * g - s A [ i , i + 1 ] = f - g for k in xrange ( i + 1 , n ) : work [ k ] = A [ i , k ] / h for j in xrange ( i + 1 , m ) : s = 0 for k in xrange ( i + 1 , n ) : s += A [ j , k ] * A [ i , k ] for k in xrange ( i + 1 , n ) : A [ j , k ] += s * work [ k ] for k in xrange ( i + 1 , n ) : A [ i , k ] *= scale anorm = max ( anorm , ctx . fabs ( S [ i ] ) + ctx . fabs ( work [ i ] ) ) if not isinstance ( V , bool ) : for i in xrange ( n - 2 , - 1 , - 1 ) : # accumulation of right hand transformations V [ i + 1 , i + 1 ] = 1 if work [ i + 1 ] != 0 : for j in xrange ( i + 1 , n ) : V [ i , j ] = ( A [ i , j ] / A [ i , i + 1 ] ) / work [ i + 1 ] for j in xrange ( i + 1 , n ) : s = 0 for k in xrange ( i + 1 , n ) : s += A [ i , k ] * V [ j , k ] for k in xrange ( i + 1 , n ) : V [ j , k ] += s * V [ i , k ] for j in xrange ( i + 1 , n ) : V [ j , i ] = V [ i , j ] = 0 V [ 0 , 0 ] = 1 if m < n : minnm = m else : minnm = n if calc_u : for i in xrange ( minnm - 1 , - 1 , - 1 ) : # accumulation of left hand transformations g = S [ i ] for j in xrange ( i + 1 , n ) : A [ i , j ] = 0 if g != 0 : g = 1 / g for j in xrange ( i + 1 , n ) : s = 0 for k in xrange ( i + 1 , m ) : s += A [ k , i ] * A [ k , j ] f = ( s / A [ i , i ] ) * g for k in xrange ( i , m ) : A [ k , j ] += f * A [ k , i ] for j in xrange ( i , m ) : A [ j , i ] *= g else : for j in xrange ( i , m ) : A [ j , i ] = 0 A [ i , i ] += 1 for k in xrange ( n - 1 , - 1 , - 1 ) : # diagonalization of the bidiagonal form: #   loop over singular values, and over allowed itations its = 0 while 1 : its += 1 flag = True for l in xrange ( k , - 1 , - 1 ) : nm = l - 1 if ctx . fabs ( work [ l ] ) + anorm == anorm : flag = False break if ctx . fabs ( S [ nm ] ) + anorm == anorm : break if flag : c = 0 s = 1 for i in xrange ( l , k + 1 ) : f = s * work [ i ] work [ i ] *= c if ctx . fabs ( f ) + anorm == anorm : break g = S [ i ] h = ctx . hypot ( f , g ) S [ i ] = h h = 1 / h c = g * h s = - f * h if calc_u : for j in xrange ( m ) : y = A [ j , nm ] z = A [ j , i ] A [ j , nm ] = y * c + z * s A [ j , i ] = z * c - y * s z = S [ k ] if l == k : # convergence if z < 0 : # singular value is made nonnegative S [ k ] = - z if not isinstance ( V , bool ) : for j in xrange ( n ) : V [ k , j ] = - V [ k , j ] break if its >= maxits : raise RuntimeError ( "svd: no convergence to an eigenvalue after %d iterations" % its ) x = S [ l ] # shift from bottom 2 by 2 minor nm = k - 1 y = S [ nm ] g = work [ nm ] h = work [ k ] f = ( ( y - z ) * ( y + z ) + ( g - h ) * ( g + h ) ) / ( 2 * h * y ) g = ctx . hypot ( f , 1 ) if f >= 0 : f = ( ( x - z ) * ( x + z ) + h * ( ( y / ( f + g ) ) - h ) ) / x else : f = ( ( x - z ) * ( x + z ) + h * ( ( y / ( f - g ) ) - h ) ) / x c = s = 1 # next qt transformation for j in xrange ( l , nm + 1 ) : g = work [ j + 1 ] y = S [ j + 1 ] h = s * g g = c * g z = ctx . hypot ( f , h ) work [ j ] = z c = f / z s = h / z f = x * c + g * s g = g * c - x * s h = y * s y *= c if not isinstance ( V , bool ) : for jj in xrange ( n ) : x = V [ j , jj ] z = V [ j + 1 , jj ] V [ j , jj ] = x * c + z * s V [ j + 1 , jj ] = z * c - x * s z = ctx . hypot ( f , h ) S [ j ] = z if z != 0 : # rotation can be arbitray if z=0 z = 1 / z c = f * z s = h * z f = c * g + s * y x = c * y - s * g if calc_u : for jj in xrange ( m ) : y = A [ jj , j ] z = A [ jj , j + 1 ] A [ jj , j ] = y * c + z * s A [ jj , j + 1 ] = z * c - y * s work [ l ] = 0 work [ k ] = f S [ k ] = x ########################## # Sort singular values into decreasing order (bubble-sort) for i in xrange ( n ) : imax = i s = ctx . fabs ( S [ i ] ) # s is the current maximal element for j in xrange ( i + 1 , n ) : c = ctx . fabs ( S [ j ] ) if c > s : s = c imax = j if imax != i : # swap singular values z = S [ i ] S [ i ] = S [ imax ] S [ imax ] = z if calc_u : for j in xrange ( m ) : z = A [ j , i ] A [ j , i ] = A [ j , imax ] A [ j , imax ] = z if not isinstance ( V , bool ) : for j in xrange ( n ) : z = V [ i , j ] V [ i , j ] = V [ imax , j ] V [ imax , j ] = z return S