module my_data
  implicit none
  INTEGER, ALLOCATABLE:: IAP(:), JA(:)
  Double precision, ALLOCATABLE:: A(:), DA(:,:), ax(:)
  Integer mode
  integer, allocatable:: start_row(:)
  Integer m, n
  Double precision s0,s1
end module my_data

program dsvd_file
  Use my_data
  use omp_lib
  IMPLICIT NONE

  integer          ldv, ldu, iparam(11), ipntr(11)
  Double precision t1
  Double precision, ALLOCATABLE:: &
       v(:,:), u(:,:), workl(:), workd(:), &
       s(:,:), resid(:)
  logical, ALLOCATABLE:: select(:)
  integer, allocatable:: seed(:)
!     
!     %---------------%
!     | Local Scalars |
!     %---------------%
!
  CHARACTER*1000   ARG1, ARG2, ARG3, ARG4, ARG5, ARG6
  character        bmat*1, which*2
  integer          ido, nev, ncv, lworkl, info, ierr, &
       accuracy, matrixseed, seedsize, each, &
       j, ishfts, maxitr, mode1, nconv, i, k, w, y
  integer          maxthreads, n0, t, k1, k0, l0, l1
  logical          rvec
  Double precision tol, sigma, temp, temp2, sum0, sum1, alpha
  Double precision epsilon
!
!     %------------%
!     | Parameters |
!     %------------%
!
  Double precision two, one, zero
  integer          lwork
  parameter        (two = 2.0D+0, one = 1.0D+0, zero = 0.0D+0)
  parameter        (lwork = 2000000000)
!  
!     %-----------------------------%
!     | BLAS & LAPACK routines used |
!     %-----------------------------%
!
  Double precision dnrm2, dlamch
  external         dnrm2, daxpy, dcopy, dscal, dlamch
  LOGICAL          LSAME
  EXTERNAL         LSAME
  integer          irand
!     
!     %-----------------------%
!     | Executable Statements |
!     %-----------------------%
!
!     %-------------------------------------------------%
!     | The following include statement and assignments |
!     | initiate trace output from the internal         |
!     | actions of ARPACK.  See debug.doc in the        |
!     | DOCUMENTS directory for usage.  Initially, the  |
!     | most useful information will be a breakdown of  |
!     | time spent in the various stages of computation |
!     | given by setting msaupd = 1.                    |
!     %-------------------------------------------------%
!
  include 'debug.h'
  ndigit = -3
  logfil = 6
  msgets = 0
  msaitr = 0 
  msapps = 0
  msaupd = 1
  msaup2 = 0
  mseigt = 0
  mseupd = 0
!
!     %-------------------------------------------------%
!     | The following sets dimensions for this problem. |
!     %-------------------------------------------------%
!
  CALL GETARG(1,ARG1)
  CALL GETARG(2,ARG2)
  CALL GETARG(3,ARG3)
  CALL GETARG(4,ARG4)
  CALL GETARG(5,ARG5)
  CALL GETARG(6,ARG6)
  
  READ (ARG1,*) nev
  READ (ARG2,*) ncv 

  WRITE (*,*) nev, ncv

  if (lsame(ARG3,'s')) THEN
     mode = 0
     write (*,*) 'runnning: sparse mode'
  else if(lsame(ARG3,'d')) THEN
     mode = 1
     write (*,*) 'runnning: dense mode'
  else
     write (*,*) 'error: select d[dense] or s[sparse]'
  end if
  
  READ (ARG4,*) accuracy
  
  tol = 10.0D0**(-accuracy)
  
  open(18,file=ARG6,status='old')
  
  READ (18,*) m
  READ (18,*) n
  READ (18,*) w
  
  ldu = m
  ldv = n
  
  WRITE (*,*) m, n, w
  
  ALLOCATE( s(ncv,2) )
  
  ALLOCATE( select(ncv) )
  
  if (.FALSE.) then
     
     read (ARG5,*) matrixseed
     call random_seed(size=seedsize)
     allocate(seed(seedsize))
     DO i=1,seedsize
        seed(i)=matrixseed
     end do
     call random_seed(put=seed)
     
     if(mode==0) then
        
        ALLOCATE( IAP( m+1 ), JA( w ), A( w ) )
        
        each = w / m
        write (*,*) 'each row has',each
        
        DO k=1,m+1
           IAP(k)=1+(k-1)*each
        ENDDO
        
        call random_number(A(1:w))
        JA(1:w)=int(A(1:w)*n)+1
        
        call random_number(A(1:w))
        
     else
        
        ALLOCATE( DA( n, m ) )
        ALLOCATE( U(M,N), V(N,N), AX(N), RESID(N), &
             WORKD(N), workl(lwork) )
        
        CALL RANDOM_NUMBER(U)
        CALL RANDOM_NUMBER(V)
        
        EPSILON = DLAMCH( 'Precision' )
        
        DO I=0,N-1
