PROGRAM ZUNTH
      IMPLICIT NONE
      INTEGER I,J,K,L,M,N,Z,INFO,LWORK,matrixseed
      PARAMETER (LWORK=1000000000)
      DOUBLE PRECISION DONE,DZERO,DTWO,EPS,DTMP,DTMP2
      PARAMETER (DONE=1.0D0,DZERO=0.0D0,DTWO=2.0D0)
      COMPLEX*16 ZONE,ZZERO,TMP,TMP2
      PARAMETER (ZONE=(1.0D0,0.0D0),ZZERO=(0.0D0,0.0D0))
      COMPLEX*16,ALLOCATABLE::A(:,:),F(:,:),B(:,:),U(:,:),U2(:,:),
     $     VT(:,:),TAU(:),G(:),H(:),C(:,:),WORK(:)
      DOUBLE PRECISION,ALLOCATABLE::WORK2(:),D(:),E(:),Y(:),S(:),T(:)
      DOUBLE PRECISION DLAMCH
      integer          irand
      
      integer seedsize
      integer,allocatable:: seed(:)

      EXTERNAL DLAMCH,ZDOTC
      COMPLEX*16 ZDOTC

      integer t1, t2, t_rate, t_max, diff

      CHARACTER*1000   ARG1, ARG2

      EPS = DLAMCH( 'Precision' )

      CALL GETARG(1,ARG1)
      CALL GETARG(2,ARG2)

      ALLOCATE (WORK(LWORK),WORK2(LWORK))

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

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

         read (ARG1,*) 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 ), Y(N) )

c         DO I=1,N
c            DO J=1,M
c               CALL RANDOM_NUMBER(DTMP)
c               CALL RANDOM_NUMBER(DTMP2)
c               A(J,I)=dcmplx(DTMP,DTMP2)
c            ENDDO
c         ENDDO
c         GO TO 999

         ALLOCATE (U(M,N),D(N),E(N),VT(N,N),TAU(N))

         DO I=1,N
            DO J=1,M
               CALL RANDOM_NUMBER(DTMP)
               CALL RANDOM_NUMBER(DTMP2)
               U(J,I)=dcmplx(DTMP,DTMP2)
            ENDDO
         ENDDO

         DO I=1,N
            DO J=1,N
               CALL RANDOM_NUMBER(DTMP)
               CALL RANDOM_NUMBER(DTMP2)
               VT(J,I)=dcmplx(DTMP,DTMP2)
            ENDDO
         ENDDO

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

         DO I=N-10,N-1
c              E(I+1)=ONE/DBLE(I+1)**101
            Y(I+1)=(EPS**(DBLE(I)/DBLE(N-1)))**19
            D(I+1)=Y(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
            DTMP=D(J)
            D(J)=D(K)
            D(K)=DTMP
         ENDDO

         INFO=0

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

         INFO=0

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

         INFO=0

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

         INFO=0

         CALL ZUNGQR(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 ZGEMM('N','C',M,N,N,ZONE,U,M,VT,N,ZZERO,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
       
 999  close(18)

      ALLOCATE (B(N,N),U(M,N),U2(M,N),D(N),E(N),VT(N,N),
     $     G(N),H(N),C(N,N),F(M,N),S(N),T(N))

      F = A

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

      CALL ZSCQR( M, N, A, M, B, N, WORK, WORK2, INFO )

      CALL ZGEBRD( N, N, B, N, D, E, G, H, WORK, LWORK, INFO)

      CALL ZUNGBR( 'Q', N, N, N, B, N, G, WORK, LWORK, INFO)

      CALL DCOPY( N, D, 1, S, 1)
      CALL DCOPY( N, E, 1, T, 1)

      CALL DDQDS( N, S, T, WORK, INFO)
      
      DTMP = S(2)/S(1)
      L = 1
      DO I = 2, N-1
         DTMP2 = S(I+1)/S(I)
         IF (DTMP2 .LE. DTMP) THEN
            DTMP = DTMP2
            L = I
         ENDIF
      ENDDO
      L = N - L
      
      CALL ZOQDSQ( L, N, D, E, B, N, WORK, INFO)

      WRITE (*,*) 'Singular Values'
      
      DO I = N-L+1, N
         WRITE (*,*) D(I),S(I)
      ENDDO

      CALL ZGEMM( 'N', 'N', M, N, N, ZONE, A, M, B, N, ZZERO, 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:N,1:N) = ZZERO
      DO I = 1, N
         B(I,I) = -ZONE
      ENDDO
      CALL ZHERK('U','C',N,M,DONE,U,M,DONE,B,N)

      TMP = ZZERO
!$OMP PARALLEL DO PRIVATE( J ) REDUCTION( +:TMP )      
      DO I = 2, N
         TMP = TMP+ZDOTC(I-1,B(1,I),1,B(1,I),1)
      ENDDO
!$OMP END PARALLEL DO
      TMP = DTWO*TMP
      DO I = 1, N
         TMP = TMP+CONJG(B(I,I))*B(I,I)
      ENDDO
      WRITE (*,*) 'ERROR',SQRT(DBLE(DTMP2))
      
      print "(A, F10.3)", "time it took was:", diff/dble(t_rate)
      
      END PROGRAM ZUNTH