MODULE SVD
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

! Based upon routines from the NSWC (Naval Surface Warfare Center),
! which were based upon LAPACK routines.

! Code converted using TO_F90 by Alan Miller
! Date: 2003-11-11  Time: 17:50:44


CONTAINS


SUBROUTINE drotg(da, db, dc, ds)
 
!     DESIGNED BY C.L.LAWSON, JPL, 1977 SEPT 08

!     CONSTRUCT THE GIVENS TRANSFORMATION

!         ( DC  DS )
!     G = (        ) ,    DC**2 + DS**2 = 1 ,
!         (-DS  DC )

!     WHICH ZEROS THE SECOND ENTRY OF THE 2-VECTOR  (DA,DB)**T .

!     THE QUANTITY R = (+/-)SQRT(DA**2 + DB**2) OVERWRITES DA IN
!     STORAGE.  THE VALUE OF DB IS OVERWRITTEN BY A VALUE Z WHICH
!     ALLOWS DC AND DS TO BE RECOVERED BY THE FOLLOWING ALGORITHM:
!           IF Z=1  SET  DC=0.D0  AND  DS=1.D0
!           IF DABS(Z) < 1  SET  DC=SQRT(1-Z**2)  AND  DS=Z
!           IF DABS(Z) > 1  SET  DC=1/Z  AND  DS=SQRT(1-DC**2)

!     NORMALLY, THE SUBPROGRAM DROT(N,DX,INCX,DY,INCY,DC,DS) WILL
!     NEXT BE CALLED TO APPLY THE TRANSFORMATION TO A 2 BY N MATRIX.

! ------------------------------------------------------------------

REAL (dp), INTENT(IN OUT)  :: da
REAL (dp), INTENT(IN OUT)  :: db
REAL (dp), INTENT(OUT)     :: dc
REAL (dp), INTENT(OUT)     :: ds

REAL (dp)  :: u, v, r
IF (ABS(da) <= ABS(db)) GO TO 10

! *** HERE ABS(DA) > ABS(DB) ***

u = da + da
v = db / u

!     NOTE THAT U AND R HAVE THE SIGN OF DA

r = SQRT(.25D0 + v**2) * u

!     NOTE THAT DC IS POSITIVE

dc = da / r
ds = v * (dc + dc)
db = ds
da = r
RETURN

! *** HERE ABS(DA) <= ABS(DB) ***

10 IF (db == 0.d0) GO TO 20
u = db + db
v = da / u

!     NOTE THAT U AND R HAVE THE SIGN OF DB
!     (R IS IMMEDIATELY STORED IN DA)

da = SQRT(.25D0 + v**2) * u

!     NOTE THAT DS IS POSITIVE

ds = db / da
dc = v * (ds + ds)
IF (dc == 0.d0) GO TO 15
db = 1.d0 / dc
RETURN
15 db = 1.d0
RETURN

! *** HERE DA = DB = 0.D0 ***

20 dc = 1.d0
ds = 0.d0
RETURN

END SUBROUTINE drotg


SUBROUTINE dswap1 (n, dx, dy)

!     INTERCHANGES TWO VECTORS.
!     USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
!     JACK DONGARRA, LINPACK, 3/11/78.
!     This version is for increments = 1.

INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: dx(*)
REAL (dp), INTENT(IN OUT)  :: dy(*)

REAL (dp)  :: dtemp
INTEGER    :: i, m, mp1

IF(n <= 0) RETURN

!       CODE FOR BOTH INCREMENTS EQUAL TO 1

!       CLEAN-UP LOOP

m = MOD(n,3)
IF( m == 0 ) GO TO 40
DO  i = 1,m
  dtemp = dx(i)
  dx(i) = dy(i)
  dy(i) = dtemp
END DO
IF( n < 3 ) RETURN
40 mp1 = m + 1
DO  i = mp1,n,3
  dtemp = dx(i)
  dx(i) = dy(i)
  dy(i) = dtemp
  dtemp = dx(i + 1)
  dx(i + 1) = dy(i + 1)
  dy(i + 1) = dtemp
  dtemp = dx(i + 2)
  dx(i + 2) = dy(i + 2)
  dy(i + 2) = dtemp
END DO
RETURN
END SUBROUTINE  dswap1


SUBROUTINE  drot1 (n, dx, dy, c, s)

!     APPLIES A PLANE ROTATION.
!     JACK DONGARRA, LINPACK, 3/11/78.
!     This version is for increments = 1.

INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: dx(*)
REAL (dp), INTENT(IN OUT)  :: dy(*)
REAL (dp), INTENT(IN)      :: c
REAL (dp), INTENT(IN)      :: s

REAL (dp)  :: dtemp
INTEGER    :: i

IF(n <= 0) RETURN
!       CODE FOR BOTH INCREMENTS EQUAL TO 1

DO  i = 1,n
  dtemp = c*dx(i) + s*dy(i)
  dy(i) = c*dy(i) - s*dx(i)
  dx(i) = dtemp
END DO
RETURN
END SUBROUTINE  drot1


SUBROUTINE dsvdc(x, n, p, s, e, u, v, job, info)

INTEGER, INTENT(IN)        :: n
INTEGER, INTENT(IN)        :: p
REAL (dp), INTENT(IN OUT)  :: x(:,:)
REAL (dp), INTENT(OUT)     :: s(:)
REAL (dp), INTENT(OUT)     :: e(:)
REAL (dp), INTENT(OUT)     :: u(:,:)
REAL (dp), INTENT(OUT)     :: v(:,:)
INTEGER, INTENT(IN)        :: job
INTEGER, INTENT(OUT)       :: info

!     DSVDC IS A SUBROUTINE TO REDUCE A DOUBLE PRECISION NXP MATRIX X
!     BY ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM.  THE
!     DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X.  THE
!     COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
!     AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.

!     ON ENTRY

!         X         DOUBLE PRECISION(LDX,P), WHERE LDX.GE.N.
!                   X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
!                   DECOMPOSITION IS TO BE COMPUTED.  X IS
!                   DESTROYED BY DSVDC.

!         LDX       INTEGER.
!                   LDX IS THE LEADING DIMENSION OF THE ARRAY X.

!         N         INTEGER.
!                   N IS THE NUMBER OF ROWS OF THE MATRIX X.

!         P         INTEGER.
!                   P IS THE NUMBER OF COLUMNS OF THE MATRIX X.

!         LDU       INTEGER.
!                   LDU IS THE LEADING DIMENSION OF THE ARRAY U.
!                   (SEE BELOW).

!         LDV       INTEGER.
!                   LDV IS THE LEADING DIMENSION OF THE ARRAY V.
!                   (SEE BELOW).

!         JOB       INTEGER.
!                   JOB CONTROLS THE COMPUTATION OF THE SINGULAR
!                   VECTORS.  IT HAS THE DECIMAL EXPANSION AB
!                   WITH THE FOLLOWING MEANING

!                        A.EQ.0    DO NOT COMPUTE THE LEFT SINGULAR VECTORS.
!                        A.EQ.1    RETURN THE N LEFT SINGULAR VECTORS IN U.
!                        A.GE.2    RETURN THE FIRST MIN(N,P) SINGULAR
!                                  VECTORS IN U.
!                        B.EQ.0    DO NOT COMPUTE THE RIGHT SINGULAR VECTORS.
!                        B.EQ.1    RETURN THE RIGHT SINGULAR VECTORS IN V.

!     ON RETURN

!         S         DOUBLE PRECISION(MM), WHERE MM=MIN(N+1,P).
!                   THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE SINGULAR
!                   VALUES OF X ARRANGED IN DESCENDING ORDER OF MAGNITUDE.

!         E         DOUBLE PRECISION(P).
!                   E ORDINARILY CONTAINS ZEROS.  HOWEVER SEE THE
!                   DISCUSSION OF INFO FOR EXCEPTIONS.

!         U         DOUBLE PRECISION(LDU,K), WHERE LDU.GE.N.  IF
!                                   JOBA.EQ.1 THEN K.EQ.N, IF JOBA.GE.2
!                                   THEN K.EQ.MIN(N,P).
!                   U CONTAINS THE MATRIX OF LEFT SINGULAR VECTORS.
!                   U IS NOT REFERENCED IF JOBA.EQ.0.  IF N.LE.P
!                   OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X
!                   IN THE SUBROUTINE CALL.

!         V         DOUBLE PRECISION(LDV,P), WHERE LDV.GE.P.
!                   V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
!                   V IS NOT REFERENCED IF JOB.EQ.0.  IF P.LE.N,
!                   THEN V MAY BE IDENTIFIED WITH X IN THE
!                   SUBROUTINE CALL.

!         INFO      INTEGER.
!                   THE SINGULAR VALUES (AND THEIR CORRESPONDING SINGULAR
!                   VECTORS) S(INFO+1),S(INFO+2),...,S(M) ARE CORRECT
!                   (HERE M=MIN(N,P)).  THUS IF INFO.EQ.0, ALL THE
!                   SINGULAR VALUES AND THEIR VECTORS ARE CORRECT.
!                   IN ANY EVENT, THE MATRIX B = TRANS(U)*X*V IS THE
!                   BIDIAGONAL MATRIX WITH THE ELEMENTS OF S ON ITS DIAGONAL
!                   AND THE ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U)
!                   IS THE TRANSPOSE OF U).  THUS THE SINGULAR VALUES
!                   OF X AND B ARE THE SAME.

!     LINPACK. THIS VERSION DATED 03/19/79 .
!     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.

!     DSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.

!     EXTERNAL DROT
!     BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2,DROTG
!     FORTRAN DABS,DMAX1,MAX0,MIN0,MOD,DSQRT

!     INTERNAL VARIABLES

INTEGER :: iter, j, jobu, k, kase, kk, l, ll, lls, lm1, lp1, ls, lu, m, maxit,  &
    mm, mm1, mp1, nct, nctp1, ncu, nrt, nrtp1
REAL (dp) :: t, work(n)
REAL (dp) :: b, c, cs, el, emm1, f, g, scale, shift, sl, sm, sn,  &
    smm1, t1, test, ztest
LOGICAL :: wantu, wantv

!     SET THE MAXIMUM NUMBER OF ITERATIONS.

maxit = 30

!     DETERMINE WHAT IS TO BE COMPUTED.

wantu = .false.
wantv = .false.
jobu = MOD(job,100)/10
ncu = n
IF (jobu > 1) ncu = MIN(n,p)
IF (jobu /= 0) wantu = .true.
IF (MOD(job,10) /= 0) wantv = .true.

!     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
!     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.

info = 0
nct = MIN(n-1, p)
s(1:nct+1) = 0.0_dp
nrt = MAX(0, MIN(p-2,n))
lu = MAX(nct,nrt)
IF (lu < 1) GO TO 170
DO  l = 1, lu
  lp1 = l + 1
  IF (l > nct) GO TO 20
  
!           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
!           PLACE THE L-TH DIAGONAL IN S(L).
  
  s(l) = SQRT( SUM( x(l:n,l)**2 ) )
  IF (s(l) == 0.0D0) GO TO 10
  IF (x(l,l) /= 0.0D0) s(l) = SIGN(s(l), x(l,l))
  x(l:n,l) = x(l:n,l) / s(l)
  x(l,l) = 1.0D0 + x(l,l)

  10 s(l) = -s(l)

  20 IF (p < lp1) GO TO 50
  DO  j = lp1, p
    IF (l > nct) GO TO 30
    IF (s(l) == 0.0D0) GO TO 30
    
!              APPLY THE TRANSFORMATION.
    
    t = -DOT_PRODUCT(x(l:n,l), x(l:n,j)) / x(l,l)
    x(l:n,j) = x(l:n,j) + t * x(l:n,l)
    
!           PLACE THE L-TH ROW OF X INTO  E FOR THE
!           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
    
    30 e(j) = x(l,j)
  END DO

  50 IF (.NOT.wantu .OR. l > nct) GO TO 70
  
!           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK MULTIPLICATION.
  
  u(l:n,l) = x(l:n,l)

  70 IF (l > nrt) CYCLE
  
!           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
!           L-TH SUPER-DIAGONAL IN E(L).
  
  e(l) = SQRT( SUM( e(lp1:p)**2 ) )
  IF (e(l) == 0.0D0) GO TO 80
  IF (e(lp1) /= 0.0D0) e(l) = SIGN(e(l), e(lp1))
  e(lp1:lp1+p-l-1) = e(lp1:p) / e(l)
  e(lp1) = 1.0D0 + e(lp1)

  80 e(l) = -e(l)
  IF (lp1 > n .OR. e(l) == 0.0D0) GO TO 120
  
!              APPLY THE TRANSFORMATION.
  
  work(lp1:n) = 0.0D0
  DO  j = lp1, p
    work(lp1:lp1+n-l-1) = work(lp1:lp1+n-l-1) + e(j) * x(lp1:lp1+n-l-1,j)
  END DO
  DO  j = lp1, p
    x(lp1:lp1+n-l-1,j) = x(lp1:lp1+n-l-1,j) - (e(j)/e(lp1)) * work(lp1:lp1+n-l-1)
  END DO

  120 IF (.NOT.wantv) CYCLE
  
!              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
!              BACK MULTIPLICATION.
  
  v(lp1:p,l) = e(lp1:p)
END DO

!     SET UP THE FINAL BIDIAGONAL MATRIX OF ORDER M.

170 m = MIN(p,n+1)
nctp1 = nct + 1
nrtp1 = nrt + 1
IF (nct < p) s(nctp1) = x(nctp1,nctp1)
IF (n < m) s(m) = 0.0D0
IF (nrtp1 < m) e(nrtp1) = x(nrtp1,m)
e(m) = 0.0D0

!     IF REQUIRED, GENERATE U.

IF (.NOT.wantu) GO TO 300
IF (ncu < nctp1) GO TO 200
DO  j = nctp1, ncu
  u(1:n,j) = 0.0_dp
  u(j,j) = 1.0_dp
END DO

200 DO  ll = 1, nct
  l = nct - ll + 1
  IF (s(l) == 0.0D0) GO TO 250
  lp1 = l + 1
  IF (ncu < lp1) GO TO 220
  DO  j = lp1, ncu
    t = -DOT_PRODUCT(u(l:n,l), u(l:n,j)) / u(l,l)
    u(l:n,j) = u(l:n,j) + t * u(l:n,l)
  END DO

  220 u(l:n,l) = -u(l:n,l)
  u(l,l) = 1.0D0 + u(l,l)
  lm1 = l - 1
  IF (lm1 < 1) CYCLE
  u(1:lm1,l) = 0.0_dp
  CYCLE

  250 u(1:n,l) = 0.0_dp
  u(l,l) = 1.0_dp
END DO

!     IF IT IS REQUIRED, GENERATE V.

300 IF (.NOT.wantv) GO TO 350
DO  ll = 1, p
  l = p - ll + 1
  lp1 = l + 1
  IF (l > nrt) GO TO 320
  IF (e(l) == 0.0D0) GO TO 320
  DO  j = lp1, p
    t = -DOT_PRODUCT(v(lp1:lp1+p-l-1,l), v(lp1:lp1+p-l-1,j)) / v(lp1,l)
    v(lp1:lp1+p-l-1,j) = v(lp1:lp1+p-l-1,j) + t * v(lp1:lp1+p-l-1,l)
  END DO

  320 v(1:p,l) = 0.0D0
  v(l,l) = 1.0D0
END DO

!     MAIN ITERATION LOOP FOR THE SINGULAR VALUES.

350 mm = m
iter = 0

!        QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.

!     ...EXIT
360 IF (m == 0) GO TO 620

!        IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET FLAG AND RETURN.

IF (iter < maxit) GO TO 370
info = m
!     ......EXIT
GO TO 620

!        THIS SECTION OF THE PROGRAM INSPECTS FOR NEGLIGIBLE ELEMENTS
!        IN THE S AND E ARRAYS.  ON COMPLETION
!        THE VARIABLES KASE AND L ARE SET AS FOLLOWS.

!           KASE = 1     IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L < M
!           KASE = 2     IF S(L) IS NEGLIGIBLE AND L < M
!           KASE = 3     IF E(L-1) IS NEGLIGIBLE, L < M, AND
!                        S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP).
!           KASE = 4     IF E(M-1) IS NEGLIGIBLE (CONVERGENCE).

370 DO  ll = 1, m
  l = m - ll
!        ...EXIT
  IF (l == 0) EXIT
  test = ABS(s(l)) + ABS(s(l+1))
  ztest = test + ABS(e(l))
  IF (ztest /= test) CYCLE
  e(l) = 0.0D0
!        ......EXIT
  EXIT
END DO

IF (l /= m - 1) GO TO 410
kase = 4
GO TO 480

410 lp1 = l + 1
mp1 = m + 1
DO  lls = lp1, mp1
  ls = m - lls + lp1
!           ...EXIT
  IF (ls == l) EXIT
  test = 0.0D0
  IF (ls /= m) test = test + ABS(e(ls))
  IF (ls /= l + 1) test = test + ABS(e(ls-1))
  ztest = test + ABS(s(ls))
  IF (ztest /= test) CYCLE
  s(ls) = 0.0D0
!           ......EXIT
  EXIT
END DO

IF (ls /= l) GO TO 450
kase = 3
GO TO 480

450 IF (ls /= m) GO TO 460
kase = 1
GO TO 480

460 kase = 2
l = ls
480 l = l + 1

!        PERFORM THE TASK INDICATED BY KASE.

SELECT CASE ( kase )
  CASE (    1)
    GO TO 490
  CASE (    2)
    GO TO 520
  CASE (    3)
    GO TO 540
  CASE (    4)
    GO TO 570
END SELECT

!        DEFLATE NEGLIGIBLE S(M).

490 mm1 = m - 1
f = e(m-1)
e(m-1) = 0.0D0
DO  kk = l, mm1
  k = mm1 - kk + l
  t1 = s(k)
  CALL drotg(t1, f, cs, sn)
  s(k) = t1
  IF (k == l) GO TO 500
  f = -sn*e(k-1)
  e(k-1) = cs*e(k-1)

  500 IF (wantv) CALL drot1(p, v(1:,k), v(1:,m), cs, sn)
END DO
GO TO 610

!        SPLIT AT NEGLIGIBLE S(L).

520 f = e(l-1)
e(l-1) = 0.0D0
DO  k = l, m
  t1 = s(k)
  CALL drotg(t1, f, cs, sn)
  s(k) = t1
  f = -sn*e(k)
  e(k) = cs*e(k)
  IF (wantu) CALL drot1(n, u(1:,k), u(1:,l-1), cs, sn)
END DO
GO TO 610

!        PERFORM ONE QR STEP.

!           CALCULATE THE SHIFT.

540 scale = MAX(ABS(s(m)), ABS(s(m-1)), ABS(e(m-1)), ABS(s(l)), ABS(e(l)))
sm = s(m)/scale
smm1 = s(m-1)/scale
emm1 = e(m-1)/scale
sl = s(l)/scale
el = e(l)/scale
b = ((smm1 + sm)*(smm1 - sm) + emm1**2)/2.0D0
c = (sm*emm1)**2
shift = 0.0D0
IF (b == 0.0D0 .AND. c == 0.0D0) GO TO 550
shift = SQRT(b**2+c)
IF (b < 0.0D0) shift = -shift
shift = c/(b + shift)

550 f = (sl + sm)*(sl - sm) - shift
g = sl*el

!           CHASE ZEROS.

mm1 = m - 1
DO  k = l, mm1
  CALL drotg(f, g, cs, sn)
  IF (k /= l) e(k-1) = f
  f = cs*s(k) + sn*e(k)
  e(k) = cs*e(k) - sn*s(k)
  g = sn*s(k+1)
  s(k+1) = cs*s(k+1)
  IF (wantv) CALL drot1(p, v(1:,k), v(1:,k+1), cs, sn)
  CALL drotg(f, g, cs, sn)
  s(k) = f
  f = cs*e(k) + sn*s(k+1)
  s(k+1) = -sn*e(k) + cs*s(k+1)
  g = sn*e(k+1)
  e(k+1) = cs*e(k+1)
  IF (wantu .AND. k < n) CALL drot1(n, u(1:,k), u(1:,k+1), cs, sn)
END DO
e(m-1) = f
iter = iter + 1
GO TO 610

!        CONVERGENCE.

!           MAKE THE SINGULAR VALUE  POSITIVE.

570 IF (s(l) >= 0.0D0) GO TO 590
s(l) = -s(l)
IF (wantv) v(1:p,l) = -v(1:p,l)

!           ORDER THE SINGULAR VALUE.

590 IF (l == mm) GO TO 600
!           ...EXIT
IF (s(l) >= s(l+1)) GO TO 600
t = s(l)
s(l) = s(l+1)
s(l+1) = t
IF (wantv .AND. l < p) CALL dswap1(p, v(1:,l), v(1:,l+1))
IF (wantu .AND. l < n) CALL dswap1(n, u(1:,l), u(1:,l+1))
l = l + 1
GO TO 590

600 iter = 0
m = m - 1

610 GO TO 360

620 RETURN
END SUBROUTINE dsvdc

END MODULE SVD