!__________________________________________________________________________
! Documentation
!__________________________________________________________________________
! DONEV ALEKSANDAR--PHY301
! Determining the lowest energy eigenstate of a quantum anharmonic oscillator
! in a potential V=x^2+b*x^4.
! Schrodinger's equation (-d2/dx^2+V)psi=E[0]*psi with first order perturbation
! gives E[0]=1+3/4*b, for 0<b<<1. This should be  used for comparison.
!
! I use the banded iterative power method and Pentium optimized BLAS routines
! to find the smallest eigenvalue of Schrodinger's equation
! by cutting the potential off to the interval [-L,L] and setting
! up an equally spaced 1D lattice in this interval and solving
! Schrodinger's eigenvalue equation on it.
! The problem with this approach is the need to determine the cutoff L.
!__________________________________________________________________________
!__________________________________________________________________________

!__________________________________________________________________________
! Module holding the global precision variables
!__________________________________________________________________________
MODULE GlobalData

  ! We set the two different types of precision, single and double.
  ! These definitions agree with those of BLAS:
  INTEGER ,PARAMETER, PUBLIC :: sp=KIND(0.0E0), dp=KIND(0.0D0)
  ! In this run we use double precision:
  INTEGER ,PARAMETER, PUBLIC :: pr=dp
  ! These are specific values we need in the program:
  REAL(KIND=pr), PARAMETER, PUBLIC :: precision=0.000001_pr
  INTEGER, PARAMETER, PUBLIC :: max_iteration=1000

END MODULE GlobalData
!__________________________________________________________________________

