rparesvec Subroutine

  subroutine rparesvec(q,w_l,w_r,v_l,v_r,ee,abd, &
                       ndsr,errors,tol,imax,tamm_dancoff)

    use precision, only: dp

    implicit none

    real(kind=dp), intent(inout), dimension(:,:) :: w_l, w_r, v_l, v_r
    real(kind=dp), intent(out), &
      dimension(ubound(w_l,1),ndsr,*) :: q
    real(kind=dp), intent(in), dimension(:) :: ee
    real(kind=dp), intent(in), dimension(:) :: abd
    integer, intent(in) :: ndsr
    real(kind=dp), intent(inout) :: errors(:)
    real(kind=dp), intent(in) :: tol
    integer, intent(in) :: imax
    logical :: tamm_dancoff

    real(kind=dp) :: er_l, er_r
    integer :: ivec, xvec_dim

    xvec_dim = ubound(w_l,1)

!   Residual vector W_l and W_r
    do ivec = 1, ndsr
      w_r(1:xvec_dim,ivec) = w_r(1:xvec_dim,ivec) - ee(ivec)*v_r(1:xvec_dim,ivec)
    end do
    if (.not.tamm_dancoff) then
      do ivec = 1, ndsr
        w_l(1:xvec_dim,ivec) = w_l(1:xvec_dim,ivec) - ee(ivec)*v_l(1:xvec_dim,ivec)
      end do
    end if

!   Norms of W
    errors = 0.0_dp

    if (tamm_dancoff) then
      do ivec = imax+1, ndsr
        errors(ivec) = dot_product(w_r(:,ivec),w_r(:,ivec))
        if (errors(ivec)>1.0_dp) then
          write(*,'(4X,"Large error detected")')
          write(*,'(4X,"State#=",I4)') ivec
          write(*,'(4X,"Error right =",1X,1P,E10.3)') errors(ivec)
          errors(ivec) = 0.0_dp
        end if

!       Get new vectors Q
        if (errors(ivec)>tol) then
          q(:,ivec,1) = w_r(:,ivec)/(ee(ivec)-abd(:))
        else
          q(:,ivec,1) = 0.0_dp
        end if
      end do
    else
      do ivec = imax+1, ndsr
        er_l = dot_product(w_l(:,ivec),w_l(:,ivec))
        er_r = dot_product(w_r(:,ivec),w_r(:,ivec))
        errors(ivec) = max(er_l,er_r)
        if (errors(ivec)>1.0_dp) then
          write(*,'(4X,"Large error detected")')
          write(*,'(4X,"State#=",I4)') ivec
          write(*,'(4X,"Error left/right =",1X,1P,E10.3,"/",1X,1P,E10.3)') er_l, er_r
          errors(ivec) = 0.0_dp
        end if

!       Get new vectors Q
        if (errors(ivec)>tol) then
          q(:,ivec,1) = 1.0d0/(ee(ivec)-abd(:))
          q(:,ivec,2) = q(:,ivec,1)

          q(:,ivec,1) = q(:,ivec,1)*w_l(:,ivec)
          q(:,ivec,2) = q(:,ivec,2)*w_r(:,ivec)
        else
          q(:,ivec,1:2) = 0.0_dp
        end if
      end do
    end if

  end subroutine rparesvec