SUBROUTINE CJASVD(M,N,A,RANK,SIGMA,V,INFO)
  
  use omp_lib
  IMPLICIT NONE
  
  INTEGER M, N, INFO, IINFO, FLAG, RANK
  INTEGER IPIV1(N),IPIV2(N)
  REAL SIGMA(N)
  COMPLEX A(M,N),V(N,N),WORK(M),SC(N)
  REAL CS(N),SN(N)
  
  REAL SZERO, SONE, SHALF
  PARAMETER (SONE = 1.0E0, SZERO = 0.0E0, SHALF = 0.5E0)
  COMPLEX CZERO, CONE, CHALF
  PARAMETER (CONE = (1.0E0,0.0E0), CZERO = (0.0E0,0.0E0), &
       CHALF = (0.5E0,0.0E0))
  
  INTEGER I, J, K, L, P, Q, R, ITER, MAXITER
  PARAMETER (MAXITER = 100)
  
  REAL TOL
  REAL TMP1, TMP2, TMP4, TMP5
  REAL TMP10, TMP11, TMP12, TMP13, TMP14, TMP15
  COMPLEX TMP6, TMP8, TMP13C, TMP14C, W, B
  
  REAL H1, G1, R1, X, YR, YC, ZR, ZC, S1, S2
  REAL C, S, T
  REAL EPS, SAFMIN, SCALE1, SIGMX
  
  EXTERNAL SFMA0, SLAMCH, CDOTC, SCNRM2
  REAL SFMA0, SLAMCH, SCNRM2
  COMPLEX CDOTC

  EPS = SLAMCH('P')
  SAFMIN = SLAMCH('S')
!  TOL = SQRT(SAFMIN)*SQRT(EPS)*SHALF
  TOL = SQRT(SAFMIN)
 
  SCALE1 = SONE/REAL(M)/SQRT(SAFMIN)
  SIGMX = MAXVAL(ABS(A(1:M,1:N)))
  IF (SIGMX .GT. SZERO) THEN
     CALL CLASCL( 'G', 0, 0, SIGMX, SCALE1, M, N, &
          A, M, IINFO )
  ENDIF
  
  DO I = 1,N
     SIGMA(I) = SCNRM2(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. SZERO) 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 CSWAP(M,A(1,J),1,A(1,R),1)
           CALL SSWAP(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. SZERO) THEN
           Q = 0
           GO TO 40
        ENDIF

        IF (SIGMA(J) .LT. SAFMIN) THEN
           TMP1 = SONE/SAFMIN
        ELSE
           TMP1 = SONE/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. SZERO) THEN
              SN(P)=SZERO
              SC(P)=CONE
              CYCLE
           ENDIF
           TMP8 = CDOTC(M,WORK(1),1,A(1,P),1)
           GO TO 20
        ENDDO
        Q = N
        GO TO 40

