PROGRAM DDPSVD
      IMPLICIT NONE
      INTEGER I,J,K,M,N,INFO,LWORK,W,W2,BL,Z,matrixseed
      INTEGER IL,IU
      PARAMETER (LWORK=2000000000)
      DOUBLE PRECISION ZERO,ONE,TWO,TMP,TMP2,EPS
      DOUBLE PRECISION VL,VU,ABSTOL
      PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0)
      DOUBLE PRECISION,ALLOCATABLE::A(:,:),F(:,:),B(:,:),U(:,:),U2(:,:),
     $     D(:),E(:),VT(:,:),TAU(:),C(:,:),G(:),H(:),D2(:),E2(:),TV(:,:)
      DOUBLE PRECISION,ALLOCATABLE::WORK(:)
      INTEGER,ALLOCATABLE::IWORK(:),IFAIL(:)
      DOUBLE PRECISION DLAMCH,DNRM2,DDOT,DFMA0
      integer            irand

      EXTERNAL DLAMCH,DNRM2,DDOT,DFMA0
      
      integer seedsize
      integer,allocatable:: seed(:)

      integer t1, t2, t_rate, t_max, diff

      CHARACTER*1000   ARG1, ARG2, ARG3, ARG4

      EPS = DLAMCH( 'Precision' )

      CALL GETARG(1,ARG1)
      CALL GETARG(2,ARG2)
      CALL GETARG(3,ARG3)
      CALL GETARG(4,ARG4)
      READ (ARG1,*) W
      READ (ARG2,*) BL

      ALLOCATE (WORK(LWORK))

      open(18,file=ARG4,status='old')

      READ (18,*) m
      READ (18,*) n
      READ (18,*) z
      
      if (.TRUE.) then

         read (ARG3,*) matrixseed
         call random_seed(size=seedsize)
         allocate(seed(seedsize))
         DO i=1,seedsize
            seed(i)=matrixseed
         end do
         call random_seed(put=seed)

         ALLOCATE( A( m, n ) )
         ALLOCATE (U(M,N),D(N),E(N),VT(N,N),TAU(N))

         CALL RANDOM_NUMBER(U)
         CALL RANDOM_NUMBER(VT)

         DO I=0,N-1
c              E(I+1)=ONE/DBLE(I+1)**101
            E(I+1)=(EPS**(DBLE(I)/DBLE(N-1)))**19
c            E(I+1)=(EPS**(DBLE(I)/DBLE(N-1)))
            D(I+1)=E(I+1)
c              E(I+1)=(ONE-EPS)*(DBLE(N-1-I)/DBLE(N-1))+EPS
         ENDDO

         WRITE (*,*) 'CONDITION',D(1),D(N),D(1)/D(N)

         DO I=1,N*N
            J=MOD(irand(0),N)+1
            K=MOD(irand(0),N)+1
            TMP=D(J)
            D(J)=D(K)
            D(K)=TMP
         ENDDO

         INFO=0

         CALL DGEQRF(M,N,U,M,TAU,WORK,LWORK,INFO)

         INFO=0

         CALL DORGQR(M,N,N,U,M,TAU,WORK,LWORK,INFO)

         INFO=0

         CALL DGEQRF(N,N,VT,N,TAU,WORK,LWORK,INFO)

         INFO=0

         CALL DORGQR(N,N,N,VT,N,TAU,WORK,LWORK,INFO)

         DO I=1,N
            DO J=1,N
               VT(J,I)=VT(J,I)*D(I)
            ENDDO
         ENDDO

         INFO=0

         CALL DGEMM('N','T',M,N,N,ONE,U,M,VT,N,ZERO,A,M)

         DEALLOCATE (U,D,E,VT,TAU)

      ELSE

         ALLOCATE( A( m, n ) )
         
 1000    READ (18,*, end=1001) i,j,TMP2
         
         A(i+1,j+1)=TMP2
         
         go to 1000
         
 1001    continue
         
      endif
       
      close(18)

      ALLOCATE (B(N,N),U(M,W),U2(M,W),D(N),E(N),VT(W,N),TAU(W),
     $     C(N,N),G(N),H(N),D2(2*N),E2(2*N),TV(2*N,W),
     $     IWORK(MAX(2*BL*N+2*BL+3*2*N,5*2*N)),IFAIL(2*N),F(M,N))

      ABSTOL = TWO*DLAMCH( 'S' )

      F = A

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

      CALL DSCQR( M, N, A, M, C, N, WORK, IWORK, INFO )

      CALL DGEBRD( N, N, C, N, D, E, G, H, WORK, LWORK, INFO)

      DO I=1,2*N
         D2(I)=ZERO
      ENDDO

      DO I=1,N
         E2(2*I-1)=D(I)
      ENDDO

      DO I=1,N-1
         E2(2*I)=E(I)
      ENDDO

      IL = 1
      IU = W
      
      CALL DSTEVX2('V','I',2*N,D2,E2,VL,VU,IL,IU,ABSTOL,
     $     W2,TAU,TV,2*N,BL,WORK,IWORK,IFAIL,INFO)

      TAU(1:W) = -TAU(1:W)