!              RESID(I+1)=ONE/DBLE(I+1)**101
           RESID(I+1)=(EPSILON**(DBLE(I)/DBLE(N-1)))**19
           AX(I+1)=RESID(I+1)
!              RESID(I+1)=(ONE-EPSILON)*(DBLE(N-1-I)/DBLE(N-1))+EPSILON
        ENDDO

        WRITE (*,*) 'CONDITION',AX(1),AX(N),AX(1)/AX(N)
        
        DO I=1,N*N
           J=MOD(irand(0),N)+1
           K=MOD(irand(0),N)+1
           TEMP=AX(J)
           AX(J)=AX(K)
           AX(K)=TEMP
        ENDDO
        
        INFO=0
        
        CALL DGEQRF(M,N,U,M,WORKD,WORKl,LWORK,INFO)
        
        INFO=0
        
        CALL DORGQR(M,N,N,U,M,WORKD,WORKl,LWORK,INFO)
        
        INFO=0
        
        CALL DGEQRF(N,N,V,N,WORKD,WORKl,LWORK,INFO)
        
        INFO=0
        
        CALL DORGQR(N,N,N,V,N,WORKD,WORKl,LWORK,INFO)
        
        DO I=1,N
           DO J=1,N
              V(J,I)=V(J,I)*AX(I)
           ENDDO
        ENDDO
        
        INFO=0
        
        CALL DGEMM('N','T',N,M,N,ONE,V,N,U,M,ZERO,DA,N)
        
        DEALLOCATE (U,V,AX,RESID,WORKD,workl)
        
     end if
     
  else
     
     if(mode==0) then
        
        ALLOCATE( IAP( m+1 ), JA( w ), A( w ) )
        
        k = 1
        y = 1
        IAP(k) = y
        sum0 = zero
        
998     READ (18,*, end=999) i,j,alpha  
        
        IF (i+1 .NE. k) THEN
           if (k+1 .lt. i+1) write (*,*) 'not row full rank'
           DO k=k+1,i+1
              IAP(k)=y
           enddo
           k = i+1
        END IF
        
        JA(y)=j+1
        sum0=sum0+alpha*alpha
        A(y)=alpha
        y=y+1
        go to 998
        
999     if (k+1 .lt. m+1) write (*,*) 'not row full rank'
        DO k=k+1,m+1
           IAP(k)=y
        enddo
        
     else
        
        ALLOCATE( DA( n, m ) )
        sum0 = zero
        
1000    READ (18,*, end=1001) i,j,alpha
        
        DA(j+1,i+1)=alpha
        sum0=sum0+alpha*alpha
        
        go to 1000
        
