	PROGRAM lanczos
!	driver for routine lanczos (c.f. Grosso's book)
	USE nrtype
	USE nr
	IMPLICIT NONE
	INTEGER(I4B), PARAMETER :: NP=2
	REAL(DP), PARAMETER :: TINY=1.0e-6_dp
	INTEGER(I4B) :: i,j
	REAL(DP), DIMENSION(NP) :: d,e,f
!REAL(DP), DIMENSION(NP,NP) :: a,c = reshape( (/ &
!!5.0_dp,4.3_dp,3.0_dp,2.0_dp,1.0_dp,&
!!0.0_dp,-1.0_dp,-2.0_dp,-3.0_dp,-4.0_dp,&
!!4.3_dp,5.1_dp,4.0_dp,3.0_dp,2.0_dp,&
!!1.0_dp,0.0_dp,-1.0_dp,-2.0_dp,-3.0_dp,&
!!3.0_dp,4.0_dp,5.0_dp,4.0_dp,3.0_dp,&
!!2.0_dp,1.0_dp,0.0_dp,-1.0_dp,-2.0_dp,&
!!2.0_dp,3.0_dp,4.0_dp,5.0_dp,4.0_dp,&
!!3.0_dp,2.0_dp,1.0_dp,0.0_dp,-1.0_dp,&
!!1.0_dp,2.0_dp,3.0_dp,4.0_dp,5.0_dp,&
!!4.0_dp,3.0_dp,2.0_dp,1.0_dp,0.0_dp,&
!!0.0_dp,1.0_dp,2.0_dp,3.0_dp,4.0_dp,&
!!5.0_dp,4.0_dp,3.0_dp,2.0_dp,1.0_dp,&
!!-1.0_dp,0.0_dp,1.0_dp,2.0_dp,3.0_dp,&
!!4.0_dp,5.0_dp,4.0_dp,3.0_dp,2.0_dp,&
!!-2.0_dp,-1.0_dp,0.0_dp,1.0_dp,2.0_dp,&
!!3.0_dp,4.0_dp,5.0_dp,4.0_dp,3.0_dp,&
!!-3.0_dp,-2.0_dp,-1.0_dp,0.0_dp,1.0_dp,&
!!2.0_dp,3.0_dp,4.0_dp,5.0_dp,4.0_dp,&
!!-4.0_dp,-3.0_dp,-2.0_dp,-1.0_dp,0.0_dp,&
!!1.0_dp,2.0_dp,3.0_dp,4.0_dp,5.0_dp /),&
!!(/ NP,NP /) )
        REAL(DP), DIMENSION(NP,NP) :: a,c = reshape( (/ &
        2.0_dp,1._dp,1._dp,0.0_dp /),  &
        (/ NP,NP /) )
	a(:,:)=c(:,:)
	a(:,:)=c(:,:)
	call lanc(a,d,e)
	call tqli(d,e,a)
	write(*,'(/1x,a)') 'Eigenvectors for a real symmetric matrix'
	do i=1,NP
		f=matmul(c(:,:),a(:,i))
		write(*,'(/1x,a,i3,a,f10.6)') 'Eigenvalue',i,' =',d(i)
		write(*,'(/1x,t7,a,t17,a,t31,a)') 'Vector','Mtrx*Vect.','Ratio'
		do j=1,NP
			if (abs(a(j,i)) < TINY) then
				write(*,'(1x,2f12.6,a12)') a(j,i),f(j),'div. by 0'
			else
				write(*,'(1x,2f12.6,e14.6)') a(j,i),f(j),&
					f(j)/a(j,i)
			end if
		end do
		write(*,'(/1x,a)') 'press ENTER to continue...'
		read(*,*)
	end do
	END PROGRAM lanczos


subroutine lanc(H,d,e)
   use nrtype
   implicit none
   INTEGER(I4B), PARAMETER :: NP=2
   REAL(DP), DIMENSION(NP),intent(out) :: d,e
   REAL(DP), DIMENSION(NP,NP),intent(in) :: H
   REAL(DP), DIMENSION(NP) :: fm,fmm,g 
   real(dp)::a,b
   integer::j
   ! initial seed vector, normalized
   j=1;
   !call random_number(fm); fm=fm/sqrt(dot_product(fm,fm))
   fm(1)=1.0; fm(2)=1.0; fm=fm/sqrt(dot_product(fm,fm))
   a= dot_product(fm,matmul(H,fm)); b=0.0
   d(j)=a; e(j)=b; 
   j=2;
   g=matmul(H,fm)-a*fm;
   b=sqrt(dot_product(g,g)); e(j)=b;
   fmm=fm;
   fm=g/b;
   a=dot_product(fm,matmul(H,fm))
   d(j)=a;
   do j=3,NP
      g=matmul(H,fm)-a*fm-b*fmm;
      b=sqrt(dot_product(g,g)); e(j)=b;
      fmm=fm
      fm=g/b
      a=dot_product(fm,matmul(H,fm))
      d(j)=a; 
   enddo
end subroutine