!$OMP PARALLEL DO PRIVATE( TMP )      
      DO I=1,W
         CALL DCOPY(N,TV(2,I),2,B(1,I),1)
         TMP = DNRM2(N,B(1,I),1)
         CALL DSCAL(N,ONE/TMP,B(1,I),1)
         CALL DCOPY(N,TV(1,I),2,VT(I,1),W)
         TMP = -DNRM2(N,VT(I,1),W)
         CALL DSCAL(N,ONE/TMP,VT(I,1),W)
      ENDDO
!$OMP END PARALLEL DO

      CALL DORMBR( 'Q', 'L', 'N', N, W, N, C, N,
     $     G, B, N, WORK, LWORK, INFO )
      CALL DORMBR( 'P', 'R', 'T', W, N, N, C, N,
     $     H, VT, W, WORK, LWORK, INFO )
      
      CALL DGEMM('N','N',M,W,N,ONE,A,M,B,N,ZERO,U,M)
      
      call system_clock(t2, t_rate, t_max) ! 終了時を記録
      if ( t2 < t1 ) then
         diff = (t_max - t1) + t2 + 1
      else
         diff = t2 - t1
      endif

      B(1:W,1:W) = ZERO
      DO I = 1, W
         B(I,I) = -ONE
      ENDDO
      CALL DSYRK('U','T',W,M,ONE,U,M,ONE,B,N)
      TMP = ZERO
!$OMP PARALLEL DO REDUCTION( +:TMP )      
      DO I = 2, W
         TMP = TMP+DDOT(I-1,B(1,I),1,B(1,I),1)
      ENDDO
!$OMP END PARALLEL DO
      TMP = TWO*TMP
      DO I = 1, W
         TMP = DFMA0(B(I,I),B(I,I),TMP)
      ENDDO
      WRITE (*,*) 'ERROR',SQRT(TMP)
      
      B(1:W,1:W) = ZERO
      DO I = 1, W
         B(I,I) = -ONE
      ENDDO
      CALL DSYRK('U','N',W,N,ONE,VT,W,ONE,B,N)
      TMP = ZERO
!$OMP PARALLEL DO REDUCTION( +:TMP )      
      DO I = 2, W
         TMP = TMP+DDOT(I-1,B(1,I),1,B(1,I),1)
      ENDDO
!$OMP END PARALLEL DO
      TMP = TWO*TMP
      DO I = 1, W
         TMP = DFMA0(B(I,I),B(I,I),TMP)
      ENDDO
      WRITE (*,*) 'ERROR',SQRT(TMP)

      CALL DGEMM('N','T',M,W,N,ONE,F,M,VT,W,ZERO,U2,M)
      CALL DGEMM('T','N',W,N,M,ONE,U,M,F,M,ZERO,B,N)
      TMP = ZERO
!$OMP PARALLEL DO REDUCTION( +:TMP )      
      DO I = 1, W
         CALL DAXPY(M,-TAU(I),U(1,I),1,U2(1,I),1)
         TMP = TMP+DDOT(M,U2(1,I),1,U2(1,I),1)
         CALL DAXPY(N,-TAU(I),VT(I,1),W,B(I,1),N)
         TMP = TMP+DDOT(N,B(I,1),N,B(I,1),N)
      ENDDO
!$OMP END PARALLEL DO

      WRITE (*,*) 'ERROR',SQRT(TMP)
      
      print "(A, F10.3)", "time it took was:", diff/dble(t_rate)
      
      END PROGRAM DDPSVD