SUBROUTINE DLARTGU (F1, G1, C, S, R1)

  IMPLICIT NONE

  DOUBLE PRECISION ZERO, ONE, HALF
  PARAMETER (ONE = 1.0D0, ZERO = 0.0D0, HALF = 0.5D0)
  DOUBLE PRECISION F1, G1, C, S, R1, U, V
  EXTERNAL DFMA0
  DOUBLE PRECISION DFMA0

  IF (F1 .GE. G1) THEN
     U = G1/F1
     R1 = SQRT(DFMA0(U,U,ONE))
     C = ONE/R1
     S = U/R1
     R1 = R1*F1
  ELSE
     V = F1/G1
     R1 = SQRT(DFMA0(V,V,ONE))
     S = ONE/R1
     C = V/R1
     R1 = R1*G1
  ENDIF

  RETURN

END SUBROUTINE DLARTGU

SUBROUTINE ZJASVD(M,N,A,RANK,SIGMA,V,INFO)
  
  use omp_lib
  IMPLICIT NONE
  
  INTEGER M, N, INFO, IINFO, FLAG, RANK
  INTEGER IPIV1(N),IPIV2(N)
  DOUBLE PRECISION SIGMA(N)
  COMPLEX*16 A(M,N),V(N,N),WORK(M),SC(N)
  DOUBLE PRECISION CS(N),SN(N)
  
  DOUBLE PRECISION DZERO, DONE, DHALF
  PARAMETER (DONE = 1.0D0, DZERO = 0.0D0, DHALF = 0.5D0)
  COMPLEX*16 ZZERO, ZONE, ZHALF
  PARAMETER (ZONE = (1.0D0,0.0D0), ZZERO = (0.0D0,0.0D0), &
       ZHALF = (0.5D0,0.0D0))
  
  INTEGER I, J, K, L, P, Q, R, ITER, MAXITER
  PARAMETER (MAXITER = 100)
  
  DOUBLE PRECISION T
  DOUBLE PRECISION TMP1, TMP2, TMP4, TMP5
  DOUBLE PRECISION TMP10, TMP11, TMP12, TMP13, TMP14, TMP15
  COMPLEX*16 TMP6, TMP7, TMP8, TMP13C, TMP14C, W, B
  
  DOUBLE PRECISION H1, G1, R1, X, YR, YC, ZR, ZC, S1, S2
  DOUBLE PRECISION C, S, C0, S0
  DOUBLE PRECISION EPS, SAFMIN, SCALE1, SIGMX
  
  EXTERNAL DFMA0, DLAMCH, ZDOTC, DZNRM2
  DOUBLE PRECISION DFMA0, DLAMCH, DZNRM2
  COMPLEX*16 ZDOTC

  EPS = DLAMCH('P')
  SAFMIN = DLAMCH('S')

  SCALE1 = DONE/DBLE(M)/SQRT(SAFMIN)
  SIGMX = MAXVAL(ABS(A(1:M,1:N)))
  IF (SIGMX .GT. DZERO) THEN
     CALL ZLASCL( 'G', 0, 0, SIGMX, SCALE1, M, N, &
          A, M, IINFO )
  ENDIF
  
  DO I = 1,N
     SIGMA(I) = DZNRM2(M,A(1,I),1)
  ENDDO

  DO J = 1, N
     IPIV1(J) = 1
  ENDDO
  
  DO ITER=1, MAXITER

!     WRITE (*,*) 'ITER',ITER
     
     FLAG = 0

     DO J = 1, N
        IPIV2(J) = 0
     ENDDO
     
     DO J = 1, N-1

        IF (IPIV1(J) .EQ. 0 .OR. SIGMA(J) .EQ. DZERO) CYCLE
        
        R = J
        DO L = J+1, N
           IF (IPIV1(L) .EQ. 1 .AND. SIGMA(L) .GT. SIGMA(R)) R = L
        ENDDO

        IF (R .NE. J) THEN
           CALL ZSWAP(M,A(1,J),1,A(1,R),1)
           CALL DSWAP(1,SIGMA(J),1,SIGMA(R),1)
           L = IPIV1(R)
           IPIV1(R) = IPIV1(J) 
           IPIV1(J) = L
           
           L = IPIV2(R)
           IPIV2(R) = IPIV2(J) 
           IPIV2(J) = L
        ENDIF

        IF (IPIV1(J) .EQ. 0 .OR. SIGMA(J) .EQ. DZERO) THEN
           Q = 0
           GO TO 40
        ENDIF

        IF (SIGMA(J) .LT. SAFMIN) THEN
           TMP1 = DONE/SAFMIN
        ELSE
           TMP1 = DONE/SIGMA(J)
        ENDIF
        S1 = SIGMA(J)*TMP1
        DO L = 1,M
           WORK(L) = A(L,J)*TMP1
        ENDDO

        K = J