1001    continue
        
     end if
     
  end if
  
  close(18)
  
  if(mode==0) then
     
     maxthreads=omp_get_max_threads()
     
     ALLOCATE (start_row(maxthreads+1))
     
     n0 = w/maxthreads
     
     t=0
     i=0
     
     do while (i<maxthreads .and. t<m)
        
        start_row(i+1)=t+1
        
        k1=0
        do while(k1<n0 .and. t<m)
           k0=k1
           k1=k0+(IAP(t+2)-IAP(t+1))
           t=t+1
        end do
        
        if(start_row(i+1) .eq. t) cycle
        
        l0=n0-k0
        
        l1=ABS(k1-n0)
        
        if (l0<l1) t=t-1
        
        i=i+1
        
     end do
     start_row(i+1)=m+1
     
     if(i .ne. maxthreads) then
        i=i+1
        do while(i<maxthreads+1)
           start_row(i+1)=1
           i=i+1
        enddo
     end if
     
     do i=1,maxthreads+1
        write (*,*) i,start_row(i)
     end do
     
  end if
  
  ALLOCATE( workl(ncv*(ncv+8)) )
  
  bmat  = 'I'
  which = 'LA'
  
  lworkl = ncv*(ncv+8)
  info = 0
  ido = 0
  
  ishfts = 1
  maxitr = n
  mode1 = 1
  
  iparam(1) = ishfts
  iparam(3) = maxitr
  iparam(7) = mode1
  
  s0 = zero
  t1 = omp_get_wtime()
  
  ALLOCATE( v(ldv,ncv), u(ldu, nev), workd(3*n),  &
       resid(n), ax(m) )
  
10 continue
  
  call dsaupd ( ido, bmat, n, which, nev, tol, resid,  &
       ncv, v, ldv, iparam, ipntr, workd, workl, &
       lworkl, info )
  
  if (ido .eq. -1 .or. ido .eq. 1) then
     
     call av_atv (workd(ipntr(1)), workd(ipntr(2)))
     
     go to 10
     
  end if
  
  if ( info .lt. 0 ) then
     
     print *, ' '
     print *, ' Error with _saupd, info = ', info
     print *, ' Check documentation in _saupd '
     print *, ' '
     
  else 
     
     rvec = .true.
     
     call dseupd ( rvec, 'All', select, s, v, ldv, sigma,  &
          bmat, n, which, nev, tol, resid, ncv, v, ldv, &
          iparam, ipntr, workd, workl, lworkl, ierr )
     
     if ( ierr .ne. 0) then
        
        print *, ' '
        print *, ' Error with _seupd, info = ', ierr
        print *, ' Check the documentation of _seupd. '
        print *, ' '
        
     else
        
        nconv =  iparam(5)
        temp2 = zero
        sum1 = zero
        do j=1, nconv
           
           sum1 = sum1+s(j,1)
           s(j,1) = sqrt(s(j,1))

           call av(v(1,j), ax)

           call dcopy(m, ax, 1, u(1,j), 1)
           temp = one/dnrm2(m, u(1,j), 1)
           call dscal(m, temp, u(1,j), 1)

           call daxpy(m, -s(j,1), u(1,j), 1, ax, 1)
           s(j,2) = dnrm2(m, ax, 1)**2

           call atv(u(1,j), ax)

           call daxpy(n, -s(j,1), v(1,j), 1, ax, 1)
           s(j,2) = s(j,2)+dnrm2(n, ax, 1)**2
           s(j,2) = sqrt(s(j,2)/TWO)
           temp2 = max(temp2,s(j,2))
           
        end do
           
        t1 = omp_get_wtime() - t1
        
        call dmout(6, nconv, 2, s, ncv, -6, &
             'Singular values and direct residuals')
        
        write (*,*) 'error',temp2
        
     end if
        
     if ( info .eq. 1) then
        print *, ' '
        print *, ' Maximum number of iterations reached.'
        print *, ' '
     else if ( info .eq. 3) then
        print *, ' ' 
        print *, ' No shifts could be applied during implicit', &
             ' Arnoldi update, try increasing NCV.'
        print *, ' '
     end if
     
     print *, ' '
     print *, ' _SVD '
     print *, ' ==== '
     print *, ' '
     print *, ' Size of the matrix is ', n
     print *, ' The number of Ritz values requested is ', nev
     print *, ' The number of Arnoldi vectors generated', &
          ' (NCV) is ', ncv
     print *, ' What portion of the spectrum: ', which
     print *, ' The number of converged Ritz values is ',  &
          nconv 
     print *, ' The number of Implicit Arnoldi update', &
          ' iterations taken is ', iparam(3)
     print *, ' The number of OP*x is ', iparam(9)
     print *, ' The convergence criterion is ', tol
     print *, ' Cumulative contribution ratio is ',  sum1/sum0
     print *, ' TIME MAXTRIX VECTOR ', s0
     print *, ' TIME ARPACK dsvd ', t1
     print *, ' '
     
  end if
  
