SUBROUTINE CSROTU(N, X, INCX, Y, INCY, F, S, P)

  IMPLICIT NONE
  
  INTEGER I,N,INCX,INCY
  COMPLEX X(*),Y(*),U,P
  REAL F,S,UR,UC,VR,VC

  IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN
     DO I=1,N
        U=X(I)*P
        UR=REAL(U)
        UC=AIMAG(U)
        VR=REAL(Y(I))
        VC=AIMAG(Y(I))
        X(I)=CMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC)
        Y(I)=CMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC)
     ENDDO
  ELSE
     DO I=1,N
        U=X((I-1)*INCX+1)*P
        UR=REAL(U)
        UC=AIMAG(U)
        VR=REAL(Y((I-1)*INCY+1))
        VC=AIMAG(Y((I-1)*INCY+1))
        X((I-1)*INCX+1)=CMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC)
        Y((I-1)*INCY+1)=CMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC)
     ENDDO
  ENDIF
  
  RETURN
  
END SUBROUTINE CSROTU

SUBROUTINE CSROTV(N, X, INCX, Y, INCY, C, G, P)

  IMPLICIT NONE
  INTEGER I,N,INCX,INCY
  COMPLEX X(*),Y(*),U,P
  REAL G,C,VR,VC,UR,UC

  IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN
     DO I=1,N
        U=X(I)*P
        UR=REAL(U)
        UC=AIMAG(U)
        VR=REAL(Y(I))
        VC=AIMAG(Y(I))
        X(I)=CMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC)
        Y(I)=CMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC)
     ENDDO
  ELSE
     DO I=1,N
        U=X((I-1)*INCX+1)*P
        UR=REAL(U)
        UC=AIMAG(U)
        VR=REAL(Y((I-1)*INCY+1))
        VC=AIMAG(Y((I-1)*INCY+1))
        X((I-1)*INCX+1)=CMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC)
        Y((I-1)*INCY+1)=CMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC)
     ENDDO
  ENDIF
  
  RETURN
  
END SUBROUTINE CSROTV

SUBROUTINE CSROTW(N, X, INCX, Y, INCY, C, G, P)

  IMPLICIT NONE
  INTEGER I,N,INCX,INCY
  COMPLEX X(*),Y(*),U,P
  REAL G,C,VR,VC,UR,UC

  IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN
     DO I=1,N
        U=X(I)*P
        UR=REAL(U)
        UC=AIMAG(U)
        VR=REAL(Y(I))
        VC=AIMAG(Y(I))
        X(I)=CMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC)
        Y(I)=CMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC)
     ENDDO
  ELSE
     DO I=1,N
        U=X((I-1)*INCX+1)*P
        UR=REAL(U)
        UC=AIMAG(U)
        VR=REAL(Y((I-1)*INCY+1))
        VC=AIMAG(Y((I-1)*INCY+1))
        X((I-1)*INCX+1)=CMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC)
        Y((I-1)*INCY+1)=CMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC)
     ENDDO
  ENDIF
  
  RETURN
  
END SUBROUTINE CSROTW