20      K = P

        TMP4 = REAL(TMP8)
        TMP5 = AIMAG(TMP8)
        CALL SLARTG (ABS(TMP4),ABS(TMP5),C,S,TMP2)
        IF (TMP2 .LE. EPS*SIGMA(K)*S1) THEN
           SN(K)=SZERO
           SC(K)=CONE
           GO TO 10
        ENDIF

        C = SIGN(C,TMP4)
        S = SIGN(S,TMP5)
        TMP6 = CMPLX(C,-S)
        
        S2 = SIGMA(K)*TMP1
        H1 = (SIGMA(J) - SIGMA(K))*((S1+S2)*SHALF)
        G1 = TMP2

        R1 = SHALF*ATAN2(G1,H1)
        C = COS(R1)
        S = SIN(R1)
        T = TAN(R1)

        TMP10 = SIGMA(J)*SQRT(SFMA0(T,TMP2/(SIGMA(J)*S1),SONE))
        IF (TMP10 .LT. SAFMIN) THEN
           TMP10 = SAFMIN
        ENDIF
        TMP1 = SONE/TMP10
        
        TMP15 = SFMA0(-T,TMP2/(SIGMA(K)*S2),SONE)
        IF (TMP15 .GT. SZERO) THEN
           TMP11 = SIGMA(K)*SQRT(TMP15)
           IF (TMP11 .LT. SAFMIN) THEN
              TMP11 = SAFMIN
           ENDIF
           TMP12 = SONE/TMP11
        ELSE
           TMP11 = SZERO
        ENDIF

        CS(K)=C
        SN(K)=S
        SC(K)=TMP6
        
        IF (S .NE. SZERO) THEN
           
           X = S/(SONE+C)
                   
           IPIV2(J) = 1
           IPIV2(K) = 1

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

              IF (TMP11 .GT. SZERO) THEN
                 TMP13C = CZERO
                 TMP8 = CZERO
                 TMP14C = CZERO
                 DO L = 1,M
                    YR = REAL(A(L,J))
                    YC = AIMAG(A(L,J))
                    W = TMP6*A(L,K)
                    ZR = REAL(W)
                    ZC = AIMAG(W)
                    A(L,J) = cmplx(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC)
                    WORK(L) = A(L,J)*TMP1
                    TMP13C = TMP13C+CONJG(WORK(L))*WORK(L)
                    TMP8 = TMP8 + CONJG(WORK(L))*A(L,P)
                    A(L,K) = cmplx(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC)
                    B = A(L,K)*TMP12
                    TMP14C = TMP14C+CONJG(B)*B
                 ENDDO
                 TMP13 = REAL(TMP13C)
                 TMP14 = REAL(TMP14C)
                 IF (ABS(T) .LE. TOL .AND. SIGMA(K) .LE. TOL*SIGMA(J)) THEN
                    SIGMA(K) = SZERO
                 ELSE
                    SIGMA(K) = SQRT(TMP14)/TMP12
                 ENDIF
                 S1 = SQRT(TMP13)
                 SIGMA(J) = S1/TMP1
              ELSE
                 TMP13C = CZERO
                 TMP8 = CZERO
                 DO L = 1,M
                    YR = REAL(A(L,J))
                    YC = AIMAG(A(L,J))
                    W = TMP6*A(L,K)
                    ZR = REAL(W)
                    ZC = AIMAG(W)
                    A(L,J) = cmplx(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC)
                    WORK(L) = A(L,J)*TMP1
                    TMP13C = TMP13C+CONJG(WORK(L))*WORK(L)
                    TMP8 = TMP8 + CONJG(WORK(L))*A(L,P)
                    A(L,K) = cmplx(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC)
                 ENDDO
                 TMP13 = REAL(TMP13C)
                 IF (ABS(T) .LE. TOL .AND. SIGMA(K) .LE. TOL*SIGMA(J)) THEN
                    SIGMA(K) = SZERO
                 ELSE
                    SIGMA(K) = SCNRM2(M,A(1,K),1)
                 ENDIF
                 S1 = SQRT(TMP13)
                 SIGMA(J) = S1/TMP1
              ENDIF

              GO TO 30
           ENDDO

           IF (TMP11 .GT. SZERO) THEN
              TMP13C = CZERO
              TMP14C = CZERO
              DO L = 1,M
                 YR = REAL(A(L,J))
                 YC = AIMAG(A(L,J))
                 W = TMP6*A(L,K)
                 ZR = REAL(W)
                 ZC = AIMAG(W)
                 A(L,J) = cmplx(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC)
                 WORK(L) = A(L,J)*TMP1
                 TMP13C = TMP13C+CONJG(WORK(L))*WORK(L)
                 A(L,K) = cmplx(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC)
                 B = A(L,K)*TMP12
                 TMP14C = TMP14C+CONJG(B)*B
              ENDDO
              TMP13 = REAL(TMP13C)
              TMP14 = REAL(TMP14C)
              IF (ABS(T) .LE. TOL .AND. SIGMA(K) .LE. TOL*SIGMA(J)) THEN
                 SIGMA(K) = SZERO
              ELSE
                 SIGMA(K) = SQRT(TMP14)/TMP12
              ENDIF
              S1 = SQRT(TMP13)
              SIGMA(J) = S1/TMP1
           ELSE
              TMP13C = CZERO
              DO L = 1,M
                 YR = REAL(A(L,J))
                 YC = AIMAG(A(L,J))
                 W = TMP6*A(L,K)
                 ZR = REAL(W)
                 ZC = AIMAG(W)
                 A(L,J) = cmplx(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC)
                 WORK(L) = A(L,J)*TMP1
                 TMP13C = TMP13C+CONJG(WORK(L))*WORK(L)
                 A(L,K) = cmplx(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC)
              ENDDO
              TMP13 = REAL(TMP13C)
              IF (ABS(T) .LE. TOL .AND. SIGMA(K) .LE. TOL*SIGMA(J)) THEN
                 SIGMA(K) = SZERO
              ELSE
                 SIGMA(K) = SCNRM2(M,A(1,K),1)
              ENDIF
              S1 = SQRT(TMP13)
              SIGMA(J) = S1/TMP1
           ENDIF
           
30         FLAG = 1
           
           IF (SIGMA(J) .GT. SZERO .AND. P .LE. N) GO TO 20
           
           Q = K
           
        ELSE IF(TMP6 .NE. CONE) THEN
           
           IPIV2(K) = 1
           
           IF (TMP11 .GT. SZERO) THEN
              TMP14C = CZERO
              DO L = 1,M
                 A(L,K) = TMP6*A(L,K)
                 B = A(L,K)*TMP12
                 TMP14C = TMP14C+CONJG(B)*B
              ENDDO
              TMP14 = REAL(TMP14C)
              SIGMA(K) = SQRT(TMP14)/TMP12
           ELSE
              CALL CSCAL(M,TMP6,A(1,K),1)
              SIGMA(K) = SCNRM2(M,A(1,K),1)
           ENDIF
           
           FLAG = 1
           
           GO TO 10
           
        ELSE
           
           GO TO 10
           
        ENDIF
        