end program dsvd_file
!!!
! a simple matrix-vector multiplications routine
!!!
Subroutine av_atv(xin, yout)
  Use my_data
  Use omp_lib
  Implicit None
  Integer i, j,thr_num
  Double precision, Dimension(n), Intent(in) :: xin
  Double precision, Dimension(n), Intent(out) :: yout
  Double precision  temp, zero
  parameter (zero=0.0d0)
  EXTERNAL           DDOT
  double precision   DDOT

  s1 = omp_get_wtime()
  if (mode .eq. 0) then

     yout=zero

!$OMP PARALLEL PRIVATE(i,thr_num,j,temp) REDUCTION(+:yout)
     thr_num = omp_get_thread_num()
     do i=start_row(thr_num+1),start_row(thr_num+2)-1
        temp = zero
        do j=IAP(i),IAP(i+1)-1
           temp=temp+(A(j))*xin(JA(j))
        enddo
        do j=IAP(i),IAP(i+1)-1
           yout(JA(j))=yout(JA(j))+(A(j))*temp
        enddo
     enddo
!$OMP END PARALLEL
     
  else

     yout=zero

!$OMP PARALLEL DO PRIVATE(temp) REDUCTION(+:yout)
     do i=1,M
        temp = ddot(n,DA(1,i),1,xin,1)
        call daxpy(n,temp,DA(1,i),1,yout,1)
     enddo
!$OMP END PARALLEL DO

  end if
  s0 = s0 + (omp_get_wtime() - s1)
  
  return
End Subroutine av_atv
!
Subroutine av(xin, yout)
  Use my_data
  Use omp_lib
  Implicit None
  Integer i, j,thr_num
  Double precision, Dimension(n), Intent(in) :: xin
  Double precision, Dimension(m), Intent(out) :: yout
  Double precision  zero
  parameter (zero=0.0d0)
  EXTERNAL           DDOT
  double precision   DDOT

  s1 = omp_get_wtime()
  if (mode .eq. 0) then

!$OMP PARALLEL PRIVATE(i,thr_num,j)
     thr_num = omp_get_thread_num()
     do i=start_row(thr_num+1),start_row(thr_num+2)-1
        yout(i) = zero
        do j=IAP(i),IAP(i+1)-1
           yout(i) = yout(i)+(A(j))*xin(JA(j))
        enddo
     enddo
!$OMP END PARALLEL
     
  else

!$OMP PARALLEL DO
     do i=1,M
        yout(i) = ddot(n,DA(1,i),1,xin,1)
     enddo
!$OMP END PARALLEL DO

  end if
  s0 = s0 + (omp_get_wtime() - s1)
  
  return
End Subroutine av
!
Subroutine atv(xin, yout)
  Use my_data
  Use omp_lib
  Implicit None
  Integer i, j,thr_num
  Double precision, Dimension(m), Intent(in) :: xin
  Double precision, Dimension(n), Intent(out) :: yout
  Double precision  zero
  parameter (zero=0.0d0)

  s1 = omp_get_wtime()
  if (mode .eq. 0) then

     yout=zero

!$OMP PARALLEL PRIVATE(i,thr_num,j) REDUCTION(+:yout)
     thr_num = omp_get_thread_num()
     do i=start_row(thr_num+1),start_row(thr_num+2)-1
        do j=IAP(i),IAP(i+1)-1
           yout(JA(j))=yout(JA(j))+(A(j))*xin(i)
        enddo
     enddo
!$OMP END PARALLEL
     
  else

     yout=zero

!$OMP PARALLEL DO REDUCTION(+:yout)
     do i=1,M
        call daxpy(n,xin(i),DA(1,i),1,yout,1)
     enddo
!$OMP END PARALLEL DO

  end if
  s0 = s0 + (omp_get_wtime() - s1)
  
  return
End Subroutine atv