SUBROUTINE CBIDIAG(N,A,U,V,INFO)
  
  use omp_lib
  IMPLICIT NONE
  
  INTEGER N, INFO
  COMPLEX A(N,N),V(N,N),U(N,N)
  
  REAL SONE
  PARAMETER (SONE = 1.0E0)

  COMPLEX CZERO
  PARAMETER (CZERO = (0.0E0,0.0E0))
  INTEGER I, J

  REAL TMP1, TMP2, TMP3, TMP6, TMP13
  COMPLEX TMP7(N),TMP8,TMP9,TMP11,TMP14(N)
  
  REAL SZERO
  PARAMETER (SZERO = 0.0E0)
  
  REAL C1(N), S1(N), C2(N), S2(N), C, S, F

  DO I=1, N-2
     DO J=N,I+2,-1

        IF (J .EQ. N) THEN
           TMP1 = REAL(A(I,J))
           TMP2 = AIMAG(A(I,J))
           CALL SLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP13)
           C = SIGN(C,TMP1)
           S = SIGN(S,TMP2)
           TMP8 = CMPLX(C,-S)
           CALL CSCAL(J-I,TMP8,A(I+1,J),1)
        ENDIF
           
        TMP1 = REAL(A(I,J-1))
        TMP2 = AIMAG(A(I,J-1))
        CALL SLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP3)
        C = SIGN(C,TMP1)
        S = SIGN(S,TMP2)
        TMP7(J) = CMPLX(C,-S)

        CALL SLARTG(TMP3,TMP13,C1(J),S1(J),TMP13)
        A(I,J-1) = TMP13
        A(I,J) = CZERO

        IF (C1(J) .GE. ABS(S1(J))) THEN
           F=S1(J)/(SONE+C1(J))
           CALL CSROTU(J-I,A(I+1,J-1),1,A(I+1,J),1,F,S1(J),TMP7(J))
        ELSE IF (S1(J) .GE. SZERO) THEN
           F=C1(J)/(SONE+S1(J))
           CALL CSROTV(J-I,A(I+1,J-1),1,A(I+1,J),1,C1(J),F,TMP7(J))
        ELSE IF (S1(J) .LE. SZERO) THEN
           F=C1(J)/(SONE-S1(J))
           CALL CSROTW(J-I,A(I+1,J-1),1,A(I+1,J),1,C1(J),F,TMP7(J))
        ENDIF

        IF (J .EQ. N) THEN
           TMP1 = REAL(A(J,J-1))
           TMP2 = AIMAG(A(J,J-1))
           CALL SLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP6)
           C = SIGN(C,TMP1)
           S = SIGN(S,TMP2)
           TMP9 = CMPLX(C,S)
           TMP11 = CONJG(TMP9)
           CALL CSCAL(N+1-J,TMP11,A(J,J),N)
        ELSE
           TMP6 = REAL(A(J,J-1))
        ENDIF
        
        TMP1 = REAL(A(J-1,J-1))
        TMP2 = AIMAG(A(J-1,J-1))
        CALL SLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP3)
        C = SIGN(C,TMP1)
        S = SIGN(S,TMP2)
        TMP14(J) = CMPLX(C,S)
        TMP11 = CONJG(TMP14(J))

        CALL SLARTG(TMP3,TMP6,C2(J),S2(J),TMP3)
        A(J-1,J-1) = TMP3
        A(J,J-1) = CZERO

        IF (C2(J) .GE. ABS(S2(J))) THEN
           F=S2(J)/(SONE+C2(J))
           CALL CSROTU(N-J+1,A(J-1,J),N,A(J,J),N,F,S2(J),TMP11)
        ELSE IF (S2(J) .GE. SZERO) THEN
           F=C2(J)/(SONE+S2(J))
           CALL CSROTV(N-J+1,A(J-1,J),N,A(J,J),N,C2(J),F,TMP11)
        ELSE IF (S2(J) .LE. SZERO) THEN
           F=C2(J)/(SONE-S2(J))
           CALL CSROTW(N-J+1,A(J-1,J),N,A(J,J),N,C2(J),F,TMP11)
        ENDIF

     ENDDO

     DO J=N,I+2,-1

        IF (J .EQ. N) THEN
            CALL CSCAL(N,TMP8,V(1,J),1)
        ENDIF

        IF (C1(J) .GE. ABS(S1(J))) THEN
           F=S1(J)/(SONE+C1(J))
           CALL CSROTU(N,V(1,J-1),1,V(1,J),1,F,S1(J),TMP7(J))
        ELSE IF (S1(J) .GE. SZERO) THEN
           F=C1(J)/(SONE+S1(J))
           CALL CSROTV(N,V(1,J-1),1,V(1,J),1,C1(J),F,TMP7(J))
        ELSE IF (S1(J) .LE. SZERO) THEN
           F=C1(J)/(SONE-S1(J))
           CALL CSROTW(N,V(1,J-1),1,V(1,J),1,C1(J),F,TMP7(J))
        ENDIF

     ENDDO
     
     DO J=N,I+2,-1
        
        IF (J .EQ. N) THEN
           CALL CSCAL(N,TMP9,U(1,J),1)
        ENDIF
        
        IF (C2(J) .GE. ABS(S2(J))) THEN
           F=S2(J)/(SONE+C2(J))
           CALL CSROTU(N,U(1,J-1),1,U(1,J),1,F,S2(J),TMP14(J))
        ELSE IF (S2(J) .GE. SZERO) THEN
           F=C2(J)/(SONE+S2(J))
           CALL CSROTV(N,U(1,J-1),1,U(1,J),1,C2(J),F,TMP14(J))
        ELSE IF (S2(J) .LE. SZERO) THEN
           F=C2(J)/(SONE-S2(J))
           CALL CSROTW(N,U(1,J-1),1,U(1,J),1,C2(J),F,TMP14(J))
        ENDIF

     ENDDO

  ENDDO

  TMP1 = REAL(A(1,1))
  TMP2 = AIMAG(A(1,1))
  CALL SLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP3)
  C=SIGN(C,TMP1)
  S=SIGN(S,TMP2)
  A(1,1) = TMP3
  TMP8 = CMPLX(C,-S)
  CALL CSCAL(N,TMP8,V(1,1),1)

  TMP1 = REAL(A(N-1,N))
  TMP2 = AIMAG(A(N-1,N))
  CALL SLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP3)
  C=SIGN(C,TMP1)
  S=SIGN(S,TMP2)
  A(N-1,N) = TMP3
  TMP8 = CMPLX(C,-S)
  A(N,N) = A(N,N)*TMP8
  CALL CSCAL(N,TMP8,V(1,N),1)

  TMP1 = REAL(A(N,N))
  TMP2 = AIMAG(A(N,N))
  CALL SLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP3)
  C=SIGN(C,TMP1)
  S=SIGN(S,TMP2)
  A(N,N) = TMP3
  TMP8 = CMPLX(C,S)
  CALL CSCAL(N,TMP8,U(1,N),1)

  RETURN
END SUBROUTINE CBIDIAG

PROGRAM TESTSJS
  
  use omp_lib
  IMPLICIT NONE
  
  INTEGER I, J, K, N, INFO, SEED
  REAL TMP1, TMP2
  COMPLEX E,S
  COMPLEX CZERO,CONE
  PARAMETER (CZERO = (0.0E0,0.0E0),CONE = (1.0E0,0.0E0))
  REAL SZERO, SONE, DTWO
  PARAMETER (SZERO = 0.0E0, SONE=1.0E0, DTWO=2.0E0)
  COMPLEX,ALLOCATABLE :: A(:,:), B(:,:), U(:,:), V(:,:)
  REAL,ALLOCATABLE :: C(:,:), F(:,:), G(:,:), SIGMA(:)
  EXTERNAL USERRAND
  DOUBLE PRECISION USERRAND

  EXTERNAL CDOTC
  COMPLEX CDOTC
  
  character*1000 ARG1,ARG2

  integer t1, t2, t_rate, t_max, diff
  
  CALL GETARG(1,ARG1)
  CALL GETARG(2,ARG2)

  read (ARG1,*) N
  read (ARG2,*) SEED
  
  ALLOCATE(A(N,N), B(N,N), U(N,N), V(N,N), C(N,N), SIGMA(N))
  ALLOCATE(F(N,N), G(N,N))

  DO I=1,SEED
     CALL XOR()
  ENDDO
  
  DO I=1,N
     DO J=1,I
        TMP1 = USERRAND()
        TMP2 = USERRAND()
        A(J,I) = CMPLX(TMP1,TMP2)
        B(J,I) = A(J,I)
     ENDDO
     DO J=I+1,N
        A(J,I) = CZERO
        B(J,I) = CZERO
     ENDDO
  ENDDO

  DO J=1,N
     DO K=1,N
        U(K,J) = SZERO
     ENDDO
     U(J,J) = SONE
  ENDDO
  
  DO J=1,N
     DO K=1,N
        V(K,J) = SZERO
     ENDDO
     V(J,J) = SONE
  ENDDO

  DO J=1,N
     DO K=1,N
        F(K,J) = SZERO
     ENDDO
     F(J,J) = SONE
  ENDDO

  DO J=1,N
     DO K=1,N
        G(K,J) = SZERO
     ENDDO
     G(J,J) = SONE
  ENDDO
  
  call system_clock(t1)     ! 開始時を記録

  CALL CBIDIAG(N,A,U,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

  print "(A, F10.3)", "time it took was:", diff/dble(t_rate)

  call system_clock(t1)     ! 開始時を記録

  C = SZERO
  
  DO I=1,N
     C(I,I) = REAL(A(I,I))
  ENDDO

  DO I=1,N-1
     C(I,I+1) = REAL(A(I,I+1))
  ENDDO

  CALL SJASVD(N,C,SIGMA,F,G,INFO)

  U = MATMUL(U,F)
  V = MATMUL(V,G)
  
  call system_clock(t2, t_rate, t_max) ! 終了時を記録
  if ( t2 < t1 ) then
     diff = (t_max - t1) + t2 + 1
  else
     diff = t2 - t1
  endif
  
  print "(A, F10.3)", "time it took was:", diff/dble(t_rate)

  IF (INFO .EQ. -1) THEN
     WRITE (*,*) 'not converge'
     STOP
  ENDIF
  
  CALL CHERK ('U', 'C', N, N, SONE, U, N, SZERO, A, N)
  DO I=1,N
     A(I,I) = A(I,I) - CONE
  ENDDO

  E = CZERO
  !$OMP PARALLEL DO REDUCTION(+:E)
  DO I=2,N
     E = E+CDOTC(I-1,A(1,I),1,A(1,I),1)
  ENDDO
  !$OMP END PARALLEL DO
  E = DTWO*E
  DO I=1,N
     E = E+CONJG(A(I,I))*A(I,I)
  ENDDO

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

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

  E = CZERO
  !$OMP PARALLEL DO REDUCTION(+:E)
  DO I=2,N
     E = E+CDOTC(I-1,A(1,I),1,A(1,I),1)
  ENDDO
  !$OMP END PARALLEL DO
  E = DTWO*E
  DO I=1,N
     E = E+CONJG(A(I,I))*A(I,I)
  ENDDO

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

  !$OMP PARALLEL DO
  DO I = 1, N
     CALL CSSCAL(N,SIGMA(I),U(1,I),1)
  ENDDO
  !$OMP END PARALLEL DO

  CALL CGEMM('N','C',N,N,N,CONE,U,N,V,N,CZERO,A,N)

  E = CZERO
  !$OMP PARALLEL DO PRIVATE(S,J) REDUCTION(+:E)
  DO I = 1, N
     DO J = 1, N
        S = B(J, I)-A(J, I)
        E = E+CONJG(S)*S
     ENDDO
  ENDDO
  !$OMP END PARALLEL DO
  
  WRITE (*,*) 'ERROR',SQRT(REAL(E))
  
END PROGRAM TESTSJS