!__________________________________________________________________________
! Module holding some useful matrix functions
!__________________________________________________________________________
MODULE MatrixFunctions

  USE GlobalData
  CONTAINS

  !_________________________________________________________
  ! Determining the dominant eigenvalue of a triangular symmetric or
  ! general banded matrix with n_diags superdiagonals using the
  ! power method and the BLAS routines GBMV and SBMV, which stand for:
  ! General/Symmetric Banded Matrix-Vector Multiply
  !_________________________________________________________
  SUBROUTINE PowerEigBanded(matrix,vector,lambda,matrix_type,n_diags)

    USE interface_gbmv, ONLY: gbmv
    USE interface_sbmv, ONLY: sbmv

    REAL(KIND=pr), INTENT(IN) :: matrix(:,:)
    REAL(KIND=pr), INTENT(INOUT) :: vector(:),lambda
    CHARACTER, INTENT(IN) :: matrix_type
    INTEGER, INTENT(IN), OPTIONAL :: n_diags

    REAL(KIND=pr) :: lambda_old,vector_copy(SIZE(vector))

    ! This is the actual iteration loop:
    Iterating: DO iteration=1,max_iteration
      ! We save the value of the eigenvalue:
      lambda_old=lambda
      ! We multiply the eigenvector with the matrix and do v<-normalized(Av)
      ! using either GBMV or SBMV depending on the matrix type:
      symmetric: SELECT CASE(matrix_type)
        CASE ('G')
          call GBMV(matrix,vector,vector_copy,n_diags)
        CASE ('S')
          call SBMV(matrix,vector,vector_copy)
        CASE DEFAULT
          STOP
      END SELECT symmetric

      lambda=SQRT(SUM(vector_copy**2))
      vector=vector_copy/lambda
      ! We stop the iteration if two successive estimates differ by precision:
      IF (ABS(lambda_old-lambda)<precision) THEN
        ! WRITE(*,*) "Successful convergence achieved at iteration",iteration
        EXIT Iterating
      END IF
    END DO Iterating

  END SUBROUTINE PowerEigBanded

  !_________________________________________________________
  ! Determining the smallest and largest eigenvalues
  ! of a general or symmetric banded matrix A.
  ! with POSITIVE ONLY eigenvalues using the power method for the largest eigenvalue.
  ! The smallest eigenvalue is in fact largest for the shifted matrix lambda_largest*I-A.
  !_________________________________________________________
  SUBROUTINE LowHighEigBanded(matrix,vector_low,vector_high,lambda_low,lambda_high,matrix_type,n_diags)

    REAL(KIND=pr), INTENT(IN) :: matrix(:,:)
    REAL(KIND=pr), INTENT(INOUT) :: vector_low(:),vector_high(:),lambda_low,lambda_high
    CHARACTER, INTENT(IN) :: matrix_type
    INTEGER, INTENT(IN), OPTIONAL :: n_diags

    ! The identity matrix (note: this is not the most efficient
    ! but it is the easiest to understand method.
    REAL(KIND=pr) :: identity(SIZE(matrix,1),SIZE(matrix,2))
    INTEGER :: i
    identity=0.0_pr

    ! We form the identity matrix and find the largest eigenvalues
    ! of both A and the shifted lambda_high*I-A.
    ! Where the diagonal of the identity is depends on the matrix type:
    symmetric: SELECT CASE(matrix_type)
      CASE ('G')
        identity(n_diags,:)=1.0_pr
        CALL PowerEigBanded(matrix,vector_high,lambda_high,'G',n_diags)
        CALL PowerEigBanded(lambda_high*identity-matrix,vector_low,lambda_low,'G',n_diags)
        lambda_low=lambda_high-lambda_low
      CASE ('S')
        identity(SIZE(matrix,1),:)=1.0_pr
        CALL PowerEigBanded(matrix,vector_high,lambda_high,'S')
        CALL PowerEigBanded(lambda_high*identity-matrix,vector_low,lambda_low,'S')
        lambda_low=lambda_high-lambda_low
      CASE DEFAULT
        STOP
    END SELECT symmetric

  END SUBROUTINE LowHighEigBanded

END MODULE MatrixFunctions
!__________________________________________________________________________

!__________________________________________________________________________
! Main program for determining the lowest energy state of
! a quantum oscillator in a well with potential V(x)
!__________________________________________________________________________
PROGRAM PowerMethod

  USE GlobalData
  USE MatrixFunctions
  USE Graphics

  ! This part simply declares all the variables along with several constants:
  INTEGER :: i,j,k,iteration,n,n_start,n_end,n_step
  REAL(KIND=pr),DIMENSION(:,:),ALLOCATABLE :: matrix
  REAL(KIND=pr),DIMENSION(:),ALLOCATABLE :: vector_low,vector_high,x_points
  REAL(KIND=pr) :: lambda_low=0.0_pr,lambda_high=0.0_pr,dx=0.0_pr
  REAL(KIND=pr) :: b,a,cutoff

  ! The user can enter various variables:
  WRITE(*,*) "Please enter the starting and ending number of points on the grid and the stride:"
  READ(*,*) n_start,n_end,n_step
  WRITE(*,*) "Please enter the cutoff length L and the perturbation b:"
  READ(*,*) cutoff,b

  WRITE(*,*) "Number of points=, smallest eigenvalue="
  WRITE(*,*) "_______________________________________"

  ! We change the number of points in the grid in a loop to see how E changes:
  NumberOfPoints: DO n=n_start,n_end,n_step

    IF(ALLOCATED(matrix)) DEALLOCATE(matrix,vector_low,vector_high,x_points)
    ALLOCATE(matrix(2,n),vector_low(n),vector_high(n),x_points(n))

    ! Now we initialize the grid for x with step dx:
    dx=2.0_pr*cutoff/(n+1)
    matrix=0.0_pr
    ! We assign the tridiagonal matrix for d2x/dx^2:
    FORALL(i=1:n)
      x_points(i)=(-cutoff+i*dx)
      matrix(2,i)=2.0_pr/dx**2+V(x_points(i))
      matrix(1,i)=-1/dx**2
    END FORALL

    ! We initialize the eigenvectors to random vectors:
    CALL RANDOM_NUMBER(vector_low)
    CALL RANDOM_NUMBER(vector_high)

    ! We calculate the smallest and largest eigenvalues of this symmetric matrix:
    CALL LowHighEigBanded(matrix,vector_low,vector_high,lambda_low,lambda_high,'S')

    !We print the results:
    WRITE(*,*) n,lambda_low

  END DO NumberOfPoints

  ! Now we plot the eigenvector for the largest number of grid points
  ! using the DISLIN library and my module Graphics:
  CALL InitGraphics("QHO.tif","TIFF","QSHO lowest eigenstate","x","|psi(x)|^2")
  CALL Plot2D(REAL(x_points,sp),REAL(vector_low,sp),"L-R",REAL(x_points,sp),REAL(vector_low,sp),"STB")
  CALL EndGraphics()

  CONTAINS
  !_________________________________________________________
  ! The potential is given as an internal function
  !_________________________________________________________
  PURE REAL(KIND=dp) FUNCTION V(x)
    REAL(KIND=dp), INTENT(IN) :: x
    V=x**2+b*(x**4)
  END FUNCTION V

END PROGRAM PowerMethod
!__________________________________________________________________________