PROGRAM control

IMPLICIT NONE
  INTEGER, PARAMETER :: DP = selected_real_kind(14,200)

  !DECLARATIONS
  INTEGER  :: npwx, npwq, ik, lter, npol
  INTEGER  :: nbnd, ikk, nbnd_occ(1), iw

  COMPLEX (DP), ALLOCATABLE :: dvpsi(:,:), dpsi(:,:), drhoscfs (:,:), dpsip(:,:), dpsim(:,:)
  COMPLEX (DP), ALLOCATABLE :: etc(:,:)
 
  COMPLEX(DP) :: cw

  REAL(DP), ALLOCATABLE :: fiu(:)
  
  REAL(DP), ALLOCATABLE :: et(:,:), h_diag(:,:)
  REAL(DP) :: thresh, anorm, averlt, dr2

  logical :: cgsolver, conv_root

  external :: ch_psi_all, cch_psi_all_fix, cg_psi
  cgsolver = .true.
  conv_root = .true.


  npol = 1
  nbnd = 1
  npwx = 10
  nbnd_occ (:) = 1  

!ALLOCATE VECTORS PSI and RHS:DVPSI

  allocate (dvpsi ( npwx*npol , nbnd))
  allocate (dpsi ( npwx*npol , nbnd))
  allocate (dpsim ( npwx*npol , nbnd))
  allocate (dpsip ( npwx*npol , nbnd))

  allocate (dpsip ( npwx*npol , nbnd))

!  Just want to operate these solvers on small matrices to see what happens. 
!  Could try some 10 by 10 systems with "analytic" solutions. Should give me a feel for exactly what the
!  difference between cBiCG and CG really is and what the difference is when we start looking at non-hermitian
!  matrices.

           if(cgsolver.eq..true.) then

               call cgsolve_all (ch_psi_all, cg_psi, et(1,ikk), dvpsi, dpsi, &
                             h_diag, npwx, npwq, thresh, ik, lter, conv_root, &
                             anorm, nbnd_occ(ikk), npol )

          !cw = dcmplx ( w, 0.00734981) ! 0.1 eV imaginary component

           else

               etc(:,:) = CMPLX( et(:,:), 0.0d0 , kind=DP)
               cw       =  CMPLX(0.0d0, fiu(iw), kind=DP)

               call cbcg_solve_fix(cch_psi_all_fix, cg_psi, etc(1,ikk), dvpsi, dpsip, h_diag, &
                npwx, npwq, thresh, ik, lter, conv_root, anorm, nbnd_occ(ikk), npol, cw, .true.)

               call cbcg_solve_fix(cch_psi_all_fix, cg_psi, etc(1,ikk), dvpsi, dpsim, h_diag, &
                npwx, npwq, thresh, ik, lter, conv_root, anorm, nbnd_occ(ikk), npol, -cw, .true.)

           endif
END PROGRAM control