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 TOL
  DOUBLE PRECISION TMP1, TMP2, TMP4, TMP5
  DOUBLE PRECISION TMP10, TMP11, TMP12, TMP13, TMP14, TMP15
  COMPLEX*16 TMP6, TMP8, TMP13C, TMP14C, W, B
  
  DOUBLE PRECISION H1, G1, R1, X, YR, YC, ZR, ZC, S1, S2
  DOUBLE PRECISION C, S, T
  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')
!  TOL = SQRT(SAFMIN)*SQRT(EPS)*DHALF
  TOL = SQRT(SAFMIN)

  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) .LE. 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) .LE. 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) .LE. 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 DLARTG (ABS(TMP4),ABS(TMP5),C,S,TMP2)
        IF (TMP2 .LE. EPS*SIGMA(K)*S1) THEN
           SN(K)=DZERO
           SC(K)=ZONE
           GO TO 10
        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 
        
        R1 = DHALF*DATAN2(G1,H1)
        C = DCOS(R1)
        S = DSIN(R1)
        T = DTAN(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
           TMP12 = 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) .LE. 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)
                 IF (ABS(T) .LE. TOL .AND. SIGMA(K) .LE. TOL*SIGMA(J)) THEN
                    SIGMA(K) = DZERO
                 ELSE
                    SIGMA(K) = SQRT(TMP14)/TMP12
                 ENDIF
                 S1 = SQRT(TMP13)
                 SIGMA(J) = S1/TMP1
              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)
                 IF (ABS(T) .LE. TOL .AND. SIGMA(K) .LE. TOL*SIGMA(J)) THEN
                    SIGMA(K) = DZERO
                 ELSE
                    SIGMA(K) = DZNRM2(M,A(1,K),1)
                 ENDIF
                 S1 = SQRT(TMP13)
                 SIGMA(J) = S1/TMP1
              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)
              IF (ABS(T) .LE. TOL .AND. SIGMA(K) .LE. TOL*SIGMA(J)) THEN
                 SIGMA(K) = DZERO
              ELSE
                 SIGMA(K) = SQRT(TMP14)/TMP12
              ENDIF
              S1 = SQRT(TMP13)
              SIGMA(J) = S1/TMP1
           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)
              IF (ABS(T) .LE. TOL .AND. SIGMA(K) .LE. TOL*SIGMA(J)) THEN
                 SIGMA(K) = DZERO
              ELSE
                 SIGMA(K) = DZNRM2(M,A(1,K),1)
              ENDIF
              S1 = SQRT(TMP13)
              SIGMA(J) = S1/TMP1
           ENDIF
                 
30         FLAG = 1
           
           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) = SQRT(TMP14)/TMP12
           ELSE
              CALL ZSCAL(M,TMP6,A(1,K),1)
              SIGMA(K) = DZNRM2(M,A(1,K),1)
           ENDIF
           
           FLAG = 1
           
           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) .LE. DZERO) 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,I
        TMP1 = USERRAND()
        TMP2 = USERRAND()
        A(J,I) =  DCMPLX(TMP1,TMP2)
        B(J,I) =  A(J,I)
     ENDDO
     DO J=I+1, M
        A(J,I) =  ZZERO
        B(J,I) =  ZZERO
     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'
     STOP
  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