c**************************************************************

      SUBROUTINE svd_cmp(a,m,n,w,v)

c Routine PUBLIC du module mod_numerique, décomposition A = U W Vt d'une matrice
c rectangulaire, voir "numerical receipes" p. 59

c ATTENTION en sortie V et non Vt (transposée), pour vérification de la décomposition
c il faut former A = U W TRANSPOSE(V)

c Pseudo inverse: A^(-1)(n,m) = V Diag(1/W) Ut

c Auteur: P.Morel, OCA, Observatoire de Nice, CESAM2k

c entrées/sorties:
c	a(m,n): matrice a décomposer, dimensionnée a(mp,np) dans le programme appelant
c	remplacée par U en sortie

c sorties:
c     w: matrice diagonale des valeurs singulières
c     v: matrice V et non V^T

c---------------------------------------------------------------------

      USE mod_kind

      IMPLICIT NONE

      INTEGER, INTENT(in) :: m, n
      REAL (kind=dp), INTENT(inout), DIMENSION(m,n) :: a
      REAL (kind=dp), INTENT(inout), DIMENSION(n,n) :: v
      REAL (kind=dp), INTENT(inout), DIMENSION(n) :: w

      REAL (kind=dp), DIMENSION(n):: rv1

      REAL (kind=dp):: g, scale, anorm, s, f, h, c, x, y, z

      INTEGER i, l, k, j, jj, its, nm

c---------------------------------------------------------------------

2000  FORMAT(8es10.3)

      g=0.0d0 ; scale=0.0d0 ; anorm=0.0d0
      DO i=1,n	!25
       l=i+1 ; rv1(i)=scale*g ; g=0.0d0 ; s=0.0d0 ; scale=0.0d0
       IF(i <= m) THEN
        DO k=i,m	!11
         scale=scale+ABS(a(k,i))
        ENDDO !k 11
        IF(scale /= 0.d0)THEN
         DO k=i,m	!12
          a(k,i)=a(k,i)/scale ; s=s+a(k,i)*a(k,i)
         ENDDO    !k 12
         f=a(i,i) ; g=-SIGN(SQRT(s),f) ; h=f*g-s ; a(i,i)=f-g
         DO j=l,n	!15
          s=0.0d0
          DO k=i,m	!13
           s=s+a(k,i)*a(k,j)
          ENDDO  !k 13
          f=s/h
          DO k=i,m	!14
           a(k,j)=a(k,j)+f*a(k,i)
          ENDDO  !k 14
         ENDDO   !j 15
         DO k=i,m	!16
          a(k,i)=scale*a(k,i)
         ENDDO    !k 16
        ENDIF
       ENDIF
       w(i)=scale*g ; g=0.0d0 ; s=0.0d0 ; scale=0.0d0
       IF((i <= m) .AND. (i /= n))THEN
        DO k=l,n	!17
         scale=scale+ABS(a(i,k))
        ENDDO !k 17
        IF(scale /= 0.0d0)THEN
         DO k=l,n	!18
          a(i,k)=a(i,k)/scale ; s=s+a(i,k)*a(i,k)
         ENDDO    !k 18
         f=a(i,l) ; g=-SIGN(SQRT(s),f) ; h=f*g-s ; a(i,l)=f-g
         DO k=l,n	!19
          rv1(k)=a(i,k)/h
         ENDDO    !k 19
         DO j=l,m	!23
          s=0.0d0
          DO k=l,n	!21
           s=s+a(j,k)*a(i,k)
          ENDDO  !k 21
          DO k=l,n	!22
           a(j,k)=a(j,k)+s*rv1(k)
          ENDDO  !k 22
         ENDDO   !j 23
         DO k=l,n	!24
          a(i,k)=scale*a(i,k)
         ENDDO    !k 24
        ENDIF
       ENDIF
       anorm=MAX(anorm,(ABS(w(i))+ABS(rv1(i))))
      ENDDO   !i 25

      DO i=n,1,-1	!32
       IF(i < n)THEN
        IF(g /= 0.0d0)THEN
         DO j=l,n	!26
          v(j,i)=(a(i,j)/a(i,l))/g
         ENDDO    !j 26
         DO j=l,n	!29
          s=0.0d0
          DO k=l,n	!27
           s=s+a(i,k)*v(k,j)
          ENDDO   !k 27
          DO k=l,n	!28
           v(k,j)=v(k,j)+s*v(k,i)
          ENDDO   !k 28
         ENDDO    !j 29
        ENDIF
        DO j=l,n	!31
         v(i,j)=0.0d0 ; v(j,i)=0.0d0
        ENDDO !j 31
       ENDIF
       v(i,i)=1.0d0 ; g=rv1(i) ; l=i
      ENDDO   !i 32

      DO i=MIN(m,n),1,-1	!39
       l=i+1 ; g=w(i)
       DO j=l,n	!33
        a(i,j)=0.0d0
       ENDDO !j 33
       IF(g /= 0.0d0)THEN
        g=1.0d0/g
        DO j=l,n	!36
         s=0.0d0
         DO k=l,m	!34
          s=s+a(k,i)*a(k,j)
         ENDDO   !k 34
         f=s/a(i,i)*g
         DO k=i,m	!35
          a(k,j)=a(k,j)+f*a(k,i)
         ENDDO   !k 35
        ENDDO    !j 36
        DO j=i,m	!37
         a(j,i)=a(j,i)*g
        ENDDO !j 37
       ELSE
        DO j=i,m	!38
         a(j,i)=0.0d0
        ENDDO !j 38
       ENDIF
       a(i,i)=a(i,i)+1.0d0
      ENDDO   !i 39

      DO k=n,1,-1	!49
       DO its=1,30	!48
        DO l=k,1,-1	!41
         nm=l-1
         IF(ABS(rv1(l))+anorm == anorm)GOTO 2
         IF(ABS(w(nm))+anorm== anorm)GOTO 1
        ENDDO !l 41
