Sources Esope | Cast3M

invsvd
C INVSVD    SOURCE    CB215821  26/08/24    21:16:55     12622          
      SUBROUTINE INVSVD(NM,M,N,A,W,MATU,U,MATV,V,iarr,RV1)
c
      IMPLICIT INTEGER(I-N)
      integer i,j,k,l,m,n,ii,i1,kk,k1,ll,l1,mn,nm,its,iarr
      real*8 a(nm,n),w(n),u(nm,n),v(nm,n),rv1(n)
      real*8 c,f,g,h,s,x,y,z,scale,anorm
      real*8 sqrt,max,abs,sign
      logical matu,matv
      real*8 zero, one, two
      parameter(zero=0.0D0, one=1.0D0, two=2.0D0)
c
c     Adaptation to Esope
c     A. BECCANTINI
c     DRN/DMT/SEMT/LTMF
c     18.08.00
c
c     this subroutine is a translation of the algol procedure svd,
c     num. math. 14, 403-420(1970) by golub and reinsch.
c     handbook for auto. comp., vol ii-linear algebra, 134-151(1971).
c
c     this subroutine determines the singular value decomposition
c          t
c     a=usv  of a real m by n rectangular matrix.  householder
c     bidiagonalization and a variant of the qr algorithm are used.
c
c     on input.
c
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement.  note that nm must be at least
c          as large as the maximum of m and n.
c
c        m is the number of rows of a (and u).
c
c        n is the number of columns of a (and u) and the order of v.
c
c        a contains the rectangular input matrix to be decomposed.
c
c        matu should be set to .true. if the u matrix in the
c          decomposition is desired, and to .false. otherwise.
c
c        matv should be set to .true. if the v matrix in the
c          decomposition is desired, and to .false. otherwise.
c
c     on output.
c
c        a is unaltered (unless overwritten by u or v).
c
c        w contains the n (non-negative) singular values of a (the
c          diagonal elements of s).  they are unordered.  if an
c          error exit is made, the singular values should be correct
c          for indices iarr+1,iarr+2,...,n.
c
c        u contains the matrix u (orthogonal column vectors) of the
c          decomposition if matu has been set to .true.  otherwise
c          u is used as a temporary array.  u may coincide with a.
c          if an error exit is made, the columns of u corresponding
c          to indices of correct singular values should be correct.
c
c        v contains the matrix v (orthogonal) of the decomposition if
c          matv has been set to .true.  otherwise v is not referenced.
c          v may also coincide with a if u is not needed.  if an error
c          exit is made, the columns of v corresponding to indices of
c          correct singular values should be correct.
c
c        iarr is set to
c          zero       for normal return,
c          k          if the k-th singular value has not been
c                     determined after 30 iterations.
c
c        rv1 is a temporary storage array.
c
c     this is a modified version of a routine from the eispack
c     collection by the nats project
c
c     modified to eliminate machep
c
      iarr = 0
c
      do 1002 i = 1, m
c
         do 100 j = 1, n
            u(i,j) = a(i,j)
  100 continue
 1002 CONTINUE
c     .......... householder reduction to bidiagonal form ..........
      g = zero
      scale = zero
      anorm = zero
c
      do 300 i = 1, n
         l = i + 1
         rv1(i) = scale * g
         g = zero
         s = zero
         scale = zero
         if (i .gt. m) go to 210
c
         do 120 k = i, m
         scale = scale + abs(u(k,i))
  120    CONTINUE
c
         if (scale .eq. zero) go to 210
c
         do 130 k = i, m
            u(k,i) = u(k,i) / scale
            s = s + u(k,i)**2
  130    continue
c
         f = u(i,i)
         g = -sign(sqrt(s),f)
         h = f * g - s
         u(i,i) = f - g
         if (i .eq. n) go to 190
c
         do 1003 j = l, n
            s = zero
c
            do 140 k = i, m
            s = s + u(k,i) * u(k,j)
  140       CONTINUE
c
            f = s / h
c
            do 150 k = i, m
               u(k,j) = u(k,j) + f * u(k,i)
  150    continue
 1003    CONTINUE
c
  190    do 200 k = i, m
  200    u(k,i) = scale * u(k,i)
c
  210    w(i) = scale * g
         g = zero
         s = zero
         scale = zero
         if (i .gt. m .or. i .eq. n) go to 290
c
         do 220 k = l, n
         scale = scale + abs(u(i,k))
  220    CONTINUE
c
         if (scale .eq. zero) go to 290
c
         do 230 k = l, n
            u(i,k) = u(i,k) / scale
            s = s + u(i,k)**2
  230    continue
c
         f = u(i,l)
         g = -sign(sqrt(s),f)
         h = f * g - s
         u(i,l) = f - g
c
         do 240 k = l, n
         rv1(k) = u(i,k) / h
  240    CONTINUE
c
         if (i .eq. m) go to 270
c
         do 1004 j = l, m
            s = zero
c
            do 250 k = l, n
            s = s + u(j,k) * u(i,k)
  250       CONTINUE
c
            do 260 k = l, n
               u(j,k) = u(j,k) + s * rv1(k)
  260    continue
 1004    CONTINUE
c
  270    do 280 k = l, n
  280    u(i,k) = scale * u(i,k)
c
  290    anorm = max(anorm,abs(w(i))+abs(rv1(i)))
  300 continue
c     .......... accumulation of right-hand transformations ..........
      if (.not. matv) go to 410
