program timed
C
C     by Petr Bour, Prague 1994
C
C     This subroutine finds eigenvectors and eigenvalues
C     in a frequency band for a matrix A (in file A3.SCR)
C     Diagonal terms are already divided by two in A
C     M is soft dimension of the subspace where direct diagonalization is done
C     NE is number of eigenvectors to be found
C     MV is the number of eigenvectors found
C
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      PARAMETER (NM0=99999,fp=0.1d0)
      DIMENSION axi(NM0),xi(NM0),xr(NM0),xw(NM0),axr(NM0),
     1j33(NM0),a33(NM0),tem(1000),is(50)
      logical lex
c
      do i=1,50
       is(i)=0
      enddo
c
      write(6,*)'Timed starts'
      inquire(file='A3.SCR',exist=lex)
      if(.not.lex)then
       write(6,*)'A3.SCR not found'
       stop
      endif
c
      write(6,*)'Dimension:'
      read(5,*)n
      if(n.gt.NM0)then
       write(6,*)'dimension too high'
       stop
      endif
      write(6,*)'Frequency band (wmin wmax):'
      read(5,*)wmin,wmax
      write(6,*)' Number of frequency samples nw:'
      read(5,*)nw
      write(6,*)'Number of time steps nsteps:'
      read(5,*)nsteps
      dt=2.0d0/(wmax+wmin)*fp
      write(6,6000)dt
6000  format('Time step: ',f10.5)
c
c     initialize time vectors
      do 1 i=1,nw
      fac=1.0d0/sqrt(real(n))
      do 1 j=1,n
      xr(j)=fac
1     xi(j)=0.0d0
c
c     initialize and save w vectors (real)
      open(48,file='VECTORW.SCR',form='unformatted')
      do 7 i=0,nw
      do 2 j=1,n
2     xw(j)=0.0d0
7     write(48)(xw(j),j=1,n)
      close(48)
c
c     start time propagation
c     ============================================================
      t=0.0d0
      do 3 it=1,nsteps
      t=t+dt
c
c     read A and form a sum X+iAX*dt=X(t+d)
      do 4 i=1,n
      axr(i)=xr(i)
4     axi(i)=xi(i)
      open(49,file='A3.SCR',form='unformatted')
      DO 204 I=1,n
      read(49)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn)
      DO 204 ii=1,nn
      J=j33(ii)
      a=a33(ii)*dt
      axi(I)=axi(I)+a*xr(J)
      axi(J)=axi(J)+a*xr(I)
      axr(I)=axr(I)-a*xi(J)
204   axr(J)=axr(J)-a*xi(I)
      close(49)
c
c     renormalize X(t+dt) and write to xr,xi
      an=0.0d0
      do 8 i=1,n
8     an=an+axi(i)**2+axr(i)**2
      fac=1.0d0/sqrt(an)
      do 9 i=1,n
      xi(i)=axi(i)*fac
9     xr(i)=axr(i)*fac
c
c     update w vectors
      open(48,file='VECTORW.SCR',form='unformatted')
      do 5 iw=0,nw
      w=wmin+(wmax-wmin)/nw*iw
      read(48)(xw(j),j=1,n)
      c=cos(w*t)
      s=sin(w*t)
      do 6 i=1,n
6     xw(i)=xw(i)+c*xr(i)+s*xi(i)
c
c     find norm for this w:
      an=0.0d0
      do 10 i=1,n
10    an=an+xw(i)**2
      tem(iw+1)=an
c
      backspace 48
5     write(48)(xw(j),j=1,n)
      close(48)
c
      ax=0
      do 11 iw=0,nw
11    if(tem(iw+1).gt.ax)ax=tem(iw+1)
      do 12 iw=0,nw
      w=wmin+(wmax-wmin)/nw*iw
      item=min(nint(50.0d0*tem(iw+1)/ax),50)
12    write(6,6002)w,(is(i),i=1,item)
6002  format(f10.2,50i1)
      write(6,*)it
c
3     continue
c     ============================================================
c
      END
c