1       c=0.0d0 ; s=1.0d0
        DO i=l,k	!43
         f=s*rv1(i)
         IF(ABS(f)+anorm == anorm)GOTO 2
         g=w(i) ; h=SQRT(f*f+g*g) ; w(i)=h
         h=1.0d0/h ; c=g*h ; s=-f*h
         DO j=1,m	!42
          y=a(j,nm) ; z=a(j,i) ; a(j,nm)=y*c+z*s ; a(j,i)=-y*s+z*c
         ENDDO   !j 42
        ENDDO !i 43
2       z=w(k)
        IF(l == k)THEN
         IF(z < 0.0d0)THEN
          w(k)=-z
          DO j=1,n	!44
           v(j,k)=-v(j,k)
          ENDDO   !j 44
         ENDIF
         GOTO 3
        ENDIF
        IF(its == 30)THEN
         WRITE(*,*)'dans svd_cmp, pas de convergence apres 30 iterations'
         STOP
        ENDIF
        x=w(l) ; nm=k-1 ; y=w(nm) ; g=rv1(nm) ; h=rv1(k)
        f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y) ; g=SQRT(f*f+1.0)
        f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x
        c=1.0d0 ; s=1.0d0
        DO j=l,nm	!47
         i=j+1 ; g=rv1(i) ; y=w(i) ; h=s*g ; g=c*g ; z=SQRT(f*f+h*h)
         rv1(j)=z ; c=f/z ; s=h/z ; f=x*c+g*s ; g=-x*s+g*c
         h=y*s ; y=y*c
         DO jj=1,n	!45
          x=v(jj,j) ; z=v(jj,i) ; v(jj,j)=x*c+z*s ; v(jj,i)=-x*s+z*c
         ENDDO    !jj 45
         z=SQRT(f*f+h*h) ; w(j)=z
         IF(z /= 0.d0)THEN
          z=1.0d0/z ; c=f*z ; s=h*z
         ENDIF
         f=c*g+s*y ; x=-s*g+c*y
         DO jj=1,m	!46
          y=a(jj,j) ; z=a(jj,i) ; a(jj,j)=y*c+z*s ; a(jj,i)=-y*s+z*c
         ENDDO    !jj 46
        ENDDO !j 47
        rv1(l)=0.0d0 ; rv1(k)=f ; w(k)=x
       ENDDO  !its 48
3      CONTINUE
      ENDDO   !k 49

      RETURN

      END SUBROUTINE svd_cmp