c     .......... for i=n step -1 until 1 do -- ..........
      do 400 ii = 1, n
         i = n + 1 - ii
         if (i .eq. n) go to 390
         if (g .eq. zero) go to 360
c
         do 320 j = l, n
c     .......... double division avoids possible underflow ..........
         v(j,i) = (u(i,j) / u(i,l)) / g
  320    CONTINUE
c
         do 1005 j = l, n
            s = zero
c
            do 340 k = l, n
            s = s + u(i,k) * v(k,j)
  340       CONTINUE
c
            do 350 k = l, n
               v(k,j) = v(k,j) + s * v(k,i)
  350    continue
 1005    CONTINUE
c
  360    do 380 j = l, n
            v(i,j) = zero
            v(j,i) = zero
  380    continue
c
  390    v(i,i) = one
         g = rv1(i)
         l = i
  400 continue
c     .......... accumulation of left-hand transformations ..........
  410 if (.not. matu) go to 510
c     ..........for i=min(m,n) step -1 until 1 do -- ..........
      mn = n
      if (m .lt. n) mn = m
c
      do 500 ii = 1, mn
         i = mn + 1 - ii
         l = i + 1
         g = w(i)
         if (i .eq. n) go to 430
c
         do 420 j = l, n
         u(i,j) = zero
  420    CONTINUE
c
  430    if (g .eq. zero) go to 475
         if (i .eq. mn) go to 460
c
         do 1006 j = l, n
            s = zero
c
            do 440 k = l, m
            s = s + u(k,i) * u(k,j)
  440       CONTINUE
c     .......... double division avoids possible underflow ..........
            f = (s / u(i,i)) / g
c
            do 450 k = i, m
               u(k,j) = u(k,j) + f * u(k,i)
  450    continue
 1006    CONTINUE
c
  460    do 470 j = i, m
  470    u(j,i) = u(j,i) / g
c
         go to 490
c
  475    do 480 j = i, m
  480    u(j,i) = zero
c
  490    u(i,i) = u(i,i) + one
  500 continue
c     .......... diagonalization of the bidiagonal form ..........
c     .......... for k=n step -1 until 1 do -- ..........
  510 do 700 kk = 1, n
         k1 = n - kk
         k = k1 + 1
         its = 0
c     .......... test for splitting.
c                for l=k step -1 until 1 do -- ..........
  520    do 530 ll = 1, k
            l1 = k - ll
            l = l1 + 1
            if (abs(rv1(l)) + anorm .eq. anorm) go to 565
c     .......... rv1(1) is always zero, so there is no exit
c                through the bottom of the loop ..........
            if (abs(w(l1)) + anorm .eq. anorm) go to 540
  530    continue
c     .......... cancellation of rv1(l) if l greater than 1 ..........
  540    c = zero
         s = one
c
         do 560 i = l, k
            f = s * rv1(i)
            rv1(i) = c * rv1(i)
            if (abs(f) + anorm .eq. anorm) go to 565
            g = w(i)
            h = sqrt(f*f+g*g)
            w(i) = h
            c = g / h
            s = -f / h
            if (.not. matu) go to 560
c
            do 550 j = 1, m
               y = u(j,l1)
               z = u(j,i)
               u(j,l1) = y * c + z * s
               u(j,i) = -y * s + z * c
  550       continue
c
  560    continue
c     .......... test for convergence ..........
  565    z = w(k)
         if (l .eq. k) go to 650
c     .......... shift from bottom 2 by 2 minor ..........
         if (its .eq. 30) go to 1000
         its = its + 1
         x = w(l)
         y = w(k1)
         g = rv1(k1)
         h = rv1(k)
         f = ((y - z) * (y + z) + (g - h) * (g + h)) / (two * h * y)
         g = sqrt(f*f+one)
         f = ((x - z) * (x + z) + h * (y / (f + sign(g,f)) - h)) / x
c     .......... next qr transformation ..........
         c = one
         s = one
c
         do 600 i1 = l, k1
            i = i1 + 1
            g = rv1(i)
            y = w(i)
            h = s * g
            g = c * g
            z = sqrt(f*f+h*h)
            rv1(i1) = z
            c = f / z
            s = h / z
            f = x * c + g * s
            g = -x * s + g * c
            h = y * s
            y = y * c
            if (.not. matv) go to 575
c
            do 570 j = 1, n
               x = v(j,i1)
               z = v(j,i)
               v(j,i1) = x * c + z * s
               v(j,i) = -x * s + z * c
  570       continue
c
  575       z = sqrt(f*f+h*h)
            w(i1) = z
c     .......... rotation can be arbitrary if z is zero ..........
            if (z .eq. zero) go to 580
            c = f / z
            s = h / z
  580       f = c * g + s * y
            x = -s * g + c * y
            if (.not. matu) go to 600
c
            do 590 j = 1, m
               y = u(j,i1)
               z = u(j,i)
               u(j,i1) = y * c + z * s
               u(j,i) = -y * s + z * c
  590       continue
c
  600    continue
c
         rv1(l) = zero
         rv1(k) = f
         w(k) = x
         go to 520
c     .......... convergence ..........
  650    if (z .ge. zero) go to 700
c     .......... w(k) is made non-negative ..........
         w(k) = -z
         if (.not. matv) go to 700
c
         do 690 j = 1, n
         v(j,k) = -v(j,k)
  690    CONTINUE
c
  700 continue
c
      go to 1001
c     .......... set error -- no convergence to a
c                singular value after 30 iterations ..........
 1000 iarr = k
 1001 return
      end