40      IF (R .NE. J) THEN
           CALL CSWAP(N,V(1,J),1,V(1,R),1)
        ENDIF
        DO P=J+1,Q
           S=SN(P)
           TMP6=SC(P)
           
           IF (S .NE. SZERO) THEN
              C=CS(P)
              X = S/(SONE+C)
              DO L = 1,N
                 YR = REAL(V(L,J))
                 YC = AIMAG(V(L,J))
                 W = TMP6*V(L,P)
                 ZR = REAL(W)
                 ZC = AIMAG(W)
                 V(L,J) = cmplx(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC)
                 V(L,P) = cmplx(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC)
              ENDDO
           ELSE IF(TMP6 .NE. CONE) THEN
              CALL CSCAL(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. SZERO) 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 CSWAP(M, A(1,K), 1, A(1,I), 1)
     CALL CSWAP(N, V(1,K), 1, V(1,I), 1)
        
  ENDDO

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

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

  INFO = 0
  RETURN
  
END SUBROUTINE CJASVD

PROGRAM TESTSJS
  
  use omp_lib
  IMPLICIT NONE
  
  INTEGER I, J, K, M, N, INFO, SEED, RANK
  REAL TMP1, TMP2
  COMPLEX E
  COMPLEX CZERO, CONE
  PARAMETER (CZERO = (0.0E0,0.0E0), CONE=(1.0E0,0.0E0))
  REAL SZERO, SONE, STWO, SAFMIN
  PARAMETER (SZERO = 0.0E0, SONE=1.0E0, STWO=2.0E0)
  REAL,ALLOCATABLE :: SIGMA(:)
  COMPLEX,ALLOCATABLE :: A(:,:), B(:,:), WORK(:,:), V(:,:)
  REAL SLAMCH
  EXTERNAL USERRAND, SLAMCH
  DOUBLE PRECISION USERRAND

  EXTERNAL CDOTC
  COMPLEX CDOTC
  
  character*1000 ARG1,ARG2,ARG3

  integer t1, t2, t_rate, t_max, diff

  SAFMIN = SLAMCH('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 = REAL(USERRAND())
        TMP2 = REAL(USERRAND())
        A(J,I) = CMPLX(TMP1,TMP2)
        B(J,I) = A(J,I)
     ENDDO
     DO J=I+1,M
        A(J,I) = CZERO
        B(J,I) = CZERO
     ENDDO
  ENDDO

  DO J=1,N
     DO K=1,N
        V(K,J) = CZERO
     ENDDO
     V(J,J) = CONE
  ENDDO

  call system_clock(t1)     ! 開始時を記録
  
  CALL CJASVD(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)*CONJG(V(J,K))
        ENDDO
     ENDDO
  ENDDO
!$OMP END PARALLEL DO

  E = CZERO
!$OMP PARALLEL DO PRIVATE(I) REDUCTION(+:E)
  DO J = 1, N
     DO I = 1, M
        E = E+CONJG(B(I, J))*B(I, J)
     ENDDO
  ENDDO
!$OMP END PARALLEL DO

  WRITE (*,*) 'ERROR',SQRT(REAL(E))

  CALL CHERK ('U', 'C', RANK, M, SONE, A, M, SZERO, WORK, N)
  DO I=1,RANK
     WORK(I,I) = WORK(I,I) - CONE
  ENDDO

  E = CZERO
!$OMP PARALLEL DO REDUCTION(+:E)
  DO I=2,RANK
     E = E+CDOTC(I-1,WORK(1,I),1,WORK(1,I),1)
  ENDDO
!$OMP END PARALLEL DO
  E = STWO*E
  DO I=1,RANK
     E = E+CONJG(WORK(I,I))*WORK(I,I)
  ENDDO
  
  WRITE (*,*) 'ERROR',SQRT(REAL(E))

  CALL CHERK ('U', 'C', RANK, N, SONE, V, N, SZERO, WORK, N)
  DO I=1,RANK
     WORK(I,I) = WORK(I,I) - CONE
  ENDDO

  E = CZERO
!$OMP PARALLEL DO REDUCTION(+:E)
  DO I=2,RANK
     E = E+CDOTC(I-1,WORK(1,I),1,WORK(1,I),1)
  ENDDO
!$OMP END PARALLEL DO
  E = STWO*E
  DO I=1,RANK
     E = E+CONJG(WORK(I,I))*WORK(I,I)
  ENDDO

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