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

program dpca_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(:)
!
!     %---------------%
!     | Local Scalars |
!     %---------------%
!
  CHARACTER*1000    ARG1, ARG2, ARG3, ARG4
  character        bmat*1, which*2
  integer          ido, nev, ncv, lworkl, info, ierr, &
       accuracy, 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
!
!     %------------%
!     | Parameters |
!     %------------%
!
  Double precision two, one, zero
  parameter        (two = 2.0D+0, one = 1.0D+0, zero = 0.0D+0)
!  
!     %-----------------------------%
!     | BLAS & LAPACK routines used |
!     %-----------------------------%
!
  Double precision dnrm2
  external         dnrm2, daxpy, dcopy, dscal
!
!     %-----------------------%
!     | 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)
  
  READ (ARG1,*) nev
  READ (ARG2,*) ncv 

  WRITE (*,*) nev, ncv

  READ (ARG3,*) accuracy

  tol = 10.0D0**(-accuracy)
      
  open(18,file=ARG4,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) )

  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
  ENDIF
  
  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
  
  close(18)

  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
  endif
  
  do i=1,maxthreads+1
     write (*,*) i,start_row(i)
  end do
  
  ALLOCATE( ss0(n), cc0(n), sx(n) )

  cc0=zero
         
!$OMP PARALLEL DO PRIVATE(j) REDUCTION(+:cc0)
  do i=1,m
     do j=IAP(i),IAP(i+1)-1
        cc0(JA(j))=cc0(JA(j))+(A(j))
     enddo
  enddo
!$OMP END PARALLEL DO
  
  cc0=cc0/DBLE(m)
  
  ss0=zero

!$OMP PARALLEL DO PRIVATE(j) REDUCTION(+:ss0)
  do i=1,m
     do j=IAP(i),IAP(i+1)-1
        ss0(JA(j))=ss0(JA(j))+(A(j)-cc0(JA(j)))**2
     enddo
  enddo
!$OMP END PARALLEL DO

  ss0(1:n)=sqrt(ss0(1:n))

  do i=1,n
     if (ss0(i) .eq. zero) write (*,*) 'not colm full rank'
  enddo

  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 dpca_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, csx, cx
  parameter (zero=0.0d0)
  external           ddot
  double precision   DDOT

  s1 = omp_get_wtime()
      
  do i=1,N
     sx(i)=xin(i)/ss0(i)
  enddo
  
  csx=ddot(N,cc0,1,sx,1)
  
  yout=zero
  
  cx=zero

!$OMP PARALLEL PRIVATE(i,thr_num,j,temp) REDUCTION(+:yout) REDUCTION(+:cx)
  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))*sx(JA(j))
     enddo
     temp=temp-csx
     do j=IAP(i),IAP(i+1)-1
        yout(JA(j))=yout(JA(j))+(A(j))*temp
     enddo
     cx = cx + temp
  enddo
!$OMP END PARALLEL

  call daxpy(N,-cx,cc0,1,yout,1)
  
  do i=1,N
     yout(i)=yout(i)/ss0(i)
  enddo
  
  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,csx
  parameter (zero=0.0d0)
  external           ddot
  double precision   DDOT

  s1 = omp_get_wtime()
  
  do i=1,N
     sx(i)=xin(i)/ss0(i)
  enddo

  csx=ddot(N,cc0,1,sx,1)

!$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))*sx(JA(j))
     enddo
     yout(i)=yout(i)-csx
  enddo
!$OMP END PARALLEL

  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,cx
  parameter (zero=0.0d0)

  s1 = omp_get_wtime()
  
  yout=zero

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

  do i=1,N
     yout(i)=yout(i)/ss0(i)
  enddo

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