text:
!===================================================================2
!  LANCZOS routines for BIGSTICK
!
!  versions for 'new' parallelization scheme -- FALL 2011
!
!  This code uses LAPACK ROUTINES
!
!  LAPACK copyright statements and license
!
!Copyright (c) 1992-2013 The University of Tennessee and The University
!                        of Tennessee Research Foundation.  All rights
!                        reserved.
!Copyright (c) 2000-2013 The University of California Berkeley. All
!                        rights reserved.
!Copyright (c) 2006-2013 The University of Colorado Denver.  All rights
!                        reserved.

!Additional copyrights may follow

!Redistribution and use in source and binary forms, with or without
!modification, are permitted provided that the following conditions are
!met:
!
!- Redistributions of source code must retain the above copyright
!  notice, this list of conditions and the following disclaimer.
!
!- Redistributions in binary form must reproduce the above copyright
!  notice, this list of conditions and the following disclaimer listed
!  in this license in the documentation and/or other materials
!  provided with the distribution.
!
!- Neither the name of the copyright holders nor the names of its
!  contributors may be used to endorse or promote products derived from
!  this software without specific prior written permission.
!
!The copyright holders provide no reassurances that the source code
!provided does not infringe any patent, copyright, or any other
!intellectual property rights of third parties.  The copyright holders
!disclaim any liability to any recipient for claims brought against
!recipient by any third party for infringement of that parties
!intellectual property rights.
!
!THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
!LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
!A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
!OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
!SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
!LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
!DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
!THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
!(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
!OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!===== END OF LACK COPYRIGHT STATEMENT ========================

module lanczos_util

	use lanczos_info
	use precisions

	implicit none

    integer startiter   ! # of starting iterations already

    real(kind=8), allocatable :: alpha(:),beta(:)  ! lanczos coefficients
    real(kind=egv_prec), allocatable ::  eiglvec(:,:)
    real(kind=egv_prec), allocatable ::  eiglvectx(:,:)
    real(kind=egv_prec), allocatable,target ::  e(:),etx(:)

    integer :: txchoice = 1   ! ADDED in 7.9.0

    logical :: flyeig
    integer :: nthick_add   ! # of extra states to add for thick-restart; must be >= ncheck_ex

  !--------------------- CONVERGENCE -----------------------------------------

    logical :: finished
    integer :: ncheck
    real(kind=egv_prec), allocatable :: eold(:)  ! old energies
    real    :: ediff0,ediff1

  !-------------- ALTERNATE CONVERGENCE ---------------------------
  !               based on wavefunctions
    logical :: altconverge

    real(kind=egv_prec), allocatable :: eiglvecold(:,:)
!--------------------------- GREEN'S FUNCTION-------------

	real :: Egf   !Energy in Green function 1/(Egf-H)
	real :: green_tol ! tolerance for convergence of green's function
	real(kind=8), allocatable :: gvec(:),gvec_old(:)

contains

!===========================================================================
!  SUBROUTINES IN THIS FILE
!   lanczos_menu
!	initialize_lanczos_vector
!       swap_vchar               : controls "direction" of mat-vec multiply
!                                   either H * vec1 = vec2  or  H * vec2 = vec1
!       lanczos_p                : main lanczos subroutine
!       initialize_final         : zeroes out next vector to be calculated
!       lanczos_output
!       find_lanczos_eigenvalues
!       density1b_output
!       exactdiag_p
!       thick_restart_sub_p
!       random_restart_p
!       open_lanczosfile
!       close_lanczosfile
!================================================
!
!  routine to find eigenvalues of lanczos matrix
!  if vecflag = .true. then find eigenvectors too
!  if nthick > 0  then arrange according to thick restart
!  nwindow = # of eigenvalues desired; only needed for targeted thick-restart
!
!  MODIFIED 7.6.8 so one can choose a target energy; REMOVED 7.7.4
!  MODIFIED 7.8.7 to allow for block lanczos
!
! CALLED BY:
!    lanczos_p
!    lanczos_block
!    lanczos_output
!    thick_restart_sub_p
!

subroutine find_lanczos_eigenvalues(n,vecflag,nthick)

  use precisions
  use lanczos_info,only:Etarget,lanczchar
  use localblocks
  use mod_balg
  implicit none

  integer(4):: n  ! actual dimension
!  integer :: lwork

!  real(kind = 4) :: alpha(np), beta(np)
!  real(kind = egv_prec) :: e(np), vec(np,np)
  logical :: vecflag  ! flag to get eigenvectors or not
  integer :: nthick
  real(kind = egv_prec) ::  work3(3*niter)  ! work(np)
  integer i,j,k
!  integer(4) :: nwindow, thickshift
  real(kind = egv_prec) :: etmp
  real(kind=egv_prec) ::  ediff(n),ediff0,xswap

  integer :: itarget

  integer info

if(block_flag == 0 .and. lanczchar/='tx')then
   eiglvec(:,:) = 0.0
   do i = 1,n
      eiglvec(i,i) = real(alpha(i),kind=egv_prec)
   enddo ! i
!-------------- put in column of thick-restart matrix elements

   if(nthick > 0)then
      do i = 1,nthick
         eiglvec(i,nthick+1) = real(beta(i),kind=egv_prec)
         eiglvec(nthick+1,i) = real(beta(i),kind=egv_prec)
      enddo
   endif
!-------- what about beta(nthick+1?)
   do i = nthick+1,n-1
      eiglvec(i,i+1) = real(beta(i),kind=egv_prec)
      eiglvec(i+1,i) = real(beta(i),kind=egv_prec)
   enddo
else   ! READ IN BLOCK STRUCTURE
	do i = 1,n
		do j = 1,n
			eiglvec(i,j)=lanczos_mat(i,j)
!			if(lanczchar=='tx')eiglvectx(i,j)=lanczos_mat(i,j)
		end do
	end do

end if

  call clocker('egv','sta')
  if(vecflag)then
      if(egv_prec==4)then
        call SSYEV( 'V','U', N, eiglvec, Niter, e, WORK3, 3*niter, INFO )
      else
        call DSYEV( 'V','U', N, eiglvec, Niter, e, WORK3, 3*niter, INFO )
      end if
  else
      if(egv_prec==4)then

         call SSYEV( 'N','U', N, eiglvec, Niter, e, WORK3, 3*niter, INFO )
      else
         call DSYEV( 'N','U', N, eiglvec, Niter, e, WORK3, 3*niter, INFO )

      end if
  endif
  call clocker('egv','end')

  if(kpow > 1 .and. .not. vecflag)then   ! only do for convergence test
	  do i =1,n
		  if(mod(kpow,2)==0)then
			  e(i)= e(i)**(1./float(kpow))
		  else
			  xswap=e(i) ! HERE xswap IS JUST A CONVENIENT DUMMY
			  e(i)=(abs(e(i)))**(1./float(kpow))
			  if(xswap < 0)e(i)=-e(i)
		  end if

	  end do

  end if
!--------- OPTION TO CHOOSE EXCITED STATES --------------------
!  added in 7.7.8 based in part upon work by R. Zbikowski
!  SORT EIGENVALUES, EIGENVECTORS BASED UPON PROXIMITY TO Etarget
!
  if(lanczchar=='tx')then

if(txchoice==1)then
!................  SET UP DIFFERENCES OF ENERGIES .............
  	do i = 1,n
  		ediff(i)= abs(e(i)-etarget)
  	end do
!........... NOW SIMPLE BUBBLE SORT .......................
      do i = 1,N-1
  		k = i
  		ediff0 = ediff(i)
  		do j = i+1,N
  			if(ediff(j) < ediff0)then
  				k=j
  				ediff0 = ediff(j)
  			end if
  		end do
  		if(k /= i)then
  			ediff(k)=ediff(i)
  			ediff(i)=ediff0
  			xswap   =e(k)
  			e(k)    =e(i)
  			e(i)    =xswap
  			do j = 1,N
  				xswap = eiglvec(j,k)
  				eiglvec(j,k)=eiglvec(j,i)
  				eiglvec(j,i)=xswap
  			end do
  		end if
      end do

!.......... DO A SECOND SORT JUST ON THE Nkeep "LOWEST"........
!          THIS PUTS EVERYTHING INTO A "PRETTY" ORDER
do i = 1,Nkeep-1
   k = i
   ediff0 = e(i)
   do j = i+1,Nkeep
	   	if(e(j) < ediff0)then
		   k=j
		   ediff0 = e(j)
	    end if
   end do
   if(k /= i)then
 	  e(k)=e(i)
	  e(i)=ediff0
 	  do j = 1,N
		xswap = eiglvec(j,k)
		eiglvec(j,k)=eiglvec(j,i)
		eiglvec(j,i)=xswap
	  end do
   end if
end do
!.............. instead, solve for (H-E0)^2
!   ADDED 7.9.9,
else  ! txchoice = 2
	do i = 1,n
		lanczos_mat(i,i)=lanczos_mat(i,i)-etarget
	end do
	eiglvectx= 0.0
   do i = 1,n
	   do j = 1,n
		   do k = 1,n
			   eiglvectx(i,j)=eiglvectx(i,j)+ lanczos_mat(i,k)*lanczos_mat(k,j)

		   end do
	   end do
   end do
   do i = 1,n
	   lanczos_mat(i,i)=lanczos_mat(i,i)+etarget
   end do
   if(vecflag)then
       if(egv_prec==4)then
         call SSYEV( 'V','U', N-1, eiglvectx, Niter, etx, WORK3, 3*niter, INFO )
       else
         call DSYEV( 'V','U', N-1, eiglvectx, Niter, etx, WORK3, 3*niter, INFO )
       end if
   else
       if(egv_prec==4)then

          call SSYEV( 'N','U', N-1, eiglvectx, Niter, etx, WORK3, 3*niter, INFO )
       else
          call DSYEV( 'N','U', N-1, eiglvectx, Niter, etx, WORK3, 3*niter, INFO )

       end if
   endif
   do i = 1,n-1
	   etx(i)=sqrt(abs(etx(i)))
   end do
!   print*,' in alternate ',e(1:10)

end if

  end if

  return
end subroutine find_lanczos_eigenvalues

!====================================================================
!
! routine to control construction and "improvement" of initial vector
! added in 7.6.3 by CWJ @ SDSU
!
!  CALLED BY:  lanczos_menu
!  SUBROUTINES CALLED:
!  br_open_u
!  setup_localvectors
!  intialize_lanczos_vector
!
subroutine setup_for_lanczos
	use flagger
    use mod_reorthog
	use io
	use bvectorlib_mod
    implicit none

    call setup_localvectors
    call initialize_lanczos_vector

    return

end subroutine setup_for_lanczos
!*----------------------------------------------------------------
!   subroutine read_pivot_block
!
!   Purpose: This is just the modified readpivot subroutine from
!   bwfnlib.f90. The routine loads in a number of pivot vectors into
!   history equal to dimblock
!
!   the file of eigenvectors must be previously opened
!   and the header reader;
!   reads the # of kept vectors (nkeep),
!   prints out E J, T^2 of each vector
!   then asks which to use as a pivot
!
!   OUTPUT
!      v(n), n = 1,dimbasis = vector read from disk
!      dnorm = norm of pivot vector
!
!   CALLS
!     wnf_read_nkeep
!     wfn_readeigenvec
!     wfn_close_file
!
!*----------------------------------------------------------------
subroutine read_pivot_block
   use localvectors
   use localblocks
   use precisions
   use basis
   use io
   use nodeinfo
   use bvectorlib_mod
   use mod_reorthog
   use wfn_mod
   implicit none

   real (kind = 8) :: dnorm
   integer nnkeep
   integer ikeep
   integer(4):: i,j,n, ierr
   real :: e,xj,xt2
   logical zeroflag
   integer :: aerr
   real(kind=8) :: scale

   if(iproc==0)print*, " reading in pivot block"

   call wfn_read_nkeep(oldwfnfile, nnkeep)  ! does BCAST

   if(nnkeep < dimblock)then
	   if(iproc==0)print*,' Insufficient pivot vectors, only ',nnkeep,', asked for ',dimblock
	   if(iproc==0)print*,' Block size must be <= ',nnkeep

	   call BMPI_ABORT(icomm,101,ierr)
	   stop
   end if

   do ikeep = 1, dimblock
!    call BMPI_BCAST(ikeep,1,0,icomm,ierr)
      vec1 = 0.d0
      call wfn_readeigenvec(oldwfnfile, frag1, fcomm1, vec1, ikeep, e, xj, xt2) ! KSM: updated
      call br_grab_vec1()! push vec1 -> br_reg
! DON'T NORMALIZE--THIS COMES LATER
      !call br_normalize(scale)
      call br_add2hist(1)
!write(logfile,*) "------------history matrix--------------"
!write(logfile,*) br_histbuf(:,:)
   enddo

   call wfn_close_file(oldwfnfile)

   !call exit
!print*, " reading in pivot block done"
   return
end subroutine read_pivot_block
!
!*----------------------------------------------------------------
!*----------------------------------------------------------------
!   subroutine read_pivot_block_chose
!
!   modified 7.9.4 from read_pivot_block
!
!   Purpose: This is just the modified readpivot subroutine from
!   bwfnlib.f90. The routine loads in a number of pivot vectors into
!   history equal to dimblock
!
!
!
!
!   OUTPUT
!      v(n), n = 1,dimbasis = vector read from disk
!      dnorm = norm of pivot vector
!
!   CALLS
!     wnf_read_nkeep
!     wfn_readeigenvec
!     wfn_close_file
!
!*----------------------------------------------------------------
subroutine read_pivot_block_choose
   use localvectors
   use localblocks
   use precisions
   use basis
   use io
   use nodeinfo
   use bvectorlib_mod
   use mod_reorthog
   use wfn_mod
   use mod_balg
   implicit none

   real (kind = 8) :: dnorm
   integer nnkeep
   integer ikeep
   integer(4):: i,j,n, ierr
   real :: e,xj,xt2
   logical zeroflag
   integer :: aerr
   real(kind=8) :: scale

   integer :: ilist

   if(iproc==0)print*, " reading in pivot block"
   allocate(pivot_choice(dimblock))
   allocate(dnormblock(dimblock))

   call wfn_read_nkeep(oldwfnfile, nnkeep)  ! does BCAST

   if(nnkeep < dimblock)then
	   if(iproc==0)then
		   print*,' Insufficient pivot vectors, only ',nnkeep,', asked for ',dimblock
		   print*,' (NEW)  Fill out rest with random vectors '
	   end if

!	   if(iproc==0)print*,' Block size must be <= ',nnkeep
!	   call BMPI_ABORT(icomm,101,ierr)
!	   stop
   end if

   if(iproc==0) print*,' There are ', nnkeep,' wavefunctions '

   if(iproc==0)then
        print*,nkeep,' states '
        print*,' '
        print*,' STATE      E         J          <H > '
	endif

    do ikeep = 1,nnkeep
      i = ikeep
      ! new interface, we say which vector to read - it checks
      call wfn_readeigenvec(oldwfnfile, frag1, fcomm1, vec1,i,e,xj,xt2) ! KSM: updated
      if(iproc==0)then
          write(6,101)i,e,xj,xt2

 101      format(i5,2x,4(1x,f11.4))
      end if
   end do ! ikeep

   if(iproc==0)then
	   print*,' '
	   print*,' Enter a list of ',dimblock,' states for the pivot block'
	   print*,' (Please enter in order; enter 0 to create random )'  ! modified in 7.9.11
	   if(dimblock > nnkeep)print*,' At least ',dimblock-nnkeep,' must be random '

	   if(auto_input)then
		   read(autoinputfile,*)pivot_choice(:)
		   print*,' choosing states ',pivot_choice
		   write(logfile,*)' Choosing pivot block states ',pivot_choice
	   else
 	      read(5,*)pivot_choice(:)
	      write(logfile,*)' Choosing pivot block states ',pivot_choice
		  write(autoinputfile,*)pivot_choice(:)
	  end if
   end if

   call BMPI_BCAST(pivot_choice,dimblock,0,icomm,ierr)
   call BMPI_BARRIER(icomm,ierr)

   do ilist = 1,dimblock
	   ikeep= pivot_choice(ilist)
!   if(iproc==0)call BMPI_BCAST(ikeep,1,0,icomm,ierr)
!	   pivot_choice(ilist)=ikeep

	   if(ikeep > nnkeep)then
           if(iproc==0)then
			   print *, "Selected vector is not in range (A):", nnkeep
			   print*,' Will replace with random vector '
		   end if
		   pivot_choice(ilist)=0
!		   stop 1
	   end if

   end do
   if(iproc==0)print*,pivot_choice,' chosen '

   do ilist = 1,dimblock
!    call BMPI_BCAST(ikeep,1,0,icomm,ierr)
      vec1 = 0.d0
	  ikeep = pivot_choice(ilist)

	  if(ikeep < 1)then ! create a random vector
	  	call add_pivot_br_reg(0,ilist)
	  	call br_normalize(dnorm)
	  	call br_orthogonalize(0)
	  	call br_normalize(dnorm)
	    call br_add2hist(1)
	  else
          call wfn_readeigenvec(oldwfnfile, frag1, fcomm1, vec1, ikeep, e, xj, xt2) ! KSM: updated
          call br_grab_vec1()! push vec1 -> br_reg
!.... MUST ORTHONORMALIZE HERE AND MUST KEEP!
! 	  call br_normalize(dnorm)

      !call br_normalize(scale)
          call br_add2hist(1)
	  end if

!write(logfile,*) "------------history matrix--------------"
!write(logfile,*) br_histbuf(:,:)
   enddo

   call wfn_close_file(oldwfnfile)

   !call exit
if(iproc==0)print*, " reading in pivot block done"
   return
end subroutine read_pivot_block_choose
!
!*----------------------------------------------------------------
!  Subroutine: initialize_pivot_block
!
!  Purpose:  Creates a pivot block of random orthonormal vectors within history.
!
!
!  Parameters:
!      dimblock(global IN)- size of block
!
!  Returns:
!*-------------------------------------------------------------------

subroutine initialize_pivot_block
use localblocks
use mod_reorthog
use mod_balg
implicit none
integer(kind=4) :: i
integer :: pivchoice
real(8) :: dnorm

pivchoice = 0

call br_set_histpos(0)
call br_reset_histbuf()
do i = 1,dimblock
	call add_pivot_br_reg(pivchoice,i)
	call br_normalize(dnorm)
	call br_orthogonalize(0)
	call br_normalize(dnorm)
    call br_add2hist(1)
enddo

end subroutine initialize_pivot_block
!====================================================================
!
! routine to initialize Lanczos vector
!
!  NOTE: currently very kludgy!
!
!  INPUT:
!     startiter:  where iterations start
!              needed for restart option
!
!  OUTPUT:
!     dnorm0  = if read in pivot, what is magnitue
!     vchar   = keeps track of which vectors is initial, which final
!             'n' (normal) initial vec1, final vec2
!             'r' (reverse) initial vec2, final vec1
!
!  CALLED BY:
!    setup_for_lanczos
!  SUBROUTINES CALLED:
!    readpivot (old version)  reads in pivot from file, normalizes
!    dnormvec_p  : normalizes a vector (or a fragment)
!
subroutine initialize_lanczos_vector
   use flagger
   use localvectors
!   use precisions
   use basis
   use fragments
   use io
!   use lanczos_info
   use onebodypot
   use nodeinfo
   use wfn_mod
   use bvectorlib_mod
   use menu_choices
   use localblocks
   use sporbit,only:allsameW
   use W_info
   use sectors
!   use mod_balg
   implicit none
   character(1) :: vchar
   real(kind=8) :: dnorm, dnorm0,da
   real :: rv
   integer startiter
   integer (kind = basis_prec) :: i
   integer (kind=4) :: iter
   logical smallflag
   integer :: pivotoption   ! different choices
   integer(8) :: countstates
!..... ADDED IN 7.8.7 to improve pivot....

   integer (kind=basis_prec) :: startvec,stopvec,npxsd,nnxsd
   integer :: isector,ncsector,icsector,csector
   integer :: ierr

   integer:: nseed,value(8)
   integer, allocatable:: seed(:)

   pivotoption = 0

   if(iproc==0)print*,' USING PIVOT OPTION ',pivotoption

   vchar = 'n'
   if ( strengthflag  .or. menu_char=='np' .or. greenflag) then

    if(piv_flag .eq. 'y' .or. blockstrengthflag) then

!        if(blockstrengthflag)then
			call read_pivot_block_choose   ! now default
!		else
!            call read_pivot_block   ! FIXED so it works with fragments
!		end if
    else
      call readpivot   ! FIXED so it works with fragments

	  return
   endif

   else if ( startiter == 0 ) then

      dnorm = dsqrt(real(dimbasis,kind=8))
1111 continue

      select case (pivotoption)

      case (0)
	  do i = v1s, v1e

!....... TO PREVENT ROUNDOFF ERRORS AND ALSO TO MINIMIZE
!        ACCIDENTAL STOPS, HAVE DECREASING WEIGHT
          vec1(i) = real(1.d0 /dnorm,kind=lanc_prec)
      end do
      case (1)
	  do i = v1s, v1e

!....... TO PREVENT ROUNDOFF ERRORS AND ALSO TO MINIMIZE
!        ACCIDENTAL STOPS, HAVE DECREASING WEIGHT
          call random_number(rv)
          vec1(i) = real( (1.d0+ 0.1*rv) /dnorm, kind=lanc_prec)
      end do

      case (2)
	  do i = v1s, v1e

!....... TO PREVENT ROUNDOFF ERRORS AND ALSO TO MINIMIZE
!        ACCIDENTAL STOPS, HAVE DECREASING WEIGHT

          vec1(i) = real(1.d0 /dsqrt(real(i,kind=8)+dnorm ), kind=lanc_prec)

       end do

!..... ALTERNATING SIGN SUGGESTED BY WEO..............
       case (3)

	  do i = v1s, v1e

          vec1(i) = real( 1.d0 * (-1)**i /dnorm, kind=lanc_prec)

       end do

!.................CHOOSING LOWEST W MAY MINIMIZE TRACE...

       case (4)

	   vec1 = 0.d0
	   if(allsameW)then
		   startvec=v1s
		   stopvec =v1e
	 	  do i = startvec,stopvec
	           vec1(i) = real( 1.d0 * (-1)**i /dnorm, kind=lanc_prec)
	       end do
	   else
		   countstates=0
		   do isector = 1,nsectors(1)
			   if(xsd(1)%sector(isector)%basisstart > v1e .or.  &
			     xsd(1)%sector(isector)%basisend < v1s  )cycle
			   if(xsd(1)%sector(isector)%Wx /= minW(1))cycle

			   do icsector = 1,xsd(1)%sector(isector)%ncsectors
				   csector = xsd(1)%sector(isector)%csector(icsector)
				   if(xsd(2)%sector(csector)%Wx /= minW(2))cycle
				   startvec = pstart( xsd(1)%sector(isector)%xsdstart) &
				   + nstart( xsd(2)%sector(csector)%xsdstart)
				   stopvec =startvec+ &
				    xsd(1)%sector(isector)%nxsd*xsd(2)%sector(csector)%nxsd-1

				   startvec = max(startvec,v1s)
				   stopvec  = min(stopvec,v1e)
				   if(stopvec < startvec)cycle

                   countstates=countstates+stopvec-startvec+1
!........ NOW LOOP OVER THESE SECTORS................
                   do i = startvec,stopvec
                      vec1(i) = real( 1.0 *(-1)**i /dnorm, kind=lanc_prec)
                   end do
			   end do
		   end do

		   call BMPI_ALLREDUCE(countstates,1,MPI_SUM,icomm,ierr)
		   if(countstates<10)then
			   pivotoption = 3
			   goto 1111
		   end if

	   end if

	   case (5,-5)

!	   if(nproc==1)then
	     call random_seed(size=nseed)
	     allocate(seed(nseed))
		 if(pivotoption==-5 .and. nproc==1)then
	        print*,' Enter random seed '
	        read*,seed(1)

		else
			call date_and_time(VALUES=value)  ! kluge to set random seed to time
			seed(1)=value(8)
			if(nseed>1)seed(2)=value(8)

		end if
		if(iproc==0)print*,' Creating pivot with random seeds ',seed

	     call random_seed(put=seed)
!	   end if
!	   print*,' seed ',nseed,seed
 	  do i = v1s, v1e

 !....... TO PREVENT ROUNDOFF ERRORS AND ALSO TO MINIMIZE
 !        ACCIDENTAL STOPS, HAVE DECREASING WEIGHT
           call random_number(rv)
!if(i< 11)print*,i,rv

           vec1(i) = real( (0.5d0-rv) /dnorm, kind=lanc_prec)
       end do

       end select

   end if

!----------- RESTART OPTION -------------
!          NEEDS WORK!

  if ( startiter > 0 ) then

     if(iproc==0)rewind(lvec_file)
     if(.not. useNewReorthog) call swap_vchar(vchar)
     do iter = 1,startiter
		  if(.not. useNewReorthog) call swap_vchar(vchar)
        call read_lanczos_vector_restart_a(vchar,'i',iter,lvec_file)
     end do
  endif

   return
end subroutine initialize_lanczos_vector
!====================================================================

!subroutine dynamic_vector_initializer
	! added 7.6.3, deleted 7.9.11
!====================================================================
!
! routine is mostly not used and should be phase out
! original idea was to go back and forth between vec1 and vec2
! where vchar='n' means normal:  H vec1 = vec2
!  and  vchar='r' means reverse: H vec2 = vec1
! but this idea doesn't lend itself well to MPI
!
subroutine swap_vchar(vchar)
  use fragments
  implicit none
  character(1) :: vchar

  if(useNewReorthog) then
     print *, "swap_vchar: should not call with new reorthog"
	 stop 1
  end if

  if(vchar == 'n')then
      vchar = 'r'
   else
      vchar = 'n'
   endif
   return
end subroutine swap_vchar

subroutine print_vec_head(vec, msg)
   use nodeinfo
   use precisions
   implicit none
   integer :: i
   real(kind=lanc_prec) :: vec(*)
   character (len=*) :: msg
   if(iproc /= 0) return

   ! print *, msg, ": ", vec(1), vec(2), vec(3), vec(4), vec(5)
   return
end subroutine print_vec_head

!================================================
!
!  WRITE OUT RESULTS
!
!  CALLS:
!   find_lanczos_eigenvalues
!   setup4obsmaster
!   reset_jumps4obs
!   read_lanczos_vector_a
!   block_reduce
!   applyobsbundled
!   clocker
!   wfn_writeeigenvec (formerly writeeigenvec_p)
!   close_lanczosfile
!   Wcounter
!   density1b_output
!   output_TRDENS
!   makespeme
!
subroutine lanczos_output(iter,i_thick,dnorm0)

  use sporbit
  use localvectors
  use nodeinfo
  use io
  use basis
  use obs
!  use lanczos_info
!  use precisions
  use fragments
  use mod_reorthog
  use flagger
  use coupledmatrixelements
  use system_parameters
  use wfn_mod
  use pocc_mod
  use butil_mod
  use btbme_mod
  use bvectorlib_mod
  use jump_mod
  use diagh
  use apply_obs_mod
  use mod_balg
  use localblocks
  implicit none

  integer(4) :: iter,i_thick
!  real    :: alpha(niter), beta(niter)
!  real(kind=egv_prec)    :: e(niter), eiglvec(niter,niter)
  real (kind =8 )::dnorm0,da    ! magnitude for strengths

  real(4) :: xj,xt
  logical :: flyeig
  integer(4) :: i,j
  integer(kind=basis_prec) :: k
  integer(kind=basis_prec) :: istate
!--------------------- CONVERGENCE -----------------------------

!  logical :: finished
  integer :: ncheck
  real, allocatable :: eold(:)  ! old energies
  real    :: ediff

  logical smallflag,zeroflag

  integer :: nread,jvec

  integer inode
  character(1) :: vchar
  integer :: ierr
  real(kind = 8) :: dtmp
  real(kind = 8) :: tmpvec(nkeep)
  ! real(kind=4),allocatable   :: spoccx(:,:,:)  ! single-particle occupations
                                   ! spoccx(it,istate,iorb)  species(p/n), eigenstate, orbit
  integer(4) :: iorb
  character  :: outchoice  !used for steering output options

  vchar = 'n'
  call find_lanczos_eigenvalues(iter,get_JT,i_thick)
  call clocker('tt8','end')

  call clocker('lan','end')
  call clocker('obs','start')

  call pocc_init_spoccx() ! allocate and initialize array

  call wfn_write_nkeep(nkeep) ! write number of vectors to wfn file
!---------- SET UP ARRAY FOR SINGLE-PARTICLE OCCUPATIONS --
  outchoice ='d'  ! default option
  if (spoccflag .and. .not.densityflag)then
     if(npeff(1)>0 .and. npeff(2) > 0)outchoice='b'
     if(npeff(1)>0 .and. npeff(2) == 0)outchoice='p'
     if(npeff(1)==0 .and. npeff(2) > 0)outchoice='n'
  end if
!------------ RECONSTRUCT EIGENVECTORS AND COMPUTE J^2, T^2.........

  if(get_JT)then

!------------ SET UP FOR ANG MOM................................

  call setup4obsmaster('J')
  call setup4obsmaster('T')
!.................. ADDED DEC 2010 by CWJ .........................
! IF USING 3-BODY (OR POSSIBLY SKIPPING ZEROES)
! NEED TO RESET ALL THE JUMPS

  call reset_jumps4obs
  end if
!...................................................................

  if ( iproc == 0 ) print*,' '
  if ( strengthflag ) then

	  if(.not.blockstrengthflag)then
          if ( iproc == 0 ) then

              write(resultfile,*)' Energy   Strength    '
              write(resultfile,*)' ______   ________  '

              write(6,*)' Energy   Strength  '
              write(6,*)' ______   ________  '

              do j = 1,iter
                  write(6,88)e(j),eiglvec(1,j)**2*dnorm0*dnorm0
                 if(writeout)write(resultfile,88)e(j),eiglvec(1,j)**2*dnorm0*dnorm0
              enddo
              write(resultfile,*)'_______   ________   '
			  write(resultfile,*)dnorm0**2,'= total strength, included above '
			  write(resultfile,*)' '
              write(6,*)'_______   ________  '
			  write(6,*)dnorm0**2,'= total strength, included above '
			  write(6,*)' '
              write(6,*)' '
           end if

	   else   ! BLOCK STRENGTH-- NEED TO CONVERT
!		   print*,' normalization ',orthonorm_pivot
	       do i = 1,dimblock
			   if(iproc==0)then

				   dtmp  = 0.d0

				   do j = 1,dimblock

					   dtmp = dtmp + orthonorm_pivot(j,i)**2
				   end do
				   write(6,*)'  '
				   write(6,*)'################### '
				   write(6,*)'  '
				   write(6,*)' For pivot state ',pivot_choice(i),' total strength = ',dtmp
	               write(6,*)' Energy   Strength   '
	               write(6,*)' ______   ________  '
				   write(resultfile,*)'###################'
				   write(resultfile,*)' For pivot state ',pivot_choice(i),' total strength = ',dtmp
	               write(resultfile,*)' Energy   Strength  '
	               write(resultfile,*)' ______   ________  '

				   do j = 1,iter
				      dtmp = 0.d0
				      do k = 1,dimblock
					      dtmp = dtmp +   orthonorm_pivot(k,i)*eiglvec(k,j)
				       end do

				       write(6,88)e(j),dtmp*dtmp !, dtmp*dtmp
			           if(writeout)write(resultfile,88)e(j),dtmp*dtmp !,dtmp*dtmp

		            end do ! j
				end if  ! iproc == 0
		   end do

       endif
!     return
  endif
88 format(3f10.5)

!------------------ APPLY RESOLVANT/GREEN FUNCTION AND WRITE OUT----------------------
!   (Added 7.9.10
if(greenflag)then

	eiglvec = 0.
	do j = 1,iter
		eiglvec(j,1)=gvec(j)*dnorm0

	end do
    call br_transform_basis(eiglvec, nkeep, iter)
    call br_retrieve_hist(1)
    ! It turns out that we need the vector loaded into both vec1 and vec2
    ! of course, they are different slices on each node
    call br_restore_vec1()
    call wfn_writeeigenvec(wfnfile, frag1, vec1, 1, 0.0, 0.0, 0.0)
    if(iproc==0)then
  	  print*,' '
  	  print*,' Note: output wfns not normalized  '
    end if

	call wfn_close_file(wfnfile)

	if(dotflag)then
		if(iproc==0)print*,' Finally, read in a previously computed file to take dot product '
		call overlap(.true.)
	end if

	return

end if
!--------------------------------------------------------------

!---------------- WRITE OUT SINGLE PARTICLE STATES with occupation report --------
  if(spoccflag .and. .not.densityflag)  call pocc_write_orbits()

  call pocc_write_table_header() !  E   Ex ...

  allocate(energy(nkeep), xjlist(nkeep),xtlist(nkeep), stat=aerr )
  if(aerr /= 0) call memerror("lanczos_output 1");

!----- Construct nkeep eigenvectors of H from  Lanczos eigenvectors....
!---- If lanczos vectors stored in core, then construct in core

if(get_JT)then                                           !added 7.6.8
  if(storelanczosincoreMPI .or. storelanczosincore1)then
     call clocker('eig','sta')
     ! KSM lvec(bitpos, iter) is the old way of storing old vectors
     ! It will go away with useNewReorthog
     if(useNewReorthog) then
        ! This applys the same transform as the else clause, but to the
        ! basis stored in br_histbuf

        call br_transform_basis(eiglvec, nkeep, iter)

     else

        ! KSM 17Aug2014 - my interpretation of the following code
        ! We have a sub basis stored in lvec(:, i)
        ! We extracted the \alpha and \beta coefficients against this basis
        ! This represents a projection of H onto this much smaller basis.
        ! The eigenvectors of the \alpha,\beta matrix can be used to construct
        ! approximate eigenvectors in the untruncated space by taking the linear
        ! combination of lvec(:, i) with weights being the component of the eigenvector.
        !
!$omp parallel do private(tmpvec,i,j,dtmp), shared(eiglvec, Lvec)
        do k = 1,Ldim     ! each row (labeled by kl) can be transformed independently
           dtmp = 0.d0
           do i = 1,nkeep
               dtmp = 0.d0
               do j = 1,iter
                   dtmp = dtmp +   real(lvec(k,j),kind=8)* real(eiglvec(j,i),kind=8)
               end do ! j
               tmpvec(i) = dtmp
           end do  ! i
           do i = 1,nkeep
               lvec(k,i) =real(tmpvec(i),kind=lanc_prec)
           end do
        end do ! k
!$omp end parallel do
     end if
     call clocker('eig','end')
     vchar = 'n'
!---- COMPUTE J,T AND WRITE OUT VECTORS-------------------------------
     do i = 1,nkeep
        ! KSM - the idea here is to restore a history vector into vec1
        if(useNewReorthog) then
           call br_retrieve_hist(i)
           ! It turns out that we need the vector loaded into both vec1 and vec2
           ! of course, they are different slices on each node
           call br_restore_vec1()
        else
           vec1 = 0.0
           call read_lanczos_vector_a(vchar,'i',i,lvec_file)
           if(storelanczosincoreMPI) then
              ! call block_reduce(dimbasis,vec1)
              ! each mpi process reads one slice.  allreduce is overkill but works
              call BMPI_ALLREDUCE(vec1, size(vec1), MPI_SUM, icomm, ierr) ! in place
           end if
        end if
        call br_load_vec2_from_vec1()

        energy(i) = real(e(i), kind(energy(1)))
!------------ COMPUTE J2, T2.................................
        call clocker('aob','sta')
        xj2 = 0.0e0_obs_prec
        xt2 = 0.0e0_obs_prec
        twoobsflag = .true.  ! compute both J^2 and T^2
        call applyobsbundled(1)
        call clocker('aob','end')
        xj = real( -0.5 + sqrt(xj2 + 0.25), kind=4)
        xt = real( -0.5 + sqrt(xt2 + 0.25), kind=4)

        xjlist(i) = xj
        xtlist(i) = xt
!!----------------WRITE OUT WFN..............................................
!!  KSM: do this before computing occupations because data like xt2 will be overwritten
!!       not clear why we are saving xt2 instead of xt, but I'm not changing it.
        call clocker('wev','sta')

        if ( writeout .and. write_wfn) then

   		   if(.not.strengthnorm .and. strengthflag .and. .not. blockstrengthflag)then ! added 7.9.7
 !  			  print*,eiglvec(1,i)*dnorm0,' testing '
   			  vec1=vec1*abs(eiglvec(1,i)*dnorm0)
   		   end if

           call wfn_writeeigenvec(wfnfile, frag1, vec1, i, real(e(i),kind=4), xj, real(xt2,kind=4))
        end if
        call clocker('wev','end')

!----------- COMPUTE SINGLE PARTICLE OCCUPATIONS---------- added in 7.3.7
        call pocc_compute_spocc(i, .true.)  ! true for restore J/T setup

!----------------- WRITE OUT RESULTS.............................
           if ( iproc == 0 ) then
           select case (outchoice)

               case('d')
                 call pocc_write_ejt(i, e, xj, xt)

               case('b')
                 ! write(6,12)i,e(i), e(i) - e(1),xj,xt,(spoccx(1,i,iorb),iorb=1,numorb(1))
                 call pocc_write_ejt(i, e, xj, xt)
                 call pocc_write_occvec(6, spoccx, i, 1, "    p occ:")
                 call pocc_write_occvec(6, spoccx, i, 2, "    n occ:")
                 if ( writeout ) then
                    call pocc_write_occvec(resultfile, spoccx, i, 1, "    p occ:")
                    call pocc_write_occvec(resultfile, spoccx, i, 2, "    n_occ:")
                 end if
               case('p')
                 ! write(6,12)i,e(i), e(i) - e(1),xj,xt,(spoccx(1,i,iorb),iorb=1,numorb(1))
                 call pocc_write_ejt(i, e, xj, xt)
                 call pocc_write_occvec(6, spoccx, i, 1, "    p occ:")

                 if ( writeout ) then
                    !  write(resultfile,12)i,e(i),e(i) - e(1),xj,xt,(spoccx(1,i,iorb),iorb=1,numorb(1))
                    call pocc_write_occvec(resultfile, spoccx, i, 1, "    p occ:")
                 end if

               case('n')
                 ! write(6,13)i,e(i), e(i) - e(1),xj,xt,(spoccx(2,i,iorb),iorb=1,numorb(2))
                 call pocc_write_ejt(i, e, xj, xt)
                 call pocc_write_occvec(6, spoccx, i, 2, "    n occ:")

                 if ( writeout ) then
                    ! write(resultfile,13)i,e(i),e(i) - e(1),xj,xt,(spoccx(1,i,iorb),iorb=1,numorb(2))
                    call pocc_write_ejt(i, e, xj, xt)
                    call pocc_write_occvec(resultfile, spoccx, i, 2, "    n occ:")
                 end if
            end select
        end if
     end do  ! i

!.....  else construct by writing to disk  (slow).....
  else
   if(useNewReorthog) then
      if(iproc == 0) print *, "useNewReorthog  oops"
      stop 1
   end if
   do i = 1, nkeep

     do k = 1, dimbasis
        vec1(k) = 0.0e0_lanc_prec
     end do !k

     call clocker('eig','sta')
     do jvec = 1,iter
        da = real(eiglvec(jvec,i),kind=8)
        if(.not.storelanczosincoreMPI)then
           call read_lanczos_vector_a(vchar,'f',jvec,lvec_file)
           do k = 1, dimbasis
              vec1(k) = vec1(k) + real(vec2(k),kind=8)*da
           end do  !k
           ! need data on both sides:  note that nodes don't have the same slices on both sides
           call br_load_vec2_from_vec1()
        else
!$omp parallel do private(istate, k)                               &
!$omp              shared(Lstart, Ldim, vec1, Lvec, jvec, da)
           do istate = 1, Ldim
              k = Lstart + istate - 1
              vec1(k) = vec1(k) + real(Lvec(istate,jvec),kind=8)*da
           end do  !istate
!$omp end parallel do
        end if
     end do  !j

     if(storelanczosincoreMPI) then
        ! call block_reduce(dimbasis,vec1)
        ! each mpi process reads one slice.  allreduce is overkill but works
        call BMPI_ALLREDUCE(vec1, size(vec1), MPI_SUM, icomm, ierr) ! in place
     end if
     call clocker('eig','end')

!------------ COMPUTE J2, T2.................................
     call clocker('aob','sta')
     xj2 = 0.0e0_obs_prec
     xt2 = 0.0e0_obs_prec
     twoobsflag = .true.  ! compute both J^2 and T^2
     call applyobsbundled(1)
     call clocker('aob','end')
     xj = real( -0.5 + sqrt(xj2 + 0.25), kind=4)
     xt = real( -0.5 + sqrt(xt2 + 0.25), kind=4)

!----------- COMPUTE SINGLE PARTICLE OCCUPATIONS---------- added in 7.3.7
     if(spoccflag .and. .not.densityflag)then
        do iorb = 1,bmax(numorb(1),numorb(2) )
           if(numorb(1) > 0 .and. iorb < numorb(1))then
              pspe = 0.0
              pspe(iorb) = 1.0
              call makespeme(1,'H')
        end if
           if(numorb(2) > 0 .and. iorb < numorb(2))then
              nspe = 0.0
              nspe(iorb) = 1.0
              call makespeme(2,'H')
           end if
           xj2 = 0.0e0_obs_prec
           xt2 = 0.0e0_obs_prec
           call applyspoccbundled(1)
           if(np(1) > 0)spoccx(1,i,iorb)= xj2
           if(np(2) > 0)spoccx(2,i,iorb)= xt2
        end do
     end if
!----------------- WRITE OUT RESULTS.............................
!     if ( isoflag ) then
     if ( iproc == 0 ) then
        select case (outchoice)

            case('d')
               write(6,11)i,e(i), e(i) - e(1),xj,xt
               if ( writeout ) write(resultfile,11)i,e(i),e(i) - e(1),xj,xt

            case('b')
               write(6,121)i,e(i), e(i) - e(1),xj,xt
               call pocc_write_occvec(6, spoccx, i, 1, "    p occ:")
               call pocc_write_occvec(6, spoccx, i, 2, "    n occ:")
               if ( writeout ) then
                  write(resultfile,121)i,e(i), e(i) - e(1),xj,xt
                  call pocc_write_occvec(resultfile, spoccx, i, 1, "    p occ:")
                  call pocc_write_occvec(resultfile, spoccx, i, 2, "    n occ:")
               end if
            case('p')

            case('n')
         end select
      end if
      call pocc_write_orbits()

!----------------WRITE OUT WFN..............................................

     call clocker('wev','sta')
     if ( writeout .and. write_wfn) then

		 if(.not.strengthnorm .and. strengthflag .and. .not. blockstrengthflag)then ! added 7.9.7
			 print*,eiglvec(1,i)*dnorm0,' testing '
			 vec1=vec1*abs(eiglvec(1,i)*dnorm0)
		 end if

        call wfn_writeeigenvec(wfnfile, frag1, vec1, i, real(e(i),kind=4), xj, real(xt2,kind=4))
     end if
     call clocker('wev','end')
  end do
  end if

  else    ! just print out energies for get_JT=.false.
	  if(iproc==0)then
          do i = 1,nkeep
             write(6,11)i,e(i), e(i) - e(1)
             if ( writeout ) write(resultfile,11)i,e(i),e(i) - e(1)
	      end do
      end if

  end if
  ! add blank line for formatting
  if(iproc == 0) then
     write(6,*) " "
     if(writeout) write(resultfile, *) " "
  end if

11 format(i5,3x,2f10.5,2x,2f8.3)
12 format(i5,3x,2f10.5,2x,2f8.3,' p occ: ',20f7.3)
121 format(i5,3x,2f10.5,2x,2f8.3)
13 format(i5,3x,2f10.5,2x,2f8.3,' n occ: ',20f7.3)
14 format(46x,' n occ: ',12f7.3)

if(densityflag) then
	if(iproc==0)write(denresultfile,*)' '
	 call pocc_write_orbits_alt(denresultfile)
 end if

  if(.not.allsamew .and. get_JT)call Wcounter    ! added in 7.2.5; computes and prints out truncation components

  call close_lanczosfile

!........... ADDED in 7.9.7: OPTION TO CARRY FORWARD NORMALIZATIONS....

  if(.not.strengthnorm .and. strengthflag .and. .not. blockstrengthflag .and. iproc==0)then
	  print*,' '
	  print*,' Note: output wfns not normalized;  '
	  print*,' they carry forward input normalization and strength'
	  print*,'  (To change this, use menu option (sn))'
	  print*,' '
	  write(logfile,*)' '
	  write(logfile,*)' Note: output wfns not normalized;  '
	  write(logfile,*)' they carry forward input normalization and strength'
	  write(logfile,*)'  (To change this, use menu option (sn))'
	  write(logfile,*)' '
	  write(resultfile,*)' '
	  write(resultfile,*)' Note: output wfns not normalized;  '
	  write(resultfile,*)' they carry forward input normalization and strength'
	  write(resultfile,*)'  (To change this, use menu option (sn))'
	  write(resultfile,*)' '

  end if
  if(strengthnorm .and. strengthflag .and. .not. blockstrengthflag .and. iproc==0)then
	  print*,' '
	  print*,' Note: output wfns are normalized;  '
	  print*,'  (To change this, use menu option (su))'
	  print*,' '
	  write(logfile,*)' '
	  write(logfile,*)' Note: output wfns are normalized;  '
	  write(logfile,*)'  (To change this, use menu option (su))'
	  write(logfile,*)' '
	  write(resultfile,*)' '
	  write(resultfile,*)' Note: output wfns are normalized;  '
	  write(resultfile,*)'  (To change this, use menu option (su))'
	  write(resultfile,*)' '

  end if

  if(densityflag)call density1b_output  !(e,eiglvec)

  if ( trdensout  ) then
     call clocker('trd','sta')
     call output_TRDENS
     call clocker('trd','end')
  end if

  call pocc_cleanup()

  return
end subroutine lanczos_output
!================================================
!
!  function to force conversion of unconverged xJ to integer J
!  that is, odd 2 x J for odd A, and even for even A
!
  function closest2J(evenA,xj)

  implicit none
  integer closest2J
  real xj
  logical evenA

  if(evenA)then
     closest2J = 2*nint(xj)
     if(closest2J < 0)closest2J = 0
  else
     closest2J = 2*nint(xj-0.5)+1
     if(closest2J < 1)closest2J = 1
  end if

  return
  end function closest2J

!---- NOTE: density1b_output moved to bdenslib1.f90
!=============================================================================
!
!  thick-restart lanczos:
!  after some iterations, take nthick lowest eigenvectors
!  and reintroduce them as lanczos vectors
!  then start up iterations again
!
!  NB: in normal runs, up to this point have done iter = nthick+1 iterations
!
! nthick = # of vectors to keep for thick-restart
! iter = # of lanczos iterations in total
! iter_thick = which iteration of thick_restart (starting at 0)
! maxiter = dimension of arrays
! e(:), eiglev(:,:) : eigenvectors
! alpha(),beta()
! vchar (no longer in significant use)
!
! ALSO USED but not passed in arguments: niter, dimension of truncated Hamiltonian to diagonalize
! NOT USED here: nkeep, the # of final eigenvectors
! CALLED BY:
!   lanczos_p
! CALLS:
!  find_lanczos_eigenvalues
!  swap_vchar
!  read_lanczos_vector_a
!  write_lanczos_vector_a
!  block_reduce
!
subroutine thick_restart_sub_p(nthick,iter,iter_thick,maxiter,vchar)
!	subroutine thick_restart_sub_p(nthick,iter,iter_thick,maxiter,e,eiglvec,alpha,beta,vchar)

  use localvectors
  use precisions
  use basis
  use lanczos_info
  use nodeinfo
  use fragments
  use bmpi_mod
  use butil_mod
  use bvectorlib_mod
  use mod_reorthog
  use flagger
  use mod_balg,only:lanczos_mat
  implicit none

  integer  :: nthick  ! # of vectors to keep for thick-restart
  integer  :: iter    ! # of lanczos iterations so far
  integer  :: iter_thick   ! which thick_restart iteration this is, starting with 0
  integer :: maxiter     ! dimension
!  real(kind=egv_prec) :: e(maxiter), eiglvec(maxiter,maxiter)
  real(kind=8)   :: tmpvec(maxiter)
!  real(kind=4) :: alpha(maxiter), beta(maxiter)
  character(1) :: vchar

  real (kind = lanc_prec),pointer :: v1(:), v2(:)

  integer nthicktmp
  integer(4):: i,j,k,l
  integer(4):: mymax  ! for targeted thick-restart
  integer(4):: iclosest
  real(kind=8) :: eclosest
  logical :: verbosetarget = .true.  ! used for debugging excited state thick-restart
  integer, parameter :: tmpfile = 58
  real (kind = lanc_prec) :: da
  integer(kind=basis_prec) :: kl
  integer(4) :: file_mode,file_info,ierror
  integer(4) :: inode
  real(8) :: lanczos_mat_tmp(niter,niter)

  real(kind = 8) :: dtmp

  integer :: ierr
  real(kind=8)  :: thickconverge

  if(iproc==0)print*,' In thick restart niter, nthick = ',niter,nthick

!-------------------SOLVE FOR TRUNCATED EIGENVALUES -----------

  if(iter_thick == 0)then
    nthicktmp = 0
  else
    nthicktmp = nthick
  endif
  call find_lanczos_eigenvalues(iter-1,.true.,nthicktmp)

!--------------------CREATE NEW COEFICIENTS -------------------
  do i = 1,nthick
     alpha(i) = e(i)
     beta(i) = beta(iter-1)*eiglvec(iter-1,i) ! correction 6/2011 by WEO
  enddo

!----- ALTERNATE, MORE EXACT but more time-consuming METHOD ADDED IN 7.9.0
!      recompute the matrix elements by applying the Hamiltonian

if(lanczchar=='tx' )then

	if(txchoice==1)then
		lanczos_mat = 0.d0

	   do i = 1,nthick
		  lanczos_mat(i,i)=e(i)
		  lanczos_mat(i,nthick+1)=beta(iter-1)*eiglvec(iter-1,i)
		  lanczos_mat(nthick+1,i)=beta(iter-1)*eiglvec(iter-1,i)
	   end do

   else
	   print*,' That choice is not implemented'
	   stop 909

    end if
end if

!------------ CHECK SIZE OF COUPLING -------------------------

  thickconverge = 0.d0
  if(lanczchar=='tx' .and. txchoice==2)then
      do i = 1,nkeep
  	     thickconverge = thickconverge + lanczos_mat(i,nthick+1)**2/real(nkeep,kind=8)
      end do

  else
      do i = 1,nkeep
  	     thickconverge = thickconverge + beta(i)**2/real(nkeep,kind=8)
      end do
  end if
  if(iproc==0)print*,' Alternate thick-restart convergence : ',thickconverge

!.................. REWRITE TO FILE......
  if(iproc ==0)then
     rewind(coef_file)
     print*,' rewriting to file ',coef_file,nthick
     do i = 1,nthick
       write(coef_file,'(i10,2(2x,f20.15))')-i,alpha(i), beta(i)
     end do
  end if

  if(useNewReorthog) then
     ! use eigenvectors of tridiagonal alpha-beta matrix to
     ! mix saved lanczos vectors into approximate eigenstates of H.
     ! we will use these as the first nthick saved lanczos vectors (the history)
     ! going forward.
     ! I understand alpha(i) = e(i) above
     ! I think the beta issue is related to the last vector produced in
     ! the lanczos process.   We have dropped the last beta and associated vector
     ! by truncating the alpha-beta matrix at nthick.
     if(iproc .eq. 0) then
        print *, "New Reorthog Thick Restart"
        print *, "iter=", iter, ", nthick=", nthick
        print *, "br_histpos=", br_histpos
        print *, "nthick_add=", nthick_add
     end if

	 if(lanczchar=='tx' .and. txchoice==2)then
         call br_transform_basis(eiglvectx, nthick, iter-1)  ! ? iter-1 or iter-2?

	 else

        call br_transform_basis(eiglvec, nthick, iter-1)

	end if

     ! init for next iteration.
     ! load vec1 from the top of the history.
     call br_retrieve_hist(iter)
     call br_restore_vec1()

     ! ends with nthick as the most recent vector.
     call br_set_histpos(nthick)  ! set last written position
     call br_add2hist(iter+1) ! br_retrieve_hist left in br_reg, put at new head
     return
  end if

!................IF STORING LANCZOS IN CORE .....
  if(storelanczosincore1 .or. storelanczosincoreMPI)then

!$omp parallel do private(tmpvec,i,j,dtmp), shared(eiglvec, Lvec,nthick,iter,Ldim)
     do kl = 1,Ldim     ! each row (labeled by kl) can be transformed independently
        dtmp = 0.d0
        do i = 1,nthick
            dtmp = 0.d0
            do j = 1,iter-1
                dtmp = dtmp +   real(lvec(kl,j),kind=8)* real(eiglvec(j,i),kind=8)
            end do ! j
            tmpvec(i) = dtmp
        end do  ! i
        do i = 1,nthick
            lvec(kl,i) =real(tmpvec(i),kind=lanc_prec)
        end do
     end do ! kl
!$omp end parallel do

!........ SET VCHAR...
     vchar = 'n'
     do i = 1,nthick
        call swap_vchar(vchar)
     end do
     if(vchar=='n')then
        do kl = 1,dimbasis
           vec1(kl) = 0.0
        end do
     end if
     if(vchar=='r')then
        do kl = 1,dimbasis
           vec2(kl) = 0.0
        end do
     end if
     if(useNewReorthog .and. iproc == 0) print *, "read_lanczos_vector in thick_restart_sub_p"
     call read_lanczos_vector_a(vchar,'i',iter,tmpfile)
     if(storelanczosincoreMPI) then
        if(vchar == 'n') call BMPI_ALLREDUCE(vec1, size(vec1), MPI_SUM, icomm, ierr) ! in place
        if(vchar == 'r') call BMPI_ALLREDUCE(vec2, size(vec2), MPI_SUM, icomm, ierr) ! in place
     end if
     return
  end if
!................ IF STORING LANCZOS VECTORS ON DISK.................
  vchar = 'n'
  open(unit=tmpfile,status ='scratch',form='unformatted')

  do i = 1,nthick
     if(vchar == 'n')then
        v1 => vec1
        v2 => vec2
     else
        v1 => vec2
        v2 => vec1
     end if
     do kl = 1,dimbasis
        v1(kl) = 0.0
     enddo ! j
     do j = 1,iter-1
         if(i == 1)then
              call read_lanczos_vector_a(vchar,'f',j,lvec_file)
              call write_lanczos_vector_a(vchar,'f',j,tmpfile)
         else
              call read_lanczos_vector_a(vchar,'f',j,tmpfile)
         endif
         da = real(eiglvec(j,i), kind(da))
         do kl = 1,dimbasis
            v1(kl) = v1(kl) + v2(kl)*da
         enddo  !kl
     enddo  !j
     if(i == 1)then
         call read_lanczos_vector_a(vchar,'f',iter,lvec_file)
         call write_lanczos_vector_a(vchar,'f',iter,tmpfile)
     endif
     call write_lanczos_vector_a(vchar,'i',i,lvec_file)
     call swap_vchar(vchar)

  enddo  ! i

  call read_lanczos_vector_a(vchar,'i',iter,tmpfile)
333 format(i4,4f10.6)
  if(nproc ==1)then
     close(tmpfile)
  else
     call BMPI_FILE_CLOSE(tmpfile,ierror)
  end if

  return
end subroutine thick_restart_sub_p

!================================================
!
!  caution: not thoroughly tested
!
!  CALLS:
!  dnormvec_p
!  reorthogonalize_a
!
subroutine random_restart_p(iter,vchar)

   use localvectors
   use fragments
!   use precisions
   use basis
!   use lanczos_info
   use nodeinfo
   use bmpi_mod
   use butil_mod
   use bvectorlib_mod
   implicit none

   character(1) :: vchar
   integer(4) :: iter
   real(kind=lanc_prec), pointer :: w(:)
   real(kind=8):: da
   integer i,j
   real(kind=8) :: dnorm
   logical :: smallflag
   real :: rv
   integer(kind=basis_prec) :: jl,vstart,vstop

!   integer iseed =1959   ! to be used for restarting lanczos if it stops

   if(vchar == 'n')then
      w => vec1
      vstart = basestart(frag1)
      vstop  = basestop (frag1)
   else
      w => vec2
      vstart = basestart(frag2)
      vstop  = basestop (frag2)
   endif

   dnorm = 1.d0 /dsqrt(real(dimbasis,kind=8))

   do jl = vstart,vstop
      call random_number(rv)
      w(jl) = dnorm*(rv-0.5)    !ran(iseed)       ! RANDOM FUNCTION -- may need another
   end do
   call dnormvec_p(vchar,'i',dnorm,smallflag)

!----------- orthogonalize against previous vectors

   call reorthogonalize_a(iter,vchar,dnorm)
   call dnormvec_p(vchar,'i',dnorm,smallflag)

   if(smallflag)then
       if(iproc == 0)print*,' Ooops, zero vec in reorthogonalization of restart '
       stop
   endif

   return
end subroutine random_restart_p

!===========================================================================
! Open Lanczos file
! revision 6/2012 by CWJ
! requires further modification when vectors stored in memory
!
! CALLED BY:
!	lanczos_master in BLANCZOS.f90
!
!  NOTE: in 7.8.1, opening .lvec file turned OFF
!
!===========================================================================
subroutine open_lanczosfile
  use nodeinfo
  use lanczos_info
  use precisions
  use io
  use bmpi_mod
  use butil_mod
  implicit none
  character (len=25) :: filename
  character (len=4)  :: proc_name
  integer(4)         :: ilast
  integer(4)         :: file_mode
  integer(4)         :: file_info
  integer(4)         :: ierror
  writetodisk = .true.

  if ( writeout .and. iproc == 0) then
     ilast = index(outfile,' ') - 1
!     open(unit=lvec_file,file=outfile(1:ilast)//'.lvec',status = 'unknown',form='unformatted')
     open(unit=coef_file,file=outfile(1:ilast)//'.lcoef',status = 'unknown',form='formatted')

  elseif(iproc==0)then
     filename = 'lanczosvec'
     ilast = index(filename,' ')-1
!     open(unit = lvec_file,file=filename(1:ilast)//'.lvec',status='unknown',form ='unformatted')
     open(unit = coef_file,file=filename(1:ilast)//'.lcoef',status='unknown',form ='formatted')

  end if
  return
end subroutine open_lanczosfile

!===========================================================================
!  subroutine close_lanczosfile
! revision 6/2012 by CWJ
!===========================================================================
subroutine close_lanczosfile
  use lanczos_info
  use io
  use nodeinfo
  use bmpi_mod
  use butil_mod
  implicit none

  integer(4) :: ierror
  if ( writetodisk ) then
     if ( .not. restart_enabled ) then
        if ( iproc == 0 ) then ! use normal I/O
           rewind(lvec_file)
           write(lvec_file)0
        end if
     end if
     if ( iproc == 0 ) then
        close(lvec_file)
     end if
     if (iproc == 0) close(coef_file)
  end if
  return
end subroutine close_lanczosfile

!======================================================================
!
!  routine to compute distribution of W values in final wavefunctions
!  added in 7.2.5 by CWJ SDSU  11/2013
!
subroutine Wcounter
   use W_info
   use sporbit
   use basis
   use sectors
   use io
   use nodeinfo
   use lanczos_info
   use localvectors
   use flagger
   use fragments
!   use  tribution
   use mod_reorthog
   use wfn_mod
   use bmpi_mod
   use butil_mod
   use bvectorlib_mod
   implicit none

   integer(4) :: ierr
   integer is,isc,jsc
   integer(8) :: ip,in
   integer(kind=basis_prec) :: ibasis
   integer Wp, Wn,W
   logical, allocatable :: Wallowed(:)
   integer :: nWvals,n
   integer, allocatable :: Wlist(:)
   real(8), allocatable :: Wfrac(:)
   real(8) ::  ftmp,dv
   integer istate
   integer :: idummy
   real(4) :: xj,xt,ei,ef,xtt,xjj
   integer :: aerr

   if(allsameW)return

   allocate(Wallowed(minWtot:maxWtot), stat=aerr )
   if(aerr /= 0) then
      call memerror("Wcounter 1")
      stop 5
   end if
   Wallowed(:) = .false.
!................ FIND POSSIBLE W IN BASIS.................
   do is = 1,nsectors(1)
      wp = xsd(1)%sector(is)%Wx

      do isc = 1,xsd(1)%sector(is)%ncsectors
         jsc = xsd(1)%sector(is)%csector(isc)
         wn  = xsd(2)%sector(jsc)%Wx
!................ ERROR TRAP..............
         if( Wp + Wn < minWtot .or. Wp+Wn > maxWtot)then
             if(iproc==0)print*,' wrong ws ',wp,wn, minwtot,maxwtot
             stop
         end if
         Wallowed(wp+wn)=.true.
      end do  ! isc
   end do   ! is
   nWvals = 0
!.............. COUNT HOW MANY................................
   do w = minWtot,maxWtot
      if(Wallowed(w))nWvals = nWvals+1
   end do
   if(nWvals==1)then
       if(iproc==0)then
           print*,' '
           print*,' Only one value of W in truncation encountered '
           print*,' '
           if(writeout)then
                write(resultfile,*)' '
                write(resultfile,*)' Only one value of W in truncation encountered '
                write(resultfile,*)' '
           end if
       end if
       return
   end if
   allocate(wlist(nWvals), stat=aerr)
   if(aerr /= 0) call memerror("Wcounter 2")
   n = 0
!........... CREATE LIST OF Ws............
   do w = minWtot,maxWtot
      if(Wallowed(w))then
            n=n+1
            Wlist(n) = W
      end if
   end do
   allocate(Wfrac(minWtot:maxWtot), stat=aerr )
   if(aerr /= 0) call memerror("Wcounter 3")
!................ GO THROUGH WAVEFUNCTIONS AND COUNT UP W.....
   if(iproc==0)then
       write(6,*)' '
       write(6,*)'       % occupation of W-truncation subspaces '
       write(6,'(a8,10i7)')'State/W=',(wlist(n)-minwtot,n=1,nWvals)
       if(writeout)then
           write(resultfile,*)' '
           write(resultfile,*)'       % occupation of W-truncation subspaces '
           write(resultfile,'(a8,10i7)')'State/W=',(wlist(n)-minWtot,n=1,nWvals)
       end if
   end if

!...................LOOP OVER WFNS..............................

  if(.not.storelanczosincore1 .and. .not.storelanczosincoreMPI)then
     call wfn_rewind(wfnfile)
     call read_wfn_header(wfnfile,.false.)
     call wfn_read_nkeep(wfnfile, n)  ! dummy reading in nkeep
     ! read(wfnfile)n   ! dummy reading in nkeep
  endif
  do istate = 1,nkeep
     if(storelanczosincore1 .or. storelanczosincoreMPI)then
        if(useNewReorthog) then
           call br_retrieve_hist(istate)
           ! It turns out that we need the vector loaded into both vec1 and vec2
           ! of course, they are different slices on each node
           call br_restore_vec1()
        else
           vec1 = 0.0 ! all ranks
           call read_lanczos_vector_a('n','i',istate,lvec_file) ! read in rank=0
           if(storelanczosincoreMPI) then
              ! call block_reduce(dimbasis,vec1)
              ! each mpi process reads one slice.  allreduce is overkill but works
              call BMPI_ALLREDUCE(vec1, size(vec1), MPI_SUM, icomm, ierr) ! in place
           end if
        end if
     else
        ! new interface, we say which vector to read and it checks
        call wfn_readeigenvec(wfnfile,frag1, fcomm1, vec1,istate,ei,xj,xtt)
     end if

     wfrac(:) = 0.d0

     ! figure out if this node is the first node in ifragment.
     if(nodal(iproc)%ifirst) then
        do is = 1,nsectors(1)
           wp = xsd(1)%sector(is)%Wx

           do isc = 1,xsd(1)%sector(is)%ncsectors
              jsc = xsd(1)%sector(is)%csector(isc)
              wn  = xsd(2)%sector(jsc)%Wx
              ftmp = 0.d0
              do ip = xsd(1)%sector(is)%xsdstart,xsd(1)%sector(is)%xsdend
                 do in = xsd(2)%sector(jsc)%xsdstart,xsd(2)%sector(jsc)%xsdend
                    ibasis = pstart(ip)+nstart(in)
                    if(ibasis .ge. v1s .and. ibasis .le. v1e) then
                       dv = vec1(ibasis)
                       ftmp = ftmp + dv*dv
                    end if
                 end do
              end do
              wfrac(wp+wn) = wfrac(wp+wn)+ftmp
           end do ! isc
        end do   ! is
     end if
     if(useNewReorthog) then
        ! have to reduce.  Note that condition %ifirst above suppresses nodes that
        ! are not "first" in their fragment.   This is not very efficient, but it works.
        call BMPI_ALLREDUCE(wfrac(:), SIZE(wfrac), MPI_SUM, icomm, ierr) ! in place reduce
     end if
     if(iproc==0)then
        write(6,'(i4,4x,10f7.2)')istate,(Wfrac(wlist(n))*100.,n=1,min(nWvals,10))
        if(writeout)write(resultfile,'(i4,4x,10f7.2)')istate,(Wfrac(wlist(n))*100.,n=1,min(10,nWvals))
		if(nWvals> 10)then
	        write(6,'(i4,4x,10f7.2)')istate,(Wfrac(wlist(n))*100.,n=11,nWvals)
	        if(writeout)write(resultfile,'(i4,4x,10f7.2)')istate,(Wfrac(wlist(n))*100.,n=11,nWvals)

		end if
     end if
  end do  ! istate

  return
end subroutine Wcounter

!======================================================================
!  added 7.9.10; computes 1/(Egf - H ) on pivot
!
!
subroutine resolvent(iter,ddot)

	implicit none
	integer, intent(in) :: iter
	real(kind=8) :: htmp(niter+2,niter+2),b(niter+2,1)
	integer :: i
	integer :: ipiv(niter+1)
	integer :: info
	real(8), intent(out) :: ddot
	real(8) :: norm1,norm2

	htmp = 0.d0

!........ SET UP MATRIX TO BE INVERTED....
	do i = 1,iter
		htmp(i,i)=Egf-alpha(i)
		if(i==iter)cycle
		htmp(i+1,i)=-beta(i)
		htmp(i,i+1)=-beta(i)
	end do
	b = 0.d0
	b(1,1)=1.d0
	call dgesv(iter,1,htmp,niter+2,ipiv,b,niter+2,info)

	if(info/=0)then
		print*,' error in inverting E- H in subroutine resolvent '
		stop 919
	end if

	do i = 1,iter
		gvec(i)=b(i,1)
	end do
	ddot = 0.d0
	if(iter ==1)return
!........... CHECK FOR CONVERGENCE.............
	norm1= 0.0
	norm2 = 0.0
	do i = 1,iter+2
		ddot = ddot + gvec(i)*gvec_old(i)
		norm1 = norm1+gvec(i)*gvec(i)
		norm2 = norm2+gvec_old(i)*gvec_old(i)
	end do

!	print*,' testing ',ddot,norm1,norm2
	ddot = abs(ddot) / sqrt(norm1*norm2)

!........... STORE CURRENT SOLUTION AS 'OLD' VECTOR.....
	gvec_old = 0.d0
	do i = 1,iter+2
		gvec_old(i)=gvec(i)
	end do

	return

end subroutine resolvent

!======================================================================

end module lanczos_util

!/////////////////////////////////////////////////////////////////////

module convergence
	use lanczos_info
	use lanczos_util
	use precisions
	use flagger
	implicit none
    integer :: ncheck_ex        ! # of extra states to check for convergence
    integer :: thick_restart_factor , thick_restart_min

!...... CONVERGENCE IN LANCZOS........
  logical :: ask_converge_test = .false.        !  flag to allow question about convergence
  integer :: converge_test_def = 0
  integer :: converge_test
!....... THERE ARE SEVERAL OPTIONS FOR CONVERGENCE TEST.....
!        FOR LANCZOS BETWEEN TWO ITERATIONS:
!
!    CONVERGE_TEST = 0   test on average difference in energies
!                  = 1   test on max difference     in energies
!                  = 2   test on avg difference in wavefunctions
!                  = 3   test on max difference in wavefunctions
!    NOTE: 2,3 are good for observables, but could be problematic
!    if there are degeneracies in the spectrum
!
   real,parameter :: converge_tol0_def =  1.0e-3
   real,parameter :: converge_tol1_def =  5.0e-4
   real,parameter :: converge_tol2_def =  1.0e-5
   real,parameter :: converge_tol3_def =  2.0e-5
   real :: converge_tol

!   real    :: ediff0,ediff1
   real(8)     :: vdiff2,vdiff3,vdot

contains

!======================================================================
!
!
!  UNITED CONVERGENCE TEST
!  ADDED in 7.8.8 by CWJ
!

subroutine convergence_test(iter,i_thick_restart,finished)
	use io,only:greenflag
	use nodeinfo
	use butil_mod
	use bmpi_mod
	use localblocks,only:dimblock,block_flag
	implicit none
	integer,intent(IN) :: iter
	integer,intent(IN) :: i_thick_restart

	integer :: green_skip_iter = 5 ! check convergence on every 5th one
	integer :: green_skip_print = 10
!	real :: green_tol = 1.e-5
	logical   :: finished
	integer :: i,j,k,ierr
    real(kind=egv_prec), pointer ::  epoint(:)

	real(8) :: ddot

	if(greenflag .and. mod(iter,green_skip_iter)==0 )then
		call resolvent(iter,ddot)
		if(mod(iter,green_skip_print)==0 .and. iproc==0)then
			print*,iter,' Convergence of resolvent on pivot: ',ddot

		end if
		if(abs(ddot-1.0)< green_tol)then
			finished=.true.
			if(iproc==0)print*,' Within tolerance (green_tol) of ',green_tol
		end if

		if(iter==niter)then
			if(iproc==0)print*,' Did not reach desired tolerance (green_tol) of ',green_tol
			finished=.true.
		end if
        call BMPI_BCAST(finished,1,0,icomm,ierr)

		if(finished)return

	end if

	if(lanczchar /= 'tx' .or. txchoice==1)then
		epoint =>e
	else
		epoint => etx
	end if

!----------------- NORMAL CHECK -------------------------

     if ( iter == niter ) then
        finished = .true.

        if ( .not. fixed_iter .and. iproc == 0 ) then
           write(6,*)' Did not converge fully '
        end if
		return
     end if
     if ( .not. fixed_iter ) then

        if ( iter > ncheck  .and. (.not.skipeig .or. (skipeig .and. flyeig) .or. &
		   (block_flag==1 .and. iter > dimblock))) then
!................ CHECK DIFFERENCES IN ENERGY..........................
           ediff0 = 0.0
           ediff1 = 0.0
           do j = 1,ncheck
              ediff0 = ediff0 + abs(eold(j) - epoint(j) )
              ediff1 = bmax(ediff1, abs(  real( eold(j) -epoint(j),kind=4) ) )
           end do
           ediff0 = ediff0/sqrt(float(nkeep+1))*(1+i_thick_restart)

!................ CHECK DIFFERENCES IN VECTORS......................
           if(altconverge)then
              vdiff2 = 0.0d0
              vdiff3 = 0.0d0

              do j = 1,ncheck
                  vdot = 0.d0

                  do k = 1,iter-1
                      vdot = vdot + real(eiglvec(k,j),kind=8)*real(eiglvecold(k,j),kind=8)
                  end do
                  vdiff2 = vdiff2 + (abs(vdot)-1.d0)**2
                  vdiff3 = bmax( vdiff3, dabs(abs(vdot)-1.d0) )

              end do
              vdiff2 = dsqrt(vdiff2)/float(ncheck)
           endif

!................. CONVERGE CHECK.............
           if ( (flyeig .or. (block_flag==1 .and. iter > dimblock)  ).and. iproc == 0 ) then
                select case (converge_test)
                case (0)
		write(6,777)ediff0,converge_tol
                case (1)
		write(6,777)ediff1,converge_tol
                case (2)
		write(6,778)vdiff2,converge_tol
                case (3)
		write(6,778)vdiff3,converge_tol
                end select
           end if
777        format('      (energy convergence  ',f8.5,' > criterion ',f8.5,')')
778        format('      (wavefn convergence  ',f8.5,' > criterion ',f8.5,')')
           select case (converge_test)
           case ( 0)
              if ( ediff0 < converge_tol ) finished = .true.
           case ( 1)
              if ( ediff1 < converge_tol ) finished = .true.
           case ( 2)
              if ( vdiff2 < converge_tol ) finished = .true.
           case ( 3)
              if ( vdiff3 < converge_tol ) finished = .true.
           end select
        end if

        do j = 1,ncheck
           eold(j) = epoint(j)
           if(altconverge)then
                do k = 1,iter
                  eiglvecold(k,j) = eiglvec(k,j)
                end do
           end if
        end do

     end if
     call BMPI_BCAST(finished,1,0,icomm,ierr)
     return

end subroutine convergence_test

end module convergence

max_stars_repo_name: noctildon/pyBigstick

max_issues_repo_name: noctildon/pyBigstick

max_forks_repo_name: noctildon/pyBigstick


!        Generated by TAPENADE     (INRIA, Ecuador team)
!  tapenade 3.x
!
!  Differentiation of buoy in reverse (adjoint) mode (with options noISIZE i8):
!   gradient     of useful results: *temp *propunit *pres buoy
!   with respect to varying inputs: *temp *propunit *pres
!   Plus diff mem management of: temp:in propunit:in pres:in
! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!>    @brief calculate buoyancy for head equation
!>    @param[in] i cell index, direction I0
!>    @param[in] j cell index, direction J0
!>    @param[in] k cell index, direction K0
!>    @param[in] ismpl local sample index
!>    @return buoyancy
!>    @details
!>calculate buoyancy for head equation\n
!>sign convention: negative for positive buoyancy\n
SUBROUTINE BUOY_AD(i, j, k, ismpl, buoy_adv)
  use arrays

  USE ARRAYS_AD

  USE MOD_GENRL
  USE MOD_FLOW
  IMPLICIT NONE
  double precision :: buoy_adv
  INTEGER :: ismpl
  INTEGER :: i, j, k
  DOUBLE PRECISION :: rhor, rhav, hh, h0, h1, prod, summ
  DOUBLE PRECISION :: rhor_ad, rhav_ad, hh_ad, h0_ad, h1_ad, prod_ad, &
& summ_ad
  EXTERNAL RHOF, KZ, VISF
  EXTERNAL RHOF_AD, KZ_AD, VISF_AD
  DOUBLE PRECISION :: RHOF, KZ, VISF
  INTEGER :: arg1
  DOUBLE PRECISION :: result1
  DOUBLE PRECISION :: result1_ad
  DOUBLE PRECISION :: result2
  DOUBLE PRECISION :: result2_ad
  DOUBLE PRECISION :: result3
  DOUBLE PRECISION :: result3_ad
  INTEGER :: arg2
  INTEGER :: arg3
  DOUBLE PRECISION :: temporary_ad
  INTEGER :: branch
  DOUBLE PRECISION :: buoy
  arg1 = k + 1
  result1 = RHOF(i, j, arg1, ismpl)
  result2 = RHOF(i, j, k, ismpl)
  rhav = 0.5d0*(result1+result2)
  rhor = (rhav-rref)/rref
  hh = 0.d0
  result1 = KZ(i, j, k, ismpl)
  result2 = RHOF(i, j, k, ismpl)
  result3 = VISF(i, j, k, ismpl)
  h0 = result1*result2*grav/result3
  arg1 = k + 1
  CALL PUSHREAL8(result1)
  result1 = KZ(i, j, arg1, ismpl)
  arg2 = k + 1
  CALL PUSHREAL8(result2)
  result2 = RHOF(i, j, arg2, ismpl)
  arg3 = k + 1
  CALL PUSHREAL8(result3)
  result3 = VISF(i, j, arg3, ismpl)
  h1 = result1*result2*grav/result3
  summ = h0 + h1
  prod = h0*h1
  IF (summ .GT. 0.d0) THEN
    hh = 2.0d0*prod/summ
    CALL PUSHCONTROL1B(0)
  ELSE
    CALL PUSHCONTROL1B(1)
  END IF
  hh_ad = rhor*buoy_adv
  rhor_ad = hh*buoy_adv
  CALL POPCONTROL1B(branch)
  IF (branch .EQ. 0) THEN
    temporary_ad = 2.0d0*hh_ad/summ
    prod_ad = temporary_ad
    summ_ad = -(prod*temporary_ad/summ)
  ELSE
    prod_ad = 0.D0
    summ_ad = 0.D0
  END IF
  h0_ad = h1*prod_ad + summ_ad
  h1_ad = h0*prod_ad + summ_ad
  temporary_ad = grav*h1_ad/result3
  result1_ad = result2*temporary_ad
  result2_ad = result1*temporary_ad
  result3_ad = -(result1*result2*temporary_ad/result3)
  CALL POPREAL8(result3)
  CALL VISF_AD(i, j, arg3, ismpl, result3_ad)
  CALL POPREAL8(result2)
  CALL RHOF_AD(i, j, arg2, ismpl, result2_ad)
  CALL POPREAL8(result1)
  CALL KZ_AD(i, j, arg1, ismpl, result1_ad)
  temporary_ad = grav*h0_ad/result3
  result1_ad = result2*temporary_ad
  result2_ad = result1*temporary_ad
  result3_ad = -(result1*result2*temporary_ad/result3)
  CALL VISF_AD(i, j, k, ismpl, result3_ad)
  CALL RHOF_AD(i, j, k, ismpl, result2_ad)
  CALL KZ_AD(i, j, k, ismpl, result1_ad)
  rhav_ad = rhor_ad/rref
  result1_ad = 0.5d0*rhav_ad
  result2_ad = 0.5d0*rhav_ad
  CALL RHOF_AD(i, j, k, ismpl, result2_ad)
  arg1 = k + 1
  CALL RHOF_AD(i, j, arg1, ismpl, result1_ad)
END SUBROUTINE BUOY_AD


text:
program gen_be_stage3

   use da_control, only : stderr, stdout, filename_len, vertical_ip
   use da_gen_be, only : da_eof_decomposition_test, da_eof_decomposition, &
      da_transform_vptovv, da_create_bins
   use da_reporting, only : da_error
   use da_tools_serial, only : da_get_unit, da_advance_cymdh

   implicit none

   character*10        :: start_date, end_date       ! Starting and ending dates.
   character*10        :: date, new_date             ! Current date (ccyymmddhh).
   character*10        :: variable                   ! Variable name
   character(len=filename_len)        :: filename                   ! Input filename.
   character*2         :: ck                         ! Level index -> character.
   character*3         :: ce                         ! Member index -> character.
   integer             :: ni, nj, nk                 ! Dimensions read in.
   integer             :: sdate, cdate, edate        ! Starting, current ending dates.
   integer             :: interval                   ! Period between dates (hours).
   integer             :: ne                         ! Number of ensemble members.
   integer             :: i, j, k, k1, k2, b, member ! Loop counters.
   integer             :: bin_type                   ! Type of bin to average over.
   integer             :: num_bins                   ! Number of bins (3D fields).
   integer             :: num_bins2d                 ! Number of bins (2D fields).
   integer             :: ios                        ! I/O status for file read
   real                :: lat_min, lat_max           ! Used if bin_type = 2 (degrees).
   real                :: binwidth_lat               ! Used if bin_type = 2 (degrees).
   real                :: hgt_min, hgt_max           ! Used if bin_type = 2 (m).
   real                :: binwidth_hgt               ! Used if bin_type = 2 (m).
   real                :: inv_nij                    ! 1 / (ni*nj).
   real                :: mean_field                 ! Mean field.
   real                :: coeffa, coeffb             ! Accumulating mean coefficients.
   logical             :: first_time                 ! True if first file.
   logical             :: testing_eofs               ! True if testing EOF decomposition.
   logical             :: use_global_eofs            ! True if projected data uses global EOFs.
   logical             :: data_on_levels             ! True if output level data (diagnostic).
   logical             :: twod_field                 ! True if 2D field.
   logical             :: allow_missing_dates        ! If data from stage 1 is not contiguous, attempt to continue

   integer, allocatable:: bin(:,:,:)                 ! Bin assigned to each 3D point.
   integer, allocatable:: bin2d(:,:)                 ! Bin assigned to each 2D point.
   integer, allocatable:: bin_pts2d(:)               ! Number of points in bin (2D fields).
   real, allocatable   :: latitude(:,:)              ! Latitude (degrees, from south).
   real, allocatable   :: height(:,:,:)              ! Height field.
   real, allocatable   :: field(:,:,:)               ! Input field.
   real, allocatable   :: field_out(:,:,:)          ! Field projected into EOF space.
   real, allocatable   :: vertical_wgt(:,:,:)        ! Inner product for vertical transform.
   real, allocatable   :: bv(:,:,:)                  ! Vertical BE for this time.
   real, allocatable   :: work(:,:)                  ! EOF work array.
   real, allocatable   :: e_vec_loc(:,:,:)           ! Latitudinally varying eigenvectors.
   real, allocatable   :: e_val_loc(:,:)             ! Latitudinally varying eigenvalues.
   real*8, allocatable :: e_vec(:,:)                 ! Domain-averaged eigenvectors.
   real*8, allocatable :: e_val(:)                   ! Domain-averaged eigenvalues.
   real*8, allocatable :: evec(:,:,:)                ! Gridpoint eigenvectors.
   real*8, allocatable :: eval(:,:)                  ! Gridpoint sqrt(eigenvalues).

   namelist / gen_be_stage3_nl / start_date, end_date, interval, variable, &
                                 ne, bin_type, &
                                 lat_min, lat_max, binwidth_lat, &
                                 hgt_min, hgt_max, binwidth_hgt, &
                                 testing_eofs, use_global_eofs, data_on_levels, &
                                 allow_missing_dates

   integer :: ounit,iunit,namelist_unit

   stderr = 0
   stdout = 6

   write(6,'(a)')' [1] Initialize namelist variables and other scalars.'

   call da_get_unit(ounit)
   call da_get_unit(iunit)
   call da_get_unit(namelist_unit)

   vertical_ip = 0
   twod_field = .false.

   start_date = '2004030312'
   end_date = '2004033112'
   interval = 24
   variable = 'psi'
   ne = 1
   bin_type = 5
   lat_min = -90.0
   lat_max = 90.0
   binwidth_lat = 10.0
   hgt_min = 0.0
   hgt_max = 20000.0
   binwidth_hgt = 1000.0
   testing_eofs = .true.
   use_global_eofs = .true.
   data_on_levels = .false.
   allow_missing_dates = .false.

   open(unit=namelist_unit, file='gen_be_stage3_nl.nl', &
        form='formatted', status='old', action='read')
   read(namelist_unit, gen_be_stage3_nl)
   close(namelist_unit)

   read(start_date(1:10), fmt='(i10)')sdate
   read(end_date(1:10), fmt='(i10)')edate
   write(6,'(4a)')' Computing vertical error statistics for dates ', start_date, ' to ', end_date
   write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
   write(6,'(a,i8)')' Number of ensemble members at each time = ', ne

   date = start_date
   cdate = sdate
   first_time = .true.

   if ( .not. data_on_levels ) then ! Can bypass [2] and [3] if outputing on levels.

      write(6,'(2a)')' [2] Process fields for variable ', variable

      do while ( cdate <= edate )
         do member = 1, ne
            write(6,'(5a,i4)')'    Processing data for date ', date, ', variable ', trim(variable), &
                              ' and member ', member

            write(ce,'(i3.3)')member

            ! Read Full-fields:
            filename = 'fullflds'//'/'//date(1:10)//'.'//'fullflds'//'.e'//ce
            open (iunit, file = trim(filename), form = 'unformatted')
            read(iunit, iostat=ios)ni, nj, nk
            if (ios /= 0) then
               if (allow_missing_dates) then
                  write(6,'(a,a)')' WARNING: CAN NOT OPEN ',filename
                  write(6,'(a)')' Attempting to continue since allow_missing_dates = .true.'
                  cycle
               else
                  call da_error(__FILE__,__LINE__,(/"Could not open "//trim(filename)/))
               endif
            endif
            if ( first_time ) then
               write(6,'(a,3i8)')'    i, j, k dimensions are ', ni, nj, nk
               allocate( latitude(1:ni,1:nj) )
               allocate( height(1:ni,1:nj,1:nk) )
               allocate( bin(1:ni,1:nj,1:nk) )
               allocate( bin2d(1:ni,1:nj) )
            end if
            read(iunit)latitude
            read(iunit)height

            close(iunit)

            if ( first_time ) then
               ! Create and sort into bins:
               call da_create_bins( ni, nj, nk, bin_type, num_bins, num_bins2d, bin, bin2d, &
                                    lat_min, lat_max, binwidth_lat, &
                                    hgt_min, hgt_max, binwidth_hgt, latitude, height )

               allocate( bin_pts2d(1:num_bins2d) )
               bin_pts2d(:) = 0
            end if

            filename = trim(variable)//'/'//date(1:10)
            filename = trim(filename)//'.'//trim(variable)//'.e'//ce
            if ( trim(variable) == 'ps_u' .or. trim(variable) == 'ps' ) then ! 2D field
               twod_field = .true.
               filename = trim(filename)//'.01'
            end if

            open (iunit, file = filename, form='unformatted')
            read(iunit)ni, nj, nk

            if ( first_time ) then
               inv_nij = 1.0 / real(ni*nj)
               allocate( field(1:ni,1:nj,1:nk) )
               allocate( bv(1:nk,1:nk,1:num_bins2d) )
               bv(:,:,:) = 0.0
               first_time = .false.
            end if

            read(iunit)field
            close(iunit)

            ! Remove mean field:
            do k = 1, nk
               mean_field = sum(field(1:ni,1:nj,k)) * inv_nij
               field(1:ni,1:nj,k) = field(1:ni,1:nj,k) - mean_field
            end do

            do j = 1, nj
               do i = 1, ni
                  b = bin2d(i,j)
                  bin_pts2d(b) = bin_pts2d(b) + 1
                  coeffa = 1.0 / real(bin_pts2d(b))
                  coeffb = real(bin_pts2d(b)-1) * coeffa
                  do k1 = 1, nk
                     do k2 = 1, k1
                        bv(k1,k2,b) = coeffb * bv(k1,k2,b) + coeffa * field(i,j,k1) * field(i,j,k2)
                     end do
                  end do
               end do
            end do
         end do  ! End loop over ensemble members.

         ! Calculate next date:
         call da_advance_cymdh( date, interval, new_date )
         date = new_date
         read(date(1:10), fmt='(i10)')cdate
      end do     ! End loop over times.

      !  Fill in upper-right part of BE matrix by symmetry:

      do b = 1, num_bins2d
         do k1 = 1, nk
            do k2 = k1+1, nk ! Symmetry.
               bv(k1,k2,b) = bv(k2,k1,b)
            end do
         end do
      end do

      write(6,'(2a)')' [3] Calculate eigenvectors and eigenvalues for variable ', variable

      allocate( work(1:nk,1:nk) )
      allocate( e_vec_loc(1:nk,1:nk,1:num_bins2d) )
      allocate( e_val_loc(1:nk,1:num_bins2d) )
      allocate( e_vec(1:nk,1:nk) )
      allocate( e_val(1:nk) )

      !  Latitudinally varying BE decomposition:
      do b = 1, num_bins2d
         write(6,'(2(a,i6))')' Calculate eigenvectors and eigenvalues for bin ', b, &
                             ' of ', num_bins2d

         work(1:nk,1:nk) = bv(1:nk,1:nk,b)
         call da_eof_decomposition( nk, work, e_vec, e_val )
         e_vec_loc(1:nk,1:nk,b) = e_vec(1:nk,1:nk)
         e_val_loc(1:nk,b) = e_val(1:nk)
      end do

      !  Domain-averaged BE decomposition:
      work(1:nk,1:nk) = 0.0
      do b = 1, num_bins2d
         work(1:nk,1:nk) = work(1:nk,1:nk) + bv(1:nk,1:nk,b)
      end do
      work(1:nk,1:nk) = work(1:nk,1:nk) / real( num_bins2d )

      call da_eof_decomposition( nk, work, e_vec, e_val )

      if ( testing_eofs ) then
         call da_eof_decomposition_test( nk, work, e_vec, e_val )
      end if

      !  Output eigenvectors, eigenvalues for use in WRF_Var:
      filename = 'gen_be_stage3.'//trim(variable)//'.dat'
      open (ounit, file = filename, form='unformatted')
      write(ounit)variable
      write(ounit)nk, num_bins2d
      write(ounit)e_vec
      write(ounit)e_val
      write(ounit)e_vec_loc
      write(ounit)e_val_loc
      close(ounit)

      !  Decide on local or domain-averaged EOFs for horizontal correlations:
      if ( use_global_eofs ) then
         do b = 1, num_bins2d
            e_vec_loc(1:nk,1:nk,b) = e_vec(1:nk,1:nk)
            e_val_loc(1:nk,b) = e_val(1:nk)
         end do
      end if

      !  Map binned eigenvectors to x, y grid, and take sqrt(this is used in WRF-Var):
      allocate( evec(1:nj,1:nk,1:nk) )
      allocate( eval(1:nj,1:nk) )

      do j = 1, nj
         do i = 1, ni
            b = bin2d(i,j)
            evec(j,1:nk,1:nk) = e_vec_loc(1:nk,1:nk,b)
            eval(j,1:nk) = sqrt(e_val_loc(1:nk,b))
         end do
      end do

   end if ! End bypass [2] and [3] if outputing on levels.

   write(6,'(2a)')' [4] Transform perturbations (or not), and output.'

   date = start_date
   cdate = sdate
   first_time = .true.

   if ( .not. twod_field ) then
      do while ( cdate <= edate )
         do member = 1, ne
            write(6,'(5a,i4)')'    Date = ', date, ', variable ', trim(variable), &
                              ' and member ', member

            write(ce,'(i3.3)')member

            filename = trim(variable)//'/'//date(1:10)
            filename = trim(filename)//'.'//trim(variable)//'.e'//ce

            open (iunit, file = trim(filename), form = 'unformatted')
            read(iunit, iostat = ios)ni, nj, nk
            if (ios /= 0) then
               if (allow_missing_dates) then
                  write(6,'(a,a)')' WARNING: CAN NOT OPEN ',filename
                  write(6,'(a)')' Attempting to continue since allow_missing_dates = .true.'
                  cycle
               else
                  call da_error(__FILE__,__LINE__,(/"Could not open "//trim(filename)/))
               endif
            endif

            if ( first_time ) then
               if ( data_on_levels) allocate( latitude(1:ni,1:nj) )   ! Not allocated earlier.
               if ( data_on_levels) allocate( field(1:ni,1:nj,1:nk) ) ! Not allocated earlier.
               allocate( field_out(1:ni,1:nj,1:nk) )
               allocate( vertical_wgt(1:ni,1:nj,1:nk) )
               vertical_wgt(1:ni,1:nj,1:nk) = 1.0 ! vertical_ip = 0 hardwired.
               first_time = .false.
            end if

            read(iunit)field
            close(iunit)

            if ( data_on_levels ) then
               ! Keep data on vertical levels:
               field_out(:,:,:) = field(:,:,:)
            else
               ! Project fields onto vertical modes:
               call da_transform_vptovv( evec, eval, vertical_wgt, &
                  field, field_out, nk, &
                  1, nk, & ! WRF ids, ide etc.
                  1, ni, 1, nj, 1, nk, & ! WRF ims, ime etc.
                  1, ni, 1, nj, 1, nk )  ! WRF its, ite etc.
            end if

            ! Output fields (split into 2D files to allow parallel horizontal treatment):

            do k = 1, nk
               write(ck,'(i2)')k
               if ( k < 10 ) ck = '0'//ck(2:2)
               ! Assumes variable directory has been created by script:
               filename = trim(variable)//'/'//date(1:10)//'.'//trim(variable)
               filename = trim(filename)//'.e'//ce//'.'//ck
               open (ounit, file = filename, form='unformatted')
               write(ounit)ni, nj, k
               write(ounit)field_out(1:ni,1:nj,k)
               close(ounit)
            end do
         end do ! End loop over members.

         ! Calculate next date:
         call da_advance_cymdh( date, interval, new_date )
         date = new_date
         read(date(1:10), fmt='(i10)')cdate
      end do
   end if

end program gen_be_stage3

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb