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

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

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

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

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

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

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

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

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

SUBROUTINE ZBIDIAG(N,A,U,V,INFO)
  
  use omp_lib
  IMPLICIT NONE
  
  INTEGER N, INFO
  COMPLEX*16 A(N,N),V(N,N),U(N,N)
  
  DOUBLE PRECISION DONE
  PARAMETER (DONE = 1.0D0)

  COMPLEX*16 ZZERO
  PARAMETER (ZZERO = (0.0D0,0.0D0))
  INTEGER I, J

  DOUBLE PRECISION TMP1, TMP2, TMP3, TMP6, TMP13
  COMPLEX*16 TMP7(N),TMP8,TMP9,TMP11,TMP14(N)
  
  DOUBLE PRECISION DZERO
  PARAMETER (DZERO = 0.0D0)
  
  DOUBLE PRECISION 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 = DBLE(A(I,J))
           TMP2 = DIMAG(A(I,J))
           CALL DLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP13)
           C = DSIGN(C,TMP1)
           S = DSIGN(S,TMP2)
           TMP8 = DCMPLX(C,-S)
           CALL ZSCAL(J-I,TMP8,A(I+1,J),1)
        ENDIF
           
        TMP1 = DBLE(A(I,J-1))
        TMP2 = DIMAG(A(I,J-1))
        CALL DLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP3)
        C = DSIGN(C,TMP1)
        S = DSIGN(S,TMP2)
        TMP7(J) = DCMPLX(C,-S)

        CALL DLARTG(TMP3,TMP13,C1(J),S1(J),TMP13)
        A(I,J-1) = TMP13
        A(I,J) = ZZERO

        IF (C1(J) .GE. ABS(S1(J))) THEN
           F=S1(J)/(DONE+C1(J))
           CALL ZDROTU(J-I,A(I+1,J-1),1,A(I+1,J),1,F,S1(J),TMP7(J))
        ELSE IF (S1(J) .GE. DZERO) THEN
           F=C1(J)/(DONE+S1(J))
           CALL ZDROTV(J-I,A(I+1,J-1),1,A(I+1,J),1,C1(J),F,TMP7(J))
        ELSE IF (S1(J) .LE. DZERO) THEN
           F=C1(J)/(DONE-S1(J))
           CALL ZDROTW(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 = DBLE(A(J,J-1))
           TMP2 = DIMAG(A(J,J-1))
           CALL DLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP6)
           C = DSIGN(C,TMP1)
           S = DSIGN(S,TMP2)
           TMP9 = DCMPLX(C,S)
           TMP11 = DCONJG(TMP9)
           CALL ZSCAL(N+1-J,TMP11,A(J,J),N)
        ELSE
           TMP6 = DBLE(A(J,J-1))
        ENDIF
        
        TMP1 = DBLE(A(J-1,J-1))
        TMP2 = DIMAG(A(J-1,J-1))
        CALL DLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP3)
        C = DSIGN(C,TMP1)
        S = DSIGN(S,TMP2)
        TMP14(J) = DCMPLX(C,S)
        TMP11 = DCONJG(TMP14(J))

        CALL DLARTG(TMP3,TMP6,C2(J),S2(J),TMP3)
        A(J-1,J-1) = TMP3
        A(J,J-1) = ZZERO

        IF (C2(J) .GE. ABS(S2(J))) THEN
           F=S2(J)/(DONE+C2(J))
           CALL ZDROTU(N-J+1,A(J-1,J),N,A(J,J),N,F,S2(J),TMP11)
        ELSE IF (S2(J) .GE. DZERO) THEN
           F=C2(J)/(DONE+S2(J))
           CALL ZDROTV(N-J+1,A(J-1,J),N,A(J,J),N,C2(J),F,TMP11)
        ELSE IF (S2(J) .LE. DZERO) THEN
           F=C2(J)/(DONE-S2(J))
           CALL ZDROTW(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 ZSCAL(N,TMP8,V(1,J),1)
        ENDIF
           
        IF (C1(J) .GE. ABS(S1(J))) THEN
           F=S1(J)/(DONE+C1(J))
           CALL ZDROTU(N,V(1,J-1),1,V(1,J),1,F,S1(J),TMP7(J))
        ELSE IF (S1(J) .GE. DZERO) THEN
           F=C1(J)/(DONE+S1(J))
           CALL ZDROTV(N,V(1,J-1),1,V(1,J),1,C1(J),F,TMP7(J))
        ELSE IF (S1(J) .LE. DZERO) THEN
           F=C1(J)/(DONE-S1(J))
           CALL ZDROTW(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 ZSCAL(N,TMP9,U(1,J),1)
        ENDIF
        
        IF (C2(J) .GE. ABS(S2(J))) THEN
           F=S2(J)/(DONE+C2(J))
           CALL ZDROTU(N,U(1,J-1),1,U(1,J),1,F,S2(J),TMP14(J))
        ELSE IF (S2(J) .GE. DZERO) THEN
           F=C2(J)/(DONE+S2(J))
           CALL ZDROTV(N,U(1,J-1),1,U(1,J),1,C2(J),F,TMP14(J))
        ELSE IF (S2(J) .LE. DZERO) THEN
           F=C2(J)/(DONE-S2(J))
           CALL ZDROTW(N,U(1,J-1),1,U(1,J),1,C2(J),F,TMP14(J))
        ENDIF
        
     ENDDO

  ENDDO

  TMP1 = DBLE(A(1,1))
  TMP2 = DIMAG(A(1,1))
  CALL DLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP3)
  C=SIGN(C,TMP1)
  S=SIGN(S,TMP2)
  A(1,1) = TMP3
  TMP8 = DCMPLX(C,-S)
  CALL ZSCAL(N,TMP8,V(1,1),1)

  TMP1 = DBLE(A(N-1,N))
  TMP2 = DIMAG(A(N-1,N))
  CALL DLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP3)
  C=SIGN(C,TMP1)
  S=SIGN(S,TMP2)
  A(N-1,N) = TMP3
  TMP8 = DCMPLX(C,-S)
  A(N,N) = A(N,N)*TMP8
  CALL ZSCAL(N,TMP8,V(1,N),1)

  TMP1 = DBLE(A(N,N))
  TMP2 = DIMAG(A(N,N))
  CALL DLARTG(ABS(TMP1),ABS(TMP2),C,S,TMP3)
  C=SIGN(C,TMP1)
  S=SIGN(S,TMP2)
  A(N,N) = TMP3
  TMP8 = DCMPLX(C,S)
  CALL ZSCAL(N,TMP8,U(1,N),1)
  
  RETURN
