Numérotation des lignes :

C INVSVD    SOURCE    CHAT      06/03/29    21:23:32     5360      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)cc     Adaptation to Esopec     A. BECCANTINIc     DRN/DMT/SEMT/LTMFc     18.08.00cc     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).cc     this subroutine determines the singular value decompositionc          tc     a=usv  of a real m by n rectangular matrix.  householderc     bidiagonalization and a variant of the qr algorithm are used.cc     on input.cc        nm must be set to the row dimension of two-dimensionalc          array parameters as declared in the calling programc          dimension statement.  note that nm must be at leastc          as large as the maximum of m and n.cc        m is the number of rows of a (and u).cc        n is the number of columns of a (and u) and the order of v.cc        a contains the rectangular input matrix to be decomposed.cc        matu should be set to .true. if the u matrix in thec          decomposition is desired, and to .false. otherwise.cc        matv should be set to .true. if the v matrix in thec          decomposition is desired, and to .false. otherwise.cc     on output.cc        a is unaltered (unless overwritten by u or v).cc        w contains the n (non-negative) singular values of a (thec          diagonal elements of s).  they are unordered.  if anc          error exit is made, the singular values should be correctc          for indices iarr+1,iarr+2,...,n.cc        u contains the matrix u (orthogonal column vectors) of thec          decomposition if matu has been set to .true.  otherwisec          u is used as a temporary array.  u may coincide with a.c          if an error exit is made, the columns of u correspondingc          to indices of correct singular values should be correct.cc        v contains the matrix v (orthogonal) of the decomposition ifc          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 errorc          exit is made, the columns of v corresponding to indices ofc          correct singular values should be correct.cc        iarr is set toc          zero       for normal return,c          k          if the k-th singular value has not beenc                     determined after 30 iterations.cc        rv1 is a temporary storage array.cc     this is a modified version of a routine from the eispackc     collection by the nats projectcc     modified to eliminate machepc      iarr = 0c      do 100 i = 1, mc         do 100 j = 1, n            u(i,j) = a(i,j)  100 continuec     .......... householder reduction to bidiagonal form ..........      g = zero      scale = zero      anorm = zeroc      do 300 i = 1, n         l = i + 1         rv1(i) = scale * g         g = zero         s = zero         scale = zero         if (i .gt. m) go to 210c         do 120 k = i, m  120    scale = scale + abs(u(k,i))c         if (scale .eq. zero) go to 210c         do 130 k = i, m            u(k,i) = u(k,i) / scale            s = s + u(k,i)**2  130    continuec         f = u(i,i)         g = -sign(sqrt(s),f)         h = f * g - s         u(i,i) = f - g         if (i .eq. n) go to 190c         do 150 j = l, n            s = zeroc            do 140 k = i, m  140       s = s + u(k,i) * u(k,j)c            f = s / hc            do 150 k = i, m               u(k,j) = u(k,j) + f * u(k,i)  150    continuec  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 290c         do 220 k = l, n  220    scale = scale + abs(u(i,k))c         if (scale .eq. zero) go to 290c         do 230 k = l, n            u(i,k) = u(i,k) / scale            s = s + u(i,k)**2  230    continuec         f = u(i,l)         g = -sign(sqrt(s),f)         h = f * g - s         u(i,l) = f - gc         do 240 k = l, n  240    rv1(k) = u(i,k) / hc         if (i .eq. m) go to 270c         do 260 j = l, m            s = zeroc            do 250 k = l, n  250       s = s + u(j,k) * u(i,k)c            do 260 k = l, n               u(j,k) = u(j,k) + s * rv1(k)  260    continuec  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 continuec     .......... accumulation of right-hand transformations ..........      if (.not. matv) go to 410c     .......... 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 360c         do 320 j = l, nc     .......... double division avoids possible underflow ..........  320    v(j,i) = (u(i,j) / u(i,l)) / gc         do 350 j = l, n            s = zeroc            do 340 k = l, n  340       s = s + u(i,k) * v(k,j)c            do 350 k = l, n               v(k,j) = v(k,j) + s * v(k,i)  350    continuec  360    do 380 j = l, n            v(i,j) = zero            v(j,i) = zero  380    continuec  390    v(i,i) = one         g = rv1(i)         l = i  400 continuec     .......... accumulation of left-hand transformations ..........  410 if (.not. matu) go to 510c     ..........for i=min(m,n) step -1 until 1 do -- ..........      mn = n      if (m .lt. n) mn = mc      do 500 ii = 1, mn         i = mn + 1 - ii         l = i + 1         g = w(i)         if (i .eq. n) go to 430c         do 420 j = l, n  420    u(i,j) = zeroc  430    if (g .eq. zero) go to 475         if (i .eq. mn) go to 460c         do 450 j = l, n            s = zeroc            do 440 k = l, m  440       s = s + u(k,i) * u(k,j)c     .......... double division avoids possible underflow ..........            f = (s / u(i,i)) / gc            do 450 k = i, m               u(k,j) = u(k,j) + f * u(k,i)  450    continuec  460    do 470 j = i, m  470    u(j,i) = u(j,i) / gc         go to 490c  475    do 480 j = i, m  480    u(j,i) = zeroc  490    u(i,i) = u(i,i) + one  500 continuec     .......... 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 = 0c     .......... 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 565c     .......... rv1(1) is always zero, so there is no exitc                through the bottom of the loop ..........            if (abs(w(l1)) + anorm .eq. anorm) go to 540  530    continuec     .......... cancellation of rv1(l) if l greater than 1 ..........  540    c = zero         s = onec         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 560c            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       continuec  560    continuec     .......... test for convergence ..........  565    z = w(k)         if (l .eq. k) go to 650c     .......... 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)) / xc     .......... next qr transformation ..........         c = one         s = onec         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 575c            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       continuec  575       z = sqrt(f*f+h*h)            w(i1) = zc     .......... 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 600c            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       continuec  600    continuec         rv1(l) = zero         rv1(k) = f         w(k) = x         go to 520c     .......... convergence ..........  650    if (z .ge. zero) go to 700c     .......... w(k) is made non-negative ..........         w(k) = -z         if (.not. matv) go to 700c         do 690 j = 1, n  690    v(j,k) = -v(j,k)c  700 continuec      go to 1001c     .......... set error -- no convergence to ac                singular value after 30 iterations .......... 1000 iarr = k 1001 return      end

© Cast3M 2003 - Tous droits réservés.
Mentions légales