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

  IMPLICIT NONE

  REAL ZERO, ONE, HALF
  PARAMETER (ONE = 1.0E0, ZERO = 0.0E0, HALF = 0.5E0)
  REAL F1, G1, C, S, R1, U, V
  EXTERNAL SFMA0
  REAL SFMA0

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

  RETURN

END SUBROUTINE SLARTGU

SUBROUTINE SJASVD(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 A(M,N),V(N,N),SIGMA(N),WORK(M),CS(N),SN(N)
  
  REAL ZERO, ONE, HALF
  PARAMETER (ONE = 1.0E0, ZERO = 0.0E0, HALF = 0.5E0)
  INTEGER I, J, K, L, P, Q, R, ITER, MAXITER
  PARAMETER (MAXITER = 100)
  
  REAL T
  REAL TMP1, TMP2, TMP3, TMP10, TMP11, TMP12, TMP13, TMP14, TMP15
  
  REAL H1, G1, R1, X, Y, Z, S1, S2
  REAL C, S
  REAL EPS, SAFMIN, SCALE1, SIGMX
  
  EXTERNAL SFMA0, SLAMCH, SDOT, SNRM2
  REAL SFMA0, SLAMCH, SDOT, SNRM2

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

  SCALE1 = ONE/REAL(M)/SQRT(SAFMIN)
  SIGMX = MAXVAL(ABS(A(1:M,1:N)))
  IF (SIGMX .GT. ZERO) THEN
     CALL SLASCL( 'G', 0, 0, SIGMX, SCALE1, M, N, &
          A, M, IINFO )
  ENDIF
  
  DO I = 1,N
     SIGMA(I) = SNRM2(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. ZERO) 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 SSWAP(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) .EQ. ZERO) THEN
           Q = 0
           GO TO 40
        ENDIF

        IF (SIGMA(J) .LT. SAFMIN) THEN
           TMP1 = ONE/SAFMIN
        ELSE
           TMP1 = ONE/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. ZERO) THEN
              SN(P)=ZERO
              CYCLE
           ENDIF
           TMP2 = SDOT(M,WORK(1),1,A(1,P),1)
           GO TO 20
        ENDDO
        Q = N
        GO TO 40

20      K = P

        IF (ABS(TMP2) .LT. SAFMIN) THEN
           IF (ABS(TMP2) .LE. SQRT(REAL(M))*EPS*SIGMA(K)*S1) THEN
              SN(K)=ZERO
              GO TO 10
           ENDIF
        ELSE
           IF (ABS(TMP2) .LE. EPS*SIGMA(K)*S1) THEN
              SN(K)=ZERO
              GO TO 10
           ENDIF
        ENDIF

        S2 = SIGMA(K)*TMP1
        H1 = (SIGMA(J) - SIGMA(K))*((S1+S2)*HALF)
        G1 = TMP2
        
        CALL SLARTGU (ABS(H1), ABS(G1), C, S, R1)
        T = G1 / (H1 + SIGN(R1,H1))
        R1 = SQRT(SFMA0(T,T,ONE))
        C = ONE/R1
        S = T/R1
        
        TMP10 = SIGMA(J)*SQRT(SFMA0(T,TMP2/(SIGMA(J)*S1),ONE))
        IF (TMP10 .LT. SAFMIN) THEN
           TMP10 = SAFMIN
        ENDIF
        TMP1 = ONE/TMP10
        
        TMP15 = SFMA0(-T,TMP2/(SIGMA(K)*S2),ONE)
        IF (TMP15 .GT. ZERO) THEN
           TMP11 = SIGMA(K)*SQRT(TMP15)
           IF (TMP11 .LT. SAFMIN) THEN
              TMP11 = SAFMIN
           ENDIF
           TMP12 = ONE/TMP11
        ELSE
           TMP11 = ZERO
        ENDIF
        
        CS(K)=C
        SN(K)=S
        
        IF (S .NE. ZERO) THEN
           
           X = S/(ONE+C)
        
           IPIV2(J) = 1
           IPIV2(K) = 1

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

              IF (TMP11 .GT. ZERO) THEN
                 TMP13 = ZERO
                 TMP2 = ZERO
                 TMP14 = ZERO
                 DO L = 1,M
                    Y = A(L,J)
                    Z = A(L,K)
                    A(L,J) = S*(-X*Y+Z)+Y
                    WORK(L) = A(L,J)*TMP1
                    TMP13 = TMP13+WORK(L)**2
                    TMP2 = TMP2+WORK(L)*A(L,P)
                    A(L,K) = -S*(X*Z+Y)+Z
                    TMP14 = TMP14+(A(L,K)*TMP12)**2
                 ENDDO
                 SIGMA(J) = TMP10*SQRT(TMP13)
                 S1 = SIGMA(J)*TMP1
                 SIGMA(K) = TMP11*SQRT(TMP14)
              ELSE
                 TMP13 = ZERO
                 TMP2 = ZERO
                 DO L = 1,M
                    Y = A(L,J)
                    Z = A(L,K)
                    A(L,J) = S*(-X*Y+Z)+Y
                    WORK(L) = A(L,J)*TMP1
                    TMP13 = TMP13+WORK(L)**2
                    TMP2 = TMP2+WORK(L)*A(L,P)
                    A(L,K) = -S*(X*Z+Y)+Z
                 ENDDO
                 SIGMA(J) = TMP10*SQRT(TMP13)
                 S1 = SIGMA(J)*TMP1
                 SIGMA(K) = SNRM2(M,A(1,K),1)
              ENDIF

              GO TO 30
           ENDDO
           
           IF (TMP11 .GT. ZERO) THEN
              TMP13 = ZERO
              TMP14 = ZERO
              DO L = 1,M
                 Y = A(L,J)
                 Z = A(L,K)
                 A(L,J) = S*(-X*Y+Z)+Y
                 WORK(L) = A(L,J)*TMP1
                 TMP13 = TMP13+WORK(L)**2
                 A(L,K) = -S*(X*Z+Y)+Z
                 TMP14 = TMP14+(A(L,K)*TMP12)**2
              ENDDO
              SIGMA(J) = TMP10*SQRT(TMP13)
              S1 = SIGMA(J)*TMP1
              SIGMA(K) = TMP11*SQRT(TMP14)
           ELSE
              TMP13 = ZERO
              DO L = 1,M
                 Y = A(L,J)
                 Z = A(L,K)
                 A(L,J) = S*(-X*Y+Z)+Y
                 WORK(L) = A(L,J)*TMP1
                 TMP13 = TMP13+WORK(L)**2
                 A(L,K) = -S*(X*Z+Y)+Z
              ENDDO
              SIGMA(J) = TMP10*SQRT(TMP13)
              S1 = SIGMA(J)*TMP1
              SIGMA(K) = SNRM2(M,A(1,K),1)
           ENDIF
                 
