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 *--------------------------------------------------------------- !$omp parallel private(i) shared(LV,Norm) !$omp do do i=1,ThisDim LV(i,1) = LV(i,1)/Norm enddo !$omp end parallel *--------------------------------------------------------------- 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 * !$omp parallel private(i) shared(LV,Skalarprodukt) !$omp do do i=1,ThisDim LV(i,2) = LV(i,2) - Skalarprodukt * LV(i,1) enddo !$omp end parallel *--------------------------------------------------------------- * * 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 *--------------------------------------------------------------- !$omp parallel private(i) shared(LV,Norm) !$omp do do i=1,ThisDim LV(i,2) = LV(i,2)/Norm enddo !$omp end parallel *--------------------------------------------------------------- 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