10      DO P = K+1, N
           IF (IPIV1(P) .EQ. 0 .OR. SIGMA(P) .EQ. DZERO) THEN
              SN(P)=DZERO
              SC(P)=ZONE
              CYCLE
           ENDIF
           TMP8 = ZDOTC(M,WORK(1),1,A(1,P),1)
           GO TO 20
        ENDDO
        Q = N
        GO TO 40

20      K = P

        TMP4 = DBLE(TMP8)
        TMP5 = DIMAG(TMP8)
        CALL DLARTGU (ABS(TMP4),ABS(TMP5),C,S,TMP2)
        IF (TMP2 .LT. SAFMIN) THEN
           IF (TMP2 .LE. SQRT(DBLE(M))*EPS*SIGMA(K)*S1) THEN
              SN(K)=DZERO
              SC(K)=ZONE
              GO TO 10
           ENDIF
        ELSE
           IF (TMP2 .LE. EPS*SIGMA(K)*S1) THEN
              SN(K)=DZERO
              SC(K)=ZONE
              GO TO 10
           ENDIF
        ENDIF

        C = SIGN(C,TMP4)
        S = SIGN(S,TMP5)
        TMP6 = DCMPLX(C,-S)
        
        S2 = SIGMA(K)*TMP1
        H1 = (SIGMA(J) - SIGMA(K))*((S1+S2)*DHALF)
        G1 = TMP2 
        
        CALL DLARTGU (ABS(H1), ABS(G1), C, S, R1)
        T = G1 / (H1 + SIGN(R1,H1))
        R1 = SQRT(DFMA0(T,T,DONE))
        C = DONE/R1
        S = T/R1
        
        TMP10 = SIGMA(J)*SQRT(DFMA0(T,TMP2/(SIGMA(J)*S1),DONE))
        IF (TMP10 .LT. SAFMIN) THEN
           TMP10 = SAFMIN
        ENDIF
        TMP1 = DONE/TMP10
        
        TMP15 = DFMA0(-T,TMP2/(SIGMA(K)*S2),DONE)
        IF (TMP15 .GT. DZERO) THEN
           TMP11 = SIGMA(K)*SQRT(TMP15)
           IF (TMP11 .LT. SAFMIN) THEN
              TMP11 = SAFMIN
           ENDIF
           TMP12 = DONE/TMP11
        ELSE
           TMP11 = DZERO
        ENDIF

        CS(K)=C
        SN(K)=S
        SC(K)=TMP6
        
        IF (S .NE. DZERO) THEN

           X = S/(DONE+C)
        
           IPIV2(J) = 1
           IPIV2(K) = 1

           DO P = K+1, N
              IF (IPIV1(P) .EQ. 0 .OR. SIGMA(P) .EQ. DZERO) THEN
                 SN(P)=DZERO
                 SC(P)=ZONE
                 CYCLE
              ENDIF

              IF (TMP11 .GT. DZERO) THEN
                 TMP13C = ZZERO
                 TMP8 = ZZERO
                 TMP14C = ZZERO
                 DO L = 1,M
                    YR = DBLE(A(L,J))
                    YC = DIMAG(A(L,J))
                    W = TMP6*A(L,K)
                    ZR = DBLE(W)
                    ZC = DIMAG(W)
                    A(L,J) = DCMPLX(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC)
                    WORK(L) = A(L,J)*TMP1
                    TMP13C = TMP13C+DCONJG(WORK(L))*WORK(L)
                    TMP8 = TMP8 + DCONJG(WORK(L))*A(L,P)
                    A(L,K) = DCMPLX(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC)
                    B = A(L,K)*TMP12
                    TMP14C = TMP14C+DCONJG(B)*B
                 ENDDO
                 TMP13 = DBLE(TMP13C)
                 TMP14 = DBLE(TMP14C)
                 SIGMA(J) = TMP10*SQRT(TMP13)
                 S1 = SIGMA(J)*TMP1
                 SIGMA(K) = TMP11*SQRT(TMP14)
              ELSE
                 TMP13C = ZZERO
                 TMP8 = ZZERO
                 DO L = 1,M
                    YR = DBLE(A(L,J))
                    YC = DIMAG(A(L,J))
                    W = TMP6*A(L,K)
                    ZR = DBLE(W)
                    ZC = DIMAG(W)
                    A(L,J) = DCMPLX(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC)
                    WORK(L) = A(L,J)*TMP1
                    TMP13C = TMP13C+DCONJG(WORK(L))*WORK(L)
                    TMP8 = TMP8 + DCONJG(WORK(L))*A(L,P)
                    A(L,K) = DCMPLX(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC)
                 ENDDO
                 TMP13 = DBLE(TMP13C)
                 SIGMA(J) = TMP10*SQRT(TMP13)
                 S1 = SIGMA(J)*TMP1
                 SIGMA(K) = DZNRM2(M,A(1,K),1)
              ENDIF

              GO TO 30
           ENDDO
           
           IF (TMP11 .GT. DZERO) THEN
              TMP13C = ZZERO
              TMP14C = ZZERO
              DO L = 1,M
                 YR = DBLE(A(L,J))
                 YC = DIMAG(A(L,J))
                 W = TMP6*A(L,K)
                 ZR = DBLE(W)
                 ZC = DIMAG(W)
                 A(L,J) = DCMPLX(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC)
                 WORK(L) = A(L,J)*TMP1
                 TMP13C = TMP13C+DCONJG(WORK(L))*WORK(L)
                 A(L,K) = DCMPLX(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC)
                 B = A(L,K)*TMP12
                 TMP14C = TMP14C+DCONJG(B)*B
              ENDDO
              TMP13 = DBLE(TMP13C)
              TMP14 = DBLE(TMP14C)
              SIGMA(J) = TMP10*SQRT(TMP13)
              S1 = SIGMA(J)*TMP1
              SIGMA(K) = TMP11*SQRT(TMP14)
           ELSE
              TMP13C = ZZERO
              DO L = 1,M
                 YR = DBLE(A(L,J))
                 YC = DIMAG(A(L,J))
                 W = TMP6*A(L,K)
                 ZR = DBLE(W)
                 ZC = DIMAG(W)
                 A(L,J) = DCMPLX(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC)
                 WORK(L) = A(L,J)*TMP1
                 TMP13C = TMP13C+DCONJG(WORK(L))*WORK(L)
                 A(L,K) = DCMPLX(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC)
              ENDDO
              TMP13 = DBLE(TMP13C)
              SIGMA(J) = TMP10*SQRT(TMP13)
              S1 = SIGMA(J)*TMP1
              SIGMA(K) = DZNRM2(M,A(1,K),1)
           ENDIF
                 
30         IF (FLAG .EQ. 0) THEN
              TMP7 = ZDOTC(M,WORK(1),1,A(1,K),1)
              TMP13 = DBLE(TMP7)
              TMP14 = DIMAG(TMP7)
              CALL DLARTGU (ABS(TMP13),ABS(TMP14),C0,S0,TMP4)
              IF (TMP4 .LT. SAFMIN) THEN
                 IF (TMP4 .LE. SQRT(DBLE(M))*EPS*SIGMA(K)*S1) FLAG = 1
              ELSE
                 IF (TMP4 .LE. EPS*SIGMA(K)*S1) FLAG = 1
              ENDIF
           ENDIF
           
           IF (SIGMA(J) .GT. DZERO .AND. P .LE. N) GO TO 20

           Q = K
           
        ELSE IF(TMP6 .NE. ZONE) THEN
           
           IPIV2(K) = 1

           IF (TMP11 .GT. DZERO) THEN
              TMP14C = ZZERO
              DO L = 1,M
                 A(L,K) = TMP6*A(L,K)
                 B = A(L,K)*TMP12
                 TMP14C = TMP14C+DCONJG(B)*B
              ENDDO
              TMP14 = DBLE(TMP14C)
              SIGMA(K) = TMP11*SQRT(TMP14)
           ELSE
              CALL ZSCAL(M,TMP6,A(1,K),1)
              SIGMA(K) = DZNRM2(M,A(1,K),1)
           ENDIF
           
           IF (FLAG .EQ. 0) THEN
              TMP7 = ZDOTC(M,WORK(1),1,A(1,K),1)
              TMP13 = DBLE(TMP7)
              TMP14 = DIMAG(TMP7)
              CALL DLARTGU(ABS(TMP13),ABS(TMP14),C0,S0,TMP4)
              IF (TMP4 .LT. SAFMIN) THEN
                 IF (TMP4 .LE. SQRT(DBLE(M))*EPS*SIGMA(K)*S1) FLAG = 1
              ELSE
                 IF (TMP4 .LE. EPS*SIGMA(K)*S1) FLAG = 1
              ENDIF
           ENDIF
           
           GO TO 10

        ELSE
           
           GO TO 10

        ENDIF

40      IF (R .NE. J) THEN
           CALL ZSWAP(N,V(1,J),1,V(1,R),1)
        ENDIF
        DO P=J+1,Q
           S=SN(P)
           TMP6=SC(P)

           IF (S .NE. DZERO) THEN
              C=CS(P)
              X = S/(DONE+C)
              DO L = 1,N
                 YR = DBLE(V(L,J))
                 YC = DIMAG(V(L,J))
                 W = TMP6*V(L,P)
                 ZR = DBLE(W)
                 ZC = DIMAG(W)
                 V(L,J) = dcmplx(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC)
                 V(L,P) = dcmplx(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC)
              ENDDO
           ELSE IF(TMP6 .NE. ZONE) THEN
              CALL ZSCAL(N,TMP6,V(1,P),1)
           ENDIF
        ENDDO
        
     ENDDO

     IF (FLAG .EQ. 0) EXIT

     DO J = 1, N
        IPIV1(J) = IPIV2(J)
     ENDDO
     
  ENDDO
  
  IF (ITER .GT. MAXITER) THEN
     INFO = -1
     RETURN
  ENDIF

  DO I = 1,N
     IF (SIGMA(I) .GT. DZERO) THEN
        DO L = 1, M
           A(L,I) = A(L,I)/SIGMA(I)
        ENDDO
     ENDIF
  ENDDO
  
  DO I = 1, N-1
         
     K = I
     T = SIGMA(K)
     
     DO J = I+1, N
        IF (SIGMA(J) .GT. T) THEN
           K = J
           T = SIGMA(K)
        ENDIF
     ENDDO
        
     SIGMA(K) = SIGMA(I)
     SIGMA(I) = T
        
     CALL ZSWAP(M, A(1,K), 1, A(1,I), 1)
     CALL ZSWAP(N, V(1,K), 1, V(1,I), 1)
        
  ENDDO

  DO I = 1,N
     IF(SIGMA(I) .LT. SAFMIN) THEN
        RANK = I-1
        GO TO 50
     ENDIF
  ENDDO
  RANK = N

50 IF (SIGMX .GT. DZERO) THEN
     CALL DLASCL( 'G', 0, 0, SCALE1, SIGMX, N, 1, SIGMA, N, IINFO )
  ENDIF

  INFO = 0
  RETURN
  
END SUBROUTINE ZJASVD