END SUBROUTINE ZBIDIAG

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

  EXTERNAL ZDOTC
  COMPLEX*16 ZDOTC
  
  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))

  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,N
        A(J,I) = ZZERO
        B(J,I) = ZZERO
     ENDDO
  ENDDO

  DO J=1,N
     DO K=1,N
        U(K,J) = ZZERO
     ENDDO
     U(J,J) = ZONE
  ENDDO
  
  DO J=1,N
     DO K=1,N
        V(K,J) = ZZERO
     ENDDO
     V(J,J) = ZONE
  ENDDO
  
  call system_clock(t1)     ! 開始時を記録
  
  CALL ZBIDIAG(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)

  ALLOCATE(C(N,N), SIGMA(N))
  ALLOCATE(F(N,N), G(N,N))

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

  C = DZERO
  
  DO I=1,N
     C(I,I) = DBLE(A(I,I))
  ENDDO

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

  CALL DJASVD(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 ZHERK ('U', 'C', N, N, DONE, U, N, DZERO, A, N)
  DO I=1,N
     A(I,I) = A(I,I) - ZONE
  ENDDO

  E = ZZERO
  !$OMP PARALLEL DO REDUCTION(+:E)
  DO I=2,N
     E = E+ZDOTC(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',DSQRT(DBLE(E))

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

  E = ZZERO
  !$OMP PARALLEL DO REDUCTION(+:E)
  DO I=2,N
     E = E+ZDOTC(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',DSQRT(DBLE(E))

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

  CALL ZGEMM('N','C',N,N,N,ZONE,U,N,V,N,ZZERO,A,N)

  E = ZZERO
  !$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',DSQRT(DBLE(E))
  
END PROGRAM TESTSJS