30         IF (FLAG .EQ. 0) THEN
              TMP3 = SDOT(M,WORK(1),1,A(1,K),1)
              IF (ABS(TMP3) .LT. SAFMIN) THEN
                 IF (ABS(TMP3) .LE. SQRT(REAL(M))*EPS*SIGMA(K)*S1) FLAG = 1
              ELSE
                 IF (ABS(TMP3) .LE. EPS*SIGMA(K)*S1) FLAG = 1
              ENDIF
           ENDIF
                 
           IF (SIGMA(J) .GT. ZERO .AND. P .LE. N) GO TO 20

           Q = K
           
        ELSE

           GO TO 10

        ENDIF

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

           IF (S .NE. ZERO) THEN
              C=CS(P)
              X = S/(ONE+C)
              DO L = 1,N
                 Y = V(L,J)
                 Z = V(L,P)
                 V(L,J) = S*(-X*Y+Z)+Y
                 V(L,P) = -S*(X*Z+Y)+Z
              ENDDO
           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. ZERO) 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 SSWAP(M, A(1,K), 1, A(1,I), 1)
     CALL SSWAP(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. ZERO) THEN
     CALL SLASCL( 'G', 0, 0, SCALE1, SIGMX, N, 1, SIGMA, N, IINFO )
  ENDIF
  
  INFO = 0
  RETURN
  
END SUBROUTINE SJASVD

PROGRAM TESTSJS
  
  use omp_lib
  IMPLICIT NONE
  
  INTEGER I, J, K, M, N, INFO, SEED, RANK
  REAL E, TMP
  REAL ZERO, ONE, TWO, SAFMIN
  PARAMETER (ZERO = 0.0E0, ONE=1.0E0, TWO=2.0E0)
  REAL,ALLOCATABLE :: SIGMA(:)
  REAL,ALLOCATABLE :: A(:,:), B(:,:), WORK(:,:), V(:,:)
  EXTERNAL SFMA0, SDOT
  REAL SFMA0, SDOT, SLAMCH
  EXTERNAL USERRAND, SLAMCH
  DOUBLE PRECISION USERRAND

  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,M
        TMP = REAL(USERRAND())
        A(J,I) = TMP
        B(J,I) = TMP
     ENDDO
  ENDDO

  DO J=1,N
     DO K=1,N
        V(K,J) = ZERO
     ENDDO
     V(J,J) = ONE
  ENDDO

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

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

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

  CALL SSYRK ('U', 'T', RANK, M, ONE, A, M, ZERO, WORK, N)
  DO I=1,RANK
     WORK(I,I) = WORK(I,I) - ONE
  ENDDO

  E = ZERO
  !$OMP PARALLEL DO REDUCTION(+:E)
  DO I=2,RANK
     E = E+SDOT(I-1,WORK(1,I),1,WORK(1,I),1)
  ENDDO
  !$OMP END PARALLEL DO
  E = TWO*E
  DO I=1,RANK
     E = SFMA0(WORK(I,I),WORK(I,I),E)
  ENDDO
  
  WRITE (*,*) 'ERROR',SQRT(E)

  CALL SSYRK ('U', 'T', RANK, N, ONE, V, N, ZERO, WORK, N)
  DO I=1,RANK
     WORK(I,I) = WORK(I,I) - ONE
  ENDDO

  E = ZERO
  !$OMP PARALLEL DO REDUCTION(+:E)
  DO I=2,RANK
     E = E+SDOT(I-1,WORK(1,I),1,WORK(1,I),1)
  ENDDO
  !$OMP END PARALLEL DO
  E = TWO*E
  DO I=1,RANK
     E = SFMA0(WORK(I,I),WORK(I,I),E)
  ENDDO

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