
      SUBROUTINE SVDCMP(A,M,N,Mp,Np,W,V,RV1)

      IMPLICIT NONE

      INTEGER M, N, Mp, Np
      DOUBLE PRECISION A(Mp, Np), W(Np), V(Np,Np), RV1(Np)

c Routine to do singular value decomposition of the matrix A.
c Based on the Press etal routine

      DOUBLE PRECISION anorm, c, f, g, h, s, scale, x, y, z
      INTEGER i, its, j, k, l, nm, jj, l1

      g = 0.0d0
      scale = 0.0d0
      anorm = 0.0d0
      DO i = 1 , N
         l = i + 1
         RV1(i) = scale*g
         g = 0.0d0
         s = 0.0d0
         scale = 0.0d0
         IF ( i.LE.M ) THEN
            DO k = i , M
               scale = scale + ABS(A(k,i))
            ENDDO
            IF ( scale .NE. 0.0d0 ) THEN
               DO k = i , M
                  A(k,i) = A(k,i)/scale
                  s = s + A(k,i)*A(k,i)
               ENDDO
               f = A(i,i)
               g = -SIGN(SQRT(s),f)
               h = f*g - s
               A(i,i) = f - g
               IF ( i.NE.N ) THEN
                  DO j = l , N
                     s = 0.0d0
                     DO k = i , M
                        s = s + A(k,i)*A(k,j)
                     ENDDO
                     f = s/h
                     DO k = i , M
                        A(k,j) = A(k,j) + f*A(k,i)
                     ENDDO
                  ENDDO
               ENDIF
               DO k = i , M
                  A(k,i) = scale*A(k,i)
               ENDDO
            ENDIF
         ENDIF
         W(i) = scale*g
         g = 0.0d0
         s = 0.0d0
         scale = 0.0d0
         IF ( (i.LE.M) .AND. (i.NE.N) ) THEN
            DO k = l , N
               scale = scale + ABS(A(i,k))
            ENDDO
            IF ( scale .NE. 0.0d0 ) THEN
               DO k = l , N
                  A(i,k) = A(i,k)/scale
                  s = s + A(i,k)*A(i,k)
               ENDDO
               f = A(i,l)
               g = -SIGN(SQRT(s),f)
               h = f*g - s
               A(i,l) = f - g
               DO k = l , N
                  RV1(k) = A(i,k)/h
               ENDDO
               IF ( i.NE.M ) THEN
                  DO j = l , M
                     s = 0.0d0
                     DO k = l , N
                        s = s + A(j,k)*A(i,k)
                     ENDDO
                     DO k = l , N
                        A(j,k) = A(j,k) + s*RV1(k)
                     ENDDO
                  ENDDO
               ENDIF
               DO k = l , N
                  A(i,k) = scale*A(i,k)
               ENDDO
            ENDIF
         ENDIF
         anorm = MAX(anorm,(ABS(W(i))+ABS(RV1(i))))
      ENDDO

      DO i = N , 1 , -1
         IF ( i.LT.N ) THEN
            IF ( g .NE. 0.0d0 ) THEN
               DO j = l , N
                  V(j,i) = (A(i,j)/A(i,l))/g
               ENDDO
               DO j = l , N
                  s = 0.0d0
                  DO k = l , N
                     s = s + A(i,k)*V(k,j)
                  ENDDO
                  DO k = l , N
                     V(k,j) = V(k,j) + s*V(k,i)
                  ENDDO
               ENDDO
            ENDIF
            DO j = l , N
               V(i,j) = 0.0d0
               V(j,i) = 0.0d0
            ENDDO
         ENDIF
         V(i,i) = 1.0
         g = RV1(i)
         l = i
      ENDDO


      DO i = N , 1 , -1
         l = i + 1
         g = W(i)
         IF ( i.LT.N ) THEN
            DO j = l , N
               A(i,j) = 0.0d0
            ENDDO
         ENDIF
         IF ( g .NE. 0.0d0 ) THEN
            g = 1.0/g
            IF ( i.NE.N ) THEN
               DO j = l , N
                  s = 0.0d0
                  DO k = l , M
                     s = s + A(k,i)*A(k,j)
                  ENDDO
                  f = (s/A(i,i))*g
                  DO k = i , M
                     A(k,j) = A(k,j) + f*A(k,i)
                  ENDDO
               ENDDO
            ENDIF
            DO j = i , M
               A(j,i) = A(j,i)*g
            ENDDO
         ELSE
            DO j = i , M
               A(j,i) = 0.0d0
            ENDDO
         ENDIF
         A(i,i) = A(i,i) + 1.0
      ENDDO


      DO k = N , 1 , -1
         DO its = 1 , 30
            DO l = k , 1 , -1
               l1 = l
               nm = l - 1
               IF ( (ABS(RV1(l))+anorm).EQ.anorm ) GOTO 380
               IF ( (ABS(W(nm))+anorm).EQ.anorm ) GOTO 340
            ENDDO
 340        CONTINUE
            l = l1
            c = 0.0d0
            s = 1.0d0

            IF ( nm .EQ. 0 ) THEN
               CALL xwrite(
     & 'Warning: problems in SVDCMP likely due to NaN input',2)
            ENDIF

            DO i = l , k
               f = s*RV1(i)
               RV1(i) = c*RV1(i)
               IF ( (ABS(f)+anorm).NE.anorm ) THEN
                  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
                     y = A(j,nm)
                     z = A(j,i)
                     A(j,nm) = (y*c) + (z*s)
                     A(j,i) = -(y*s) + (z*c)
                  ENDDO
               ENDIF
 360        ENDDO
 380        CONTINUE
            z = W(k)
            IF ( l.EQ.k ) THEN
               IF ( z .LT. 0.0d0 ) THEN
                  W(k) = -z
                  DO j = 1 , N
                     V(j,k) = -V(j,k)
                  ENDDO
               ENDIF
               GOTO 500
            ENDIF
            IF ( its .EQ. 30 ) CALL xwrite(
     &               'SVDCMP: No convergence in 30 iterations', 10)
            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.0d0*h*y)
            g = SQRT(f*f+1.0d0)
            f = ((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x
            c = 1.0d0
            s = 1.0d0
            DO j = l , nm
               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
                  x = V(jj,j)
                  z = V(jj,i)
                  V(jj,j) = (x*c) + (z*s)
                  V(jj,i) = -(x*s) + (z*c)
               ENDDO
               z = SQRT(f*f+h*h)
               W(j) = z
               IF ( z .NE. 0.0d0 ) 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
                  y = A(jj,j)
                  z = A(jj,i)
                  A(jj,j) = (y*c) + (z*s)
                  A(jj,i) = -(y*s) + (z*c)
               ENDDO
            ENDDO
            RV1(l) = 0.0d0
            RV1(k) = f
            W(k) = x
         ENDDO
 500     CONTINUE
      ENDDO

      RETURN
      END