PROGRAM TESTSJS
  
  use omp_lib
  IMPLICIT NONE
  
  INTEGER I, J, K, M, N, INFO, SEED, RANK
  DOUBLE PRECISION TMP1,TMP2
  COMPLEX*16 E
  COMPLEX*16 ZZERO,ZONE
  PARAMETER (ZZERO = (0.0D0,0.0D0), ZONE=(1.0D0,0.0D0))
  DOUBLE PRECISION DZERO, DONE, DTWO, SAFMIN
  PARAMETER (DZERO = 0.0D0, DONE=1.0D0, DTWO=2.0D0)
  DOUBLE PRECISION,ALLOCATABLE :: SIGMA(:)
  COMPLEX*16,ALLOCATABLE :: A(:,:), B(:,:), WORK(:,:), V(:,:)
  EXTERNAL USERRAND, DLAMCH
  DOUBLE PRECISION USERRAND, DLAMCH

  EXTERNAL ZDOTC
  COMPLEX*16 ZDOTC
  
  character*1000 ARG1,ARG2,ARG3

  integer t1, t2, t_rate, t_max, diff

  SAFMIN = DLAMCH('S')
  
  CALL GETARG(1,ARG1)
  CALL GETARG(2,ARG2)
  CALL GETARG(3,ARG3)

  read (ARG1,*) M
  read (ARG2,*) N
  read (ARG3,*) SEED
  
  ALLOCATE(SIGMA(N))
  ALLOCATE(A(M,N), B(M,N), WORK(N,N), V(N,N))

  DO I=1,SEED
     CALL XOR()
  ENDDO
  
  DO I=1,N
     DO J=1,M
        TMP1 = USERRAND()
        TMP2 = USERRAND()
        A(J,I) =  DCMPLX(TMP1,TMP2)
        B(J,I) =  A(J,I)
     ENDDO
  ENDDO

  DO J=1,N
     DO K=1,N
        V(K,J) = ZZERO
     ENDDO
     V(J,J) = ZONE
  ENDDO
  
  call system_clock(t1)     ! 開始時を記録
  
  CALL ZJASVD(M,N,A,RANK,SIGMA,V,INFO)

  call system_clock(t2, t_rate, t_max) ! 終了時を記録
  if ( t2 < t1 ) then
     diff = (t_max - t1) + t2 + 1
  else
     diff = t2 - t1
  endif

  WRITE (*,*) 'RANK',RANK

  print "(A, F10.3)", "time it took was:", diff/dble(t_rate)
  
  IF (INFO .EQ. -1) THEN
     WRITE (*,*) 'not converge'
     RETURN
  ENDIF

  !$OMP PARALLEL DO PRIVATE(K,I)
  DO J = 1, N
     DO K = 1, RANK
        DO I = 1, M
           B(I, J) = B(I, J) - SIGMA(K)*A(I,K)*DCONJG(V(J,K))
        ENDDO
     ENDDO
  ENDDO
  !$OMP END PARALLEL DO

  E = ZZERO
  !$OMP PARALLEL DO PRIVATE(I) REDUCTION(+:E)
  DO J = 1, N
     DO I = 1, M
        E = E+DCONJG(B(I, J))*B(I, J)
     ENDDO
  ENDDO
  !$OMP END PARALLEL DO
  
  WRITE (*,*) 'ERROR',SQRT(DBLE(E))

  CALL ZHERK('U', 'C', RANK, M, DONE, A, M, DZERO, WORK, N)
  DO I=1,RANK
     WORK(I,I) = WORK(I,I) - ZONE
  ENDDO

  E = ZZERO
  !$OMP PARALLEL DO REDUCTION(+:E)
  DO I=2,RANK
     E = E+ZDOTC(I-1,WORK(1,I),1,WORK(1,I),1)
  ENDDO
  !$OMP END PARALLEL DO
  E = DTWO*E
  DO I=1,RANK
     E = E+DCONJG(WORK(I,I))*WORK(I,I)
  ENDDO
  
  WRITE (*,*) 'ERROR',SQRT(DBLE(E))

  CALL ZHERK('U', 'C', RANK, N, DONE, V, N, DZERO, WORK, N)
  DO I=1,RANK
     WORK(I,I) = WORK(I,I) - ZONE
  ENDDO

  E = ZZERO
  !$OMP PARALLEL DO REDUCTION(+:E)
  DO I=2,RANK
     E = E+ZDOTC(I-1,WORK(1,I),1,WORK(1,I),1)
  ENDDO
  !$OMP END PARALLEL DO
  E = DTWO*E
  DO I=1,RANK
     E = E+DCONJG(WORK(I,I))*WORK(I,I)
  ENDDO

  WRITE (*,*) 'ERROR',SQRT(DBLE(E))
  
END PROGRAM TESTSJS