      program lanczos
*---------------------------------------------------------------
c
c     Lanczos-Verfahren für eine hermitesche Matrix
c     Das Programm bestimmt die Dimension der Matrix selbständig und
c     geht dabei davon aus, dass die letzte Zeile bzw. Spalte besetzt ist.
c
c     UNVOLLSTÄNDIG
*---------------------------------------------------------------
      implicit   none
      integer*8  i,k,L,Dim      ! Dim = max. zu erwartende Dimensionen
      integer*8  ThisDim        ! Dimension der urspr. Matrix
      integer*8  ListLen        ! Dimension der Liste der Matrixelemente
      Parameter  (Dim=10000)
*
      integer*8  kList(Dim),LList(Dim)
      complex*16 HkL,HkLList(Dim)
*     
      complex*16 LV(Dim,3)
      integer*8  LVDim          ! Zahl der Lanczos-Schritte
      Parameter  (LVDim=100)
      complex*16 LanczosMatrix(LVDim,LVDim)
      real*8     Norm
      complex*16 Skalarprodukt
*
      integer    wkDim,rkDim,Error
      Parameter  (wkDim=3*LVDim+1,rkDim=3*LVDim+1)
      complex*16 wk(wkDim)
      real*8     Eigenwerte(LVDim)
      real*8     rk(rkDim)
*     
      character  Filename*50
      character  Version*50
*
      Version= '# lanczos 01 Juergen Schnack  12. 06. 2009'
*---------------------------------------------------------------
*
*     Version
*
      write(6,'(a)') Version
      LanczosMatrix= DCMPLX(0.d0,0.d0)
*---------------------------------------------------------------
*
*     Einlesen
*
      ListLen = 0
      ThisDim = 0
      call GETARG(1,FileName)
      open(1,file=Filename)
 100  read(1,*,err=110,end=110) k,L,Hkl
      IF(k.gt.ThisDim) ThisDim=k
      IF(L.gt.ThisDim) ThisDim=L
      ListLen = ListLen + 1     
      kList(ListLen) = k
      LList(ListLen) = L
      HkLList(ListLen) = HkL
      goto 100
 110  continue
      close(1)
      write(6,*) 'Dimension der Matrix = ',ThisDim
      write(6,*) 'Laenge der Liste     = ',ListLen
*---------------------------------------------------------------
*
*     1. Lanczos-Vektor definieren
*
      LV = DCMPLX(0.d0,0.d0)

      do k=1,ThisDim
         LV(k,1) = DCMPLX(DSIN(2.718281828d0*k),DCOS(2.718281828d0/k))
      enddo

*     oder
c      LV(1,1) = DCMPLX(1.d0,0.d0) ! nur 1. Komponente ungleich Null

*---------------------------------------------------------------
*
*     1. Lanczos-Vektor normieren
*
      Norm = 0.d0
      do k=1,ThisDim
         Norm = Norm + CONJG(LV(k,1))*LV(k,1)
      enddo
      Norm = DSQRT(Norm)
      write(6,*) 'Norm des ersten Lanczos-Vektors (vorher) = ',Norm
*---------------------------------------------------------------
      do i=1,ThisDim
         LV(i,1) = LV(i,1)/Norm
      enddo
*---------------------------------------------------------------
      Norm = 0.d0
      do k=1,ThisDim
         Norm = Norm + CONJG(LV(k,1))*LV(k,1)
      enddo
      Norm = DSQRT(Norm)
      write(6,*) 'Norm des ersten Lanczos-Vektors (nachher) = ',Norm

*---------------------------------------------------------------
*
*     2. Lanczos-Vektor ausrechnen
*
      do i=1,ListLen
         k = kList(i)
         L = LList(i)
         HkL = HkLList(i)
         LV(k,2) = LV(k,2) + HkL * LV(L,1)
      enddo
*---------------------------------------------------------------
*
*     Skalarprodukt mit der ersten <1|2>; Vorsicht komplexes Skalarprodukt
*
      Skalarprodukt = DCMPLX(0.d0,0.d0)
      do k=1,ThisDim
         Skalarprodukt = Skalarprodukt + CONJG(LV(k,1))*LV(k,2)
      enddo
      write(6,*) 'Skalarprodukt zwischen 1 und 2 = ',Skalarprodukt

      LanczosMatrix(1,1) = Skalarprodukt
      
*---------------------------------------------------------------
*
*     |1> aus |2> rausprojizieren
*
      do i=1,ThisDim
         LV(i,2) = LV(i,2) - Skalarprodukt * LV(i,1)
      enddo
*---------------------------------------------------------------
*
*     2. Lanczos-Vektor normieren
*
      Norm = 0.d0
      do k=1,ThisDim
         Norm = Norm + CONJG(LV(k,2))*LV(k,2)
      enddo
      Norm = DSQRT(Norm)
      write(6,*) 'Norm des zweiten Lanczos-Vektors (vorher) = ',Norm

*---------------------------------------------------------------
      do i=1,ThisDim
         LV(i,2) = LV(i,2)/Norm
      enddo
*---------------------------------------------------------------
      Norm = 0.d0
      do k=1,ThisDim
         Norm = Norm + CONJG(LV(k,2))*LV(k,2)
      enddo
      Norm = DSQRT(Norm)
      write(6,*) 'Norm des zweiten Lanczos-Vektors (nachher) = ',Norm
*---------------------------------------------------------------
*
*     Zusatzaufgabe : ab drei geht es ja allgemein weiter ...
*








*---------------------------------------------------------------
*
      end
