text:
*----------------------------------------------------------------------*
      subroutine trace_op(xfac,casign,
     &     update,xret,type_xret,
     &     me_op,me_trop,me_troptmp,
     &     tra_op, tra_trop,
     &     iblkop,iblktrop,iblktroptmp,
     &     idoffop,idofftrop,
     &     cnt_info,reo_info,
     &     str_info,strmap_info,orb_info)
*----------------------------------------------------------------------*
*
*     calculate (partial) trace of operator:
*
*        OPtr(IC,IA) = sum_J OP(IC*JC,IA*JA)
*
*     andreas, may 2008
*
*----------------------------------------------------------------------*
      implicit none

      include 'routes.h'
      include 'contr_times.h'

      include 'opdim.h'
      include 'stdunit.h'
      include 'ioparam.h'
      include 'multd2h.h'
      include 'def_operator.h'
      include 'def_me_list.h'
      include 'def_orbinf.h'
      include 'def_graph.h'
      include 'def_strinf.h'
      include 'def_filinf.h'
      include 'def_strmapinf.h'
      include 'def_reorder_info.h'
      include 'def_contraction_info.h'
      include 'ifc_memman.h'
      include 'ifc_operators.h'
      include 'hpvxseq.h'

      integer, parameter ::
     &     ntest = 000

      logical, intent(in) ::
     &     update
      real(8), intent(in) ::
     &     xfac, casign
      real(8), intent(inout), target ::
     &     xret(1)
      type(contraction_info), target ::
     &     cnt_info
      integer, intent(in) ::
     &     type_xret,
     &     iblkop, iblktrop, iblktroptmp,
     &     idoffop,idofftrop
      logical, intent(in) ::
     &     tra_op, tra_trop
      type(me_list), intent(in) ::
     &     me_op, me_trop, me_troptmp
      type(strinf), intent(in) ::
     &     str_info
      type(strmapinf), intent(inout) ::
     &     strmap_info
      type(orbinf), intent(in) ::
     &     orb_info
      type(reorder_info), intent(in), target ::
     &     reo_info

      logical ::
     &     bufop, buftrop,
     &     first1, first2, first3, first4, first5,
     &     reo_trop, nonzero, ms_fix, fix_success
      integer ::
     &     mstop,msttrop,
     &     igamtop,igamttrop,
     &     nc_op, na_op,
     &     nc_trop, na_trop,
     &     nc_troptmp, na_troptmp,
     &     nc_cnt, na_cnt,
     &     nsym, ifree, lenscr, lenblock,
     &     idxst_op, idxst_trop,
     &     idxop, idxtrop,
     &     lenop, lentrop,
     &     mscmx_a, mscmx_c, msc_ac, msc_a, msc_c,
     &     igamc_ac, igamc_a, igamc_c,
     &     idxms, idxdis, lenmap, lblk_troptmp,
     &     idxdis_trop
      integer ::
     &     ncblk_op, nablk_op,
     &     ncblk_trop, nablk_trop, ncblk_troptmp, nablk_troptmp,
     &     ncblk_cnt, nablk_cnt,
     &     ncblk_trop_0, nablk_trop_0,
     &     ncblk_reo12,    nablk_reo12
      type(filinf), pointer ::
     &     ffop,ffop2,fftrop
      type(operator), pointer ::
     &     op, op2, trop, troptmp
      integer, pointer ::
     &     cinfo_opc(:,:),cinfo_opa(:,:),
     &     cinfo_tropc(:,:),
     &     cinfo_tropa(:,:),
     &     cinfo_troptmpc(:,:),
     &     cinfo_troptmpa(:,:),
     &     cinfo_cntc(:,:),cinfo_cnta(:,:),
     &     map_info_c(:),
     &     map_info_a(:),
     &     dmap_troptmpc(:), dmap_troptmpa(:),
     &     dmap_opc(:), dmap_opa(:)

      real(8) ::
     &     xnrm
      real(8) ::
     &     cpu, sys, cpu0, sys0, cpu00, sys00

      real(8), pointer ::
     &     xop(:), xtrop(:), xscr(:)
      real(8), pointer ::
     &     xbf(:), xbftr(:), xbftrtmp(:), xtropblk(:)

      integer ::
     &     msbnd(2,2), igambnd(2,2),
     &     ms_op_tr_a(2), ms_op_tr_c(2),
c dbg mh igam_op_tr_a(2) --> igam_op_tr_a(3)
     &     igam_op_tr_a(3), igam_op_tr_c(2)

      integer, pointer ::
     &     gmopdis_c(:), gmopdis_a(:),
     &     gmc_dis_c (:), gmc_dis_a (:),
     &     gmtropdis_c (:), gmtropdis_a (:),
     &     msopdis_c(:), msopdis_a(:),
     &     msc_dis_c (:), msc_dis_a (:),
     &     mstropdis_c(:), mstropdis_a(:),
     &     idxmsopdis_c(:), idxmsopdis_a(:),
     &     idxmsc_dis_c (:), idxmsc_dis_a (:),
     &     idxmstropdis_c(:), idxmstropdis_a(:),
     &     lstrcnt(:),
     &     lstrop(:),lstrtroptmp(:)

      integer, pointer ::
     &     map_tropcnta(:), map_tropcntc(:)

      integer, pointer ::
     &     ndis_op(:,:), d_gam_ms_op(:,:,:), gam_ms_op(:,:),
     &     ndis_troptmp(:,:), d_gam_ms_trop(:,:,:),
     &     gam_ms_trop(:,:),
     &     len_gam_ms_troptmp(:,:), len_d_gam_ms_troptmp(:,:,:)

      integer, pointer ::
     &     cinfo_reo12c(:,:), cinfo_reo12a(:,:),
     &     cinfo_trop_0c(:,:), cinfo_trop_0a(:,:),
     &     map_info_reo1c(:), map_info_reo1a(:),
     &     map_info_reo2c(:), map_info_reo2a(:)

      type(graph), pointer ::
     &     graphs(:)

      integer, external ::
     &     ielsum, ielprd, idx_msgmdst2, get_lenmap, idxlist,
     &     max_dis_blk, msa2idxms4op
      logical, external ::
     &     next_dist, next_msgamdist2,
     &     next_msgamdist_diag
      real(8), external ::
     &     ddot
c dbg
      real(8) xold
c dbg

      if (ntest.gt.0) then
        call write_title(lulog,wst_dbg_subr,
     &       'trace_op at work')
      end if

      op => me_op%op
      trop => me_trop%op
      troptmp => me_troptmp%op

      ffop => me_op%fhand
      fftrop => me_trop%fhand

      ms_fix=.false.
      if(me_op%fix_vertex_ms.or.me_trop%fix_vertex_ms)then
        ms_fix = me_op%fix_vertex_ms.and.me_trop%fix_vertex_ms
        if(.not.ms_fix) call quit(1,'trace_op',
     &                            'fix ms or not?')
      endif

      if (ntest.ge.10) then
        write(lulog,*) 'list(in):   ',trim(me_op%label),
     &       ' transp:',tra_op
        write(lulog,*) 'list(out):  ',trim(me_trop%label),
     &       ' transp:',tra_trop
        write(lulog,*) 'ffop:   ',ffop%name(1:len_trim(ffop%name))
        write(lulog,*) 'fftrop:',
     &       fftrop%name(1:len_trim(fftrop%name))
        write(lulog,*) 'xfac = ',xfac
        write(lulog,*) 'casign = ',casign
        if (type_xret.ne.0)
     &       write(lulog,*) 'xret on entry = ',xret(1)
        write(lulog,*) 'op: ',trim(op%name),
     &       ' block ',iblkop
        if (iblktrop.gt.0) then
          write(lulog,*) 'trop: ',trim(trop%name),
     &       ' block ',iblktrop
        else
          write(lulog,*) 'trop: scalar'
        end if
      end if

      ! flag whether non-zero contribution to trop occurred
      nonzero = .false.

      reo_trop = reo_info%n_op_reo.gt.0
      if (reo_trop.and..not.associated(reo_info%map_reo1c))
     &     call quit(1,'trace_op',
     &     'reo_info is not consistent')
      if (ntest.ge.10) write(lulog,*) 'reo_trop: ',reo_trop
      if (ntest.ge.10 .and. reo_trop) then
        write(lulog,*) 'troptmp: ',trim(troptmp%name),
     &       ' block ',iblktroptmp
      end if

      ! set pointers to cnt_info entries
      ncblk_op = cnt_info%ncblk_op1
      nablk_op = cnt_info%nablk_op1
      ncblk_cnt = cnt_info%ncblk_cnt
      nablk_cnt = cnt_info%nablk_cnt
      ncblk_trop = cnt_info%ncblk_op1op2
      nablk_trop = cnt_info%nablk_op1op2
      ncblk_troptmp = cnt_info%ncblk_op1op2tmp
      nablk_troptmp = cnt_info%nablk_op1op2tmp

      cinfo_opc => cnt_info%cinfo_op1c
      cinfo_opa => cnt_info%cinfo_op1a
      cinfo_tropc => cnt_info%cinfo_op1op2c
      cinfo_tropa => cnt_info%cinfo_op1op2a
      cinfo_troptmpc => cnt_info%cinfo_op1op2tmpc
      cinfo_troptmpa => cnt_info%cinfo_op1op2tmpa
      cinfo_cntc => cnt_info%cinfo_cntc
      cinfo_cnta => cnt_info%cinfo_cnta

      map_info_c => cnt_info%map_info_1c
      map_info_a => cnt_info%map_info_1a

      allocate(
     &       gmopdis_c(ncblk_op), gmopdis_a(nablk_op),
     &       gmtropdis_c(ncblk_trop), gmtropdis_a(nablk_trop),
     &       gmc_dis_c(ncblk_cnt), gmc_dis_a(nablk_cnt),
     &       msopdis_c(ncblk_op), msopdis_a(nablk_op),
     &       mstropdis_c(ncblk_trop), mstropdis_a(nablk_trop),
     &       msc_dis_c(ncblk_cnt), msc_dis_a(nablk_cnt),
     &       idxmsopdis_c(ncblk_op), idxmsopdis_a(nablk_op),
     &       idxmstropdis_c(ncblk_trop), idxmstropdis_a(nablk_trop),
     &       idxmsc_dis_c(ncblk_cnt), idxmsc_dis_a(nablk_cnt),
     &       lstrcnt(ncblk_cnt+nablk_cnt),
     &       lstrop(ncblk_op+nablk_op),
     &       lstrtroptmp(ncblk_troptmp+nablk_troptmp),
     &       dmap_troptmpc(ncblk_trop),dmap_troptmpa(nablk_trop),
     &       dmap_opc(ncblk_op),dmap_opa(nablk_op))

      call set_dis_tra_map(dmap_troptmpc,dmap_troptmpa,
     &     cinfo_troptmpc(1,3),cinfo_troptmpa(1,3),
     &                                ncblk_troptmp,nablk_troptmp)
      call set_dis_tra_map(dmap_opc,dmap_opa,
     &     cinfo_opc(1,3),cinfo_opa(1,3),ncblk_op,nablk_op)

      idxst_op = me_op%off_op_occ(iblkop) + 1
      lenop    = me_op%len_op_occ(iblkop)
      if (iblktrop.gt.0) then
        ! refers to reordered trop (if that makes a difference)
        idxst_trop = me_trop%off_op_occ(iblktrop) + 1
        lentrop    = me_trop%len_op_occ(iblktrop)
      else
        idxst_trop = 1
        lentrop = 1
      end if

      mstop = me_op%mst
      msttrop = me_trop%mst
      igamtop = me_op%gamt
      igamttrop = me_trop%gamt

      if (igamtop.ne.igamttrop)
     &     call quit(1,'trace_op','inconsistent symmetries')

      if (lenop.le.0.or.lentrop.le.0) then
        write(lulog,*)
     &       trim(op%name),' ',
     &       trim(trop%name)
        write(lulog,*) 'lenop, lentrop: ',
     &                  lenop, lentrop
ctest        call warn('trace_op','zero length for operator?')
        return
      end if

      ifree = mem_setmark('trace_op')

      if (ffop%buffered.and.ffop%incore(iblkop).ge.0) then
        bufop = .true.
        xop => ffop%buffer(idxst_op:)
      else
        bufop = .false.
        ! LOWER incore requirements:
        ! check for length of operator
        ! if length < 2 * blocks * da_reclen, try to alloc
        !   full block
        ! else get largest symmetry block
        ifree = mem_alloc_real(xbf,lenop,'xbf')
        xop => xbf
        call get_vec(ffop,xop,idoffop+idxst_op,
     &                          idoffop+idxst_op-1+lenop)
      end if

      if (ntest.ge.100) write(lulog,*) ' bufop: ',bufop

      ! get result vector as well (as we update)
      ! refers to reordered trop
      if (iblktrop.gt.0) then
        if (fftrop%buffered.and.fftrop%incore(iblktrop).ge.0) then
          buftrop = .true.
          xtrop => fftrop%buffer(idxst_trop:)
          if (.not.update) xtrop(1:lentrop) = 0d0
        else
          buftrop = .false.
          ! LOWER incore requirements:
          ! see above
          ifree = mem_alloc_real(xbftr,lentrop,'xbftr')
          xtrop => xbftr
          if (update) then
            ! read from disc
            call get_vec(fftrop,xtrop,idofftrop+idxst_trop,
     &                             idofftrop+idxst_trop-1+lentrop)
          else
            ! init with zero
            xtrop(1:lentrop) = 0d0
          end if
        end if
        if (ntest.ge.100) write(lulog,*) ' buftrop: ',buftrop
      else
        buftrop = .true.
        xtrop => xret
        if (ntest.ge.100) write(lulog,*) ' result is scalar '
      end if

      if (ntest.ge.1000) then
        ! this will work if all blocks incore, only:
        write(lulog,*) 'operator(in) (',trim(op%name),
     &                    ',list=',trim(me_op%label),')'
        call wrt_mel_buf(lulog,5,xop,me_op,iblkop,iblkop,
     &                  str_info,orb_info)
        if (iblktrop.gt.0) then
          write(lulog,*) 'operator(out) on entry (',trim(trop%name),
     &                                ',list=',trim(me_trop%label),')'

          call wrt_mel_buf(lulog,5,xtrop,me_trop,
     &                    iblktrop,iblktrop,
     &                    str_info,orb_info)
        end if
      end if

      if (reo_trop) then
        call quit(1,'trace_op','reo-route is not yet debugged')
        lblk_troptmp =
     &       max_dis_blk(0,me_troptmp,iblktroptmp,orb_info)
        ifree = mem_alloc_real(xbftrtmp,lblk_troptmp,'xbftrtmp')
      end if

      graphs => str_info%g

      ndis_op => me_op%off_op_gmox(iblkop)%ndis
      gam_ms_op => me_op%off_op_gmo(iblkop)%gam_ms
      d_gam_ms_op => me_op%off_op_gmox(iblkop)%d_gam_ms
      ndis_troptmp => me_troptmp%off_op_gmox(iblktroptmp)%ndis
      gam_ms_trop => me_trop%off_op_gmo(iblktrop)%gam_ms
      len_gam_ms_troptmp =>
     &                   me_troptmp%len_op_gmo(iblktroptmp)%gam_ms
      d_gam_ms_trop => me_trop%off_op_gmox(iblktrop)%d_gam_ms
      len_d_gam_ms_troptmp =>
     &                  me_troptmp%len_op_gmox(iblktroptmp)%d_gam_ms

      call sum_occ(nc_op,cinfo_opc,ncblk_op)
      call sum_occ(na_op,cinfo_opa,nablk_op)
      call sum_occ(nc_trop,cinfo_tropc,ncblk_trop)
      call sum_occ(na_trop,cinfo_tropa,nablk_trop)
      call sum_occ(nc_troptmp,cinfo_troptmpc,ncblk_troptmp)
      call sum_occ(na_troptmp,cinfo_troptmpa,nablk_troptmp)
c dbg
      if (na_trop.ne.na_troptmp)
     &     call quit(1,'trace_op','unexpected 1a')
      if (nc_trop.ne.nc_troptmp) then
        write(lulog,*) 'TROP (C)   : ',nc_trop,
     &       ' <- ',cinfo_tropc(1:ncblk_trop,1)
        write(lulog,*) 'TROPTMP (C): ',nc_troptmp,
     &       ' <- ',cinfo_troptmpc(1:ncblk_troptmp,1)
        call quit(1,'trace_op','unexpected 1b')
      end if
c dbg
      call sum_occ(nc_cnt,cinfo_cntc,ncblk_cnt)
      call sum_occ(na_cnt,cinfo_cnta,nablk_cnt)

c dbg
c      print *,'map(C):'
c      call print_mapinfo(lulog,map_info_c,ncblk_op)
c      print *,'map(A):'
c      call print_mapinfo(lulog,map_info_a,nablk_op)
c dbg
      ! set up maps (if necessary)
      call strmap_man_c(1,lenmap,
     &     cinfo_cntc(1,2),ncblk_cnt,
     &     cinfo_tropc(1,2),ncblk_trop,
     &     cinfo_opc(1,2),ncblk_op,map_info_c,
     &     str_info,strmap_info,orb_info)
      ifree = mem_alloc_int(map_tropcntc,lenmap,'opmap_c')
      call strmap_man_c(1,lenmap,
     &     cinfo_cnta(1,2),nablk_cnt,
     &     cinfo_tropa(1,2),nablk_trop,
     &     cinfo_opa(1,2),nablk_op,map_info_a,
     &     str_info,strmap_info,orb_info)
      ifree = mem_alloc_int(map_tropcnta,lenmap,'opmap_a')

      if (reo_trop) then
        cinfo_reo12c => reo_info%cinfo_reo_c
        cinfo_reo12a => reo_info%cinfo_reo_a
        cinfo_trop_0c => reo_info%cinfo_opreo0c
        cinfo_trop_0a => reo_info%cinfo_opreo0a
        map_info_reo1c  => reo_info%map_reo1c
        map_info_reo1a  => reo_info%map_reo1a
        map_info_reo2c  => reo_info%map_reo2c
        map_info_reo2a  => reo_info%map_reo2a
        ncblk_trop_0  = reo_info%ncblk_reo0
        nablk_trop_0  = reo_info%nablk_reo0
        ncblk_reo12  = reo_info%ncblk_reo
        nablk_reo12  = reo_info%nablk_reo
        call strmap_man_c(2,lenmap,
     &     cinfo_trop_0c(1,2),ncblk_trop_0,
     &     cinfo_reo12c(1,2),ncblk_reo12,
     &     cinfo_troptmpc(1,2),ncblk_troptmp,map_info_reo1c,
     &     str_info,strmap_info,orb_info)
        call strmap_man_c(2,lenmap,
     &     cinfo_reo12a(1,2),nablk_reo12,
     &     cinfo_trop_0a(1,2),nablk_trop_0,
     &     cinfo_troptmpa(1,2),nablk_troptmp,map_info_reo1a,
     &     str_info,strmap_info,orb_info)
        call strmap_man_c(2,lenmap,
     &     cinfo_trop_0c(1,2),ncblk_trop_0,
     &     cinfo_reo12c(1,2),ncblk_reo12,
     &     cinfo_tropc(1,2),ncblk_trop,map_info_reo2c,
     &     str_info,strmap_info,orb_info)
        call strmap_man_c(2,lenmap,
     &     cinfo_reo12a(1,2),nablk_reo12,
     &     cinfo_trop_0a(1,2),nablk_trop_0,
     &     cinfo_tropa(1,2),nablk_trop,map_info_reo2a,
     &     str_info,strmap_info,orb_info)
      end if

      ! minimum Ms(A) for ...
      msbnd(1,1) = -na_op ! operator
      msbnd(1,2) = -na_trop ! trace
      ! maximum Ms(A) for ...
      msbnd(2,1) = -msbnd(1,1)
      msbnd(2,2) = -msbnd(1,2)
      ! max |Ms| for ...
      mscmx_a = na_cnt!+nc_cnt ! C(A)
      mscmx_c = nc_cnt!+na_cnt ! C(C)
      ! minimum IRREP for operators
      igambnd(1,1) = 1
      igambnd(1,2) = 1
      ! maximum IRREP
      nsym = orb_info%nsym
      igambnd(2,1) = nsym
      igambnd(2,2) = nsym
      ! loop Ms-cases of (Op(A),TrOp(A))
      first1 = .true.
c dbg
      xold = xtrop(1)
c dbg
      ms_loop: do
        if (first1) then
          first1 = .false.
          ! initial Ms distribution
          ms_op_tr_a(1:2) = msbnd(2,1:2)
        else
          ! next Ms distribution
          if (.not.next_dist(ms_op_tr_a,2,msbnd,-2)) exit
        end if

        ms_op_tr_c(1) = ms_op_tr_a(1) + mstop
        ms_op_tr_c(2) = ms_op_tr_a(2) + msttrop
        if (abs(ms_op_tr_c(1)).gt.nc_op) cycle ms_loop
        if (abs(ms_op_tr_c(2)).gt.nc_trop) cycle ms_loop

        msc_c = ms_op_tr_c(1) - ms_op_tr_c(2)
        msc_a = ms_op_tr_a(1) - ms_op_tr_a(2)
        if (msc_c.ne.msc_a)
     &       call quit(1,'trace_op','trap 1')

        if (mscmx_a.lt.abs(msc_a)) cycle ms_loop

        ! loop IRREP cases of (Op(A),TrOp(A))
        first2 = .true.
        gam_loop: do
          if (first2) then
            first2 = .false.
            ! initial IRREP distribution
            igam_op_tr_a(1:2) = igambnd(1,1:2)
          else
            ! next IRREP distribution
            if (.not.next_dist(igam_op_tr_a,2,igambnd,+1)) exit
          end if
          ! set up start addresses
          ! need to be modified, if more than one distribution
          ! exists, see below
          idxms = msa2idxms4op(ms_op_tr_a(1),mstop,na_op,nc_op)
          idxop = gam_ms_op(igam_op_tr_a(1),idxms) + 1
     &             - idxst_op+1
          idxms = msa2idxms4op(ms_op_tr_a(2),msttrop,na_trop,nc_trop)
          ! relevant for case where no reordering necessary
          ! then we have: troptmp == trop
          if (iblktrop.gt.0)
     &           idxtrop = gam_ms_trop(igam_op_tr_a(2),idxms) + 1
     &                - idxst_trop+1
          if (reo_trop)
     &           lblk_troptmp=len_gam_ms_troptmp(igam_op_tr_a(2),idxms)
          if (iblktrop.eq.0) idxtrop = 1

          igam_op_tr_c(1) = multd2h(igam_op_tr_a(1),igamtop)
          igam_op_tr_c(2) = multd2h(igam_op_tr_a(2),igamttrop)

          igamc_a = multd2h(igam_op_tr_a(1),igam_op_tr_a(2))
          igamc_c = multd2h(igam_op_tr_c(1),igam_op_tr_c(2))
c dbg
c          print *,'msc_c, msc_a:     ',msc_c, msc_a
c          print *,'igamc_c, igamc_a: ',igamc_c, igamc_a
c dbg

          ! loop over distributions of current Ms and IRREP
          ! of Aex and Cex (=trOP(A),trOP(C)) over ngastypes
          first3 = .true.
          caex_loop: do
            if (.not.next_msgamdist2(first3,
     &          mstropdis_c,mstropdis_a,gmtropdis_c,gmtropdis_a,
     &          ncblk_trop, nablk_trop,
     &          cinfo_tropc,cinfo_tropa,
     &          ms_op_tr_c(2),ms_op_tr_a(2),
     &          igam_op_tr_c(2),igam_op_tr_a(2),nsym,
     &          ms_fix,fix_success))
     &          exit caex_loop
            first3 = .false.
            if(ms_fix.and..not.fix_success)cycle caex_loop

            call ms2idxms(idxmstropdis_c,mstropdis_c,
     &           cinfo_tropc,ncblk_trop)
            call ms2idxms(idxmstropdis_a,mstropdis_a,
     &           cinfo_tropa,nablk_trop)
c dbg
c            print *,'mstropdis_c: ',mstropdis_c(1:ncblk_trop)
c            print *,'mstropdis_a: ',mstropdis_a(1:nablk_trop)
c            print *,'gmtropdis_c: ',gmtropdis_c(1:ncblk_trop)
c            print *,'gmtropdis_a: ',gmtropdis_a(1:nablk_trop)
c dbg

            call set_len_str(lstrtroptmp,ncblk_trop,nablk_trop,
     &           graphs,
     &           cinfo_tropc(1,2),idxmstropdis_c,
     &             gmtropdis_c,cinfo_tropc(1,3),
     &           cinfo_tropa(1,2),idxmstropdis_a,
     &             gmtropdis_a,cinfo_tropa(1,3),
     &           hpvxseq,.false.)

            ! test C and A separately to avoid overflow
            if ( ncblk_trop+nablk_trop.gt.0 .and.
     &               idxlist(0,lstrtroptmp,
     &                          ncblk_trop+nablk_trop,1).gt.0)
     &               cycle

            idxms =  msa2idxms4op(ms_op_tr_a(2),msttrop,na_trop,nc_trop)
            if (iblktroptmp.gt.0.and.
     &           ndis_troptmp(igam_op_tr_a(2),idxms).gt.1) then
              idxdis =
     &             idx_msgmdst2(.true.,
     &                   iblktroptmp,idxms,igam_op_tr_a(2),
     &                   cinfo_troptmpc,idxmstropdis_c,
     &                              gmtropdis_c,ncblk_troptmp,
     &                   cinfo_troptmpa,idxmstropdis_a,
     &                              gmtropdis_a,nablk_troptmp,
     &                   tra_trop,dmap_troptmpc,dmap_troptmpa,
     &                              me_troptmp,nsym)
              idxdis_trop = idxdis

              ! relevant for case w/o reordering
              ! then we have troptmp == trop
              idxtrop =
     &             d_gam_ms_trop(idxdis,igam_op_tr_a(3),idxms)+1
     &                   - idxst_trop+1

              if (reo_trop)
     &             lblk_troptmp =
     &             len_d_gam_ms_troptmp(idxdis,igam_op_tr_a(3),idxms)

            end if

            if (.not.reo_trop) then
              ! direct update of result block
c              xtropblk => xtrop(idxtrop:)
              ! somehow, "=> xxx(i:)" may lead to seg. faults(?)
              xtropblk => xtrop(idxtrop:lentrop)
            else
              ! put result to intermediate buffer
              xtropblk => xbftrtmp
              xbftrtmp(1:lblk_troptmp) = 0d0
            end if

            ! loop over distributions of current Ms and IRREP
            ! of AC and CC over ngastypes
            first5 = .true.
            cac_loop: do
              if (.not.next_msgamdist_diag(first5,
     &              msc_dis_a,gmc_dis_a,
     &              nablk_cnt,
     &              cinfo_cnta,
     &              msc_a,igamc_a,nsym))
     &                 exit cac_loop
              first5 = .false.

              call ms2idxms(idxmsc_dis_a,msc_dis_a,
     &              cinfo_cnta,nablk_cnt)

              ! length of contraction
              call set_len_str(lstrcnt,ncblk_cnt,nablk_cnt,
     &              graphs,
     &              cinfo_cntc(1,2),idxmsc_dis_a,
     &                                  gmc_dis_a,cinfo_cntc(1,3),
     &              cinfo_cnta(1,2),idxmsc_dis_a,
     &                                  gmc_dis_a,cinfo_cnta(1,3),
     &              hpvxseq,.false.)

              if ( ncblk_cnt+nablk_cnt.gt.0 .and.
     &                   idxlist(0,lstrcnt,
     &                   ncblk_cnt+nablk_cnt,1).gt.0)
     &                   cycle cac_loop

              ! get Ms and IRREP distribution of op
              call merge_msgmdis(msopdis_c,gmopdis_c,
     &                                 ncblk_op,
     &                                 msc_dis_a,gmc_dis_a,
     &                                 mstropdis_c,gmtropdis_c,
     &                                 map_info_c)
              call merge_msgmdis(msopdis_a,gmopdis_a,
     &                                 nablk_op,
     &                                 msc_dis_a,gmc_dis_a,
     &                                 mstropdis_a,gmtropdis_a,
     &                                 map_info_a)

              call ms2idxms(idxmsopdis_c,msopdis_c,
     &              cinfo_opc,ncblk_op)
              call ms2idxms(idxmsopdis_a,msopdis_a,
     &              cinfo_opa,nablk_op)

              call set_len_str(
     &                   lstrop,ncblk_op,nablk_op,
     &                   graphs,
     &                   cinfo_opc(1,2),idxmsopdis_c,
     &                                    gmopdis_c,cinfo_opc(1,3),
     &                   cinfo_opa(1,2),idxmsopdis_a,
     &                                    gmopdis_a,cinfo_opa(1,3),
     &                   hpvxseq,.false.)

              if ( ncblk_op+nablk_op.gt.0 .and.
     &                   idxlist(0,lstrop,
     &                             ncblk_op+nablk_op,1).gt.0)
     &                   cycle cac_loop

              ! get distribution index
              idxms =  msa2idxms4op(ms_op_tr_a(1),mstop,na_op,nc_op)
              if (ndis_op(igam_op_tr_a(1),idxms).gt.1) then
                idxdis =
     &                idx_msgmdst2(.true.,
     &                     iblkop,idxms,igam_op_tr_a(1),
     &                     cinfo_opc,idxmsopdis_c,
     &                              gmopdis_c,ncblk_op,
     &                     cinfo_opa,idxmsopdis_a,
     &                              gmopdis_a,nablk_op,
     &                     tra_op,dmap_opc,dmap_opa,me_op,nsym)
c     &                     .false.,me_op,nsym)
                idxop =   d_gam_ms_op(idxdis,igam_op_tr_a(1),idxms) + 1
     &                     - idxst_op+1
              end if

              xnrm = ddot(ielprd(lstrop,
     &                   ncblk_op+nablk_op),
     &                   xop(idxop),1,xop(idxop),1)
              if (xnrm.lt.1d-28) cycle cac_loop

              ! if we get here, trop will change:
              nonzero = .true.

              ! get igrphcnt,igrphext1->igrphop map
              ! for given ms and irreps
              lenmap = get_lenmap(lstrcnt,lstrtroptmp,
     &                   map_info_c,ncblk_op)
              call get_strmap_blk_c(map_tropcntc,
     &                   ncblk_cnt,ncblk_trop,ncblk_op,
     &                   cinfo_cntc,cinfo_tropc,lstrcnt,lstrtroptmp,
     &                   cinfo_cntc(1,2),cinfo_tropc(1,2),
     &                                   cinfo_opc(1,2),
     &                   idxmsc_dis_a,idxmstropdis_c,
     &                   gmc_dis_a,gmtropdis_c,map_info_c,
     &                   strmap_info,nsym,str_info%ngraph)

              lenmap = get_lenmap(lstrcnt(ncblk_cnt+1),
     &                                  lstrtroptmp(ncblk_trop+1),
     &                                  map_info_a,nablk_op)
              call get_strmap_blk_c(map_tropcnta,
     &                   nablk_cnt,nablk_trop,nablk_op,
     &                   cinfo_cnta,cinfo_tropa,
     &                     lstrcnt(ncblk_cnt+1),
     &                             lstrtroptmp(ncblk_trop+1),
     &                   cinfo_cnta(1,2),cinfo_tropa(1,2),
     &                                   cinfo_opa(1,2),
     &                   idxmsc_dis_a,idxmstropdis_a,
     &                   gmc_dis_a,gmtropdis_a,map_info_a,
     &                   strmap_info,nsym,str_info%ngraph)

              ! calculate the trace for this block
              call trace_op_blk(xfac*casign,
     &                   xtropblk,xop(idxop),
     &                   tra_op, tra_trop,
     &                   ncblk_op,nablk_op,ncblk_troptmp,nablk_troptmp,
     &                   ncblk_cnt,nablk_cnt,
     &                   cinfo_opc(1,3),cinfo_opa(1,3),
     &                   cinfo_troptmpc(1,3),cinfo_troptmpa(1,3),
     &                   lstrop,lstrtroptmp,
     &                   lstrcnt,
     &                   map_info_c, map_info_a,
     &                   map_tropcntc, map_tropcnta
     &                   )
c dbg
c              print *,'after blk: ',xtropblk(1),xtropblk(1)-xold
c              xold = xtropblk(1)
c dbg
            end do cac_loop

            ! if necessary, reorder trop block:
            if (reo_trop.and.nonzero) then
              call reo_blk_wmaps_c(1d0,xtrop,xtropblk,
     &                  lentrop,lblk_troptmp,  ! just for checks
     &                  reo_info%sign_reo,
     &                  tra_trop,
     &                  ms_op_tr_c(2),ms_op_tr_a(2),
     &                                  igam_op_tr_c(2),igam_op_tr_a(2),
     &                  mstropdis_c,mstropdis_a,gmtropdis_c,gmtropdis_a,
     &                  ncblk_troptmp,nablk_troptmp,
     &                  cinfo_troptmpc,cinfo_troptmpa,
     &                  lstrtroptmp,
     &                  me_trop,iblktrop,
     &                  ncblk_trop,nablk_trop,
     &                  cinfo_tropc,cinfo_tropa,
     &                  reo_info%ncblk_reo,reo_info%nablk_reo,
     &                  reo_info%cinfo_reo_c,reo_info%cinfo_reo_a,
     &                  reo_info%ncblk_reo0,reo_info%nablk_reo0,
     &                  reo_info%cinfo_opreo0c,reo_info%cinfo_opreo0a,
     &                  reo_info%map_reo1c,reo_info%map_reo1a,
     &                  reo_info%map_reo2c,reo_info%map_reo2a,
     &                  nsym,str_info,strmap_info)
            end if

          end do caex_loop

        end do gam_loop
      end do ms_loop

      if (ntest.ge.1000) then
        if (iblktrop.gt.0
     &       ) then
          write(lulog,*) 'operator(out) on exit (',trim(trop%name),')'
          call wrt_mel_buf(lulog,5,xtrop,me_trop,
     &         iblktrop,iblktrop,str_info,orb_info)
        end if
      end if
c dbg
      if (trim(trop%name).eq.'Z-INT' .and. iblktrop.eq.2)
     &     print *,'jjjj ',xtrop(1),xtrop(4)
c dbg

      if (type_xret.eq.2) then
        xret(1) = xtrop(1)
      else if (type_xret.eq.1) then
        xret(1) = ddot(lentrop,xtrop,1,xtrop,1)
      end if

      ! put result to disc
      if (.not.buftrop) then
        call put_vec(fftrop,xtrop,idofftrop+idxst_trop,
     &                    idofftrop+idxst_trop-1+lentrop)
      end if

      deallocate(
     &     gmopdis_c, gmopdis_a,
     &     gmtropdis_c, gmtropdis_a,
     &     gmc_dis_c , gmc_dis_a ,
     &     msopdis_c, msopdis_a,
     &     mstropdis_c, mstropdis_a,
     &     msc_dis_c , msc_dis_a ,
     &     idxmsopdis_c, idxmsopdis_a,
     &     idxmstropdis_c, idxmstropdis_a,
     &     idxmsc_dis_c , idxmsc_dis_a ,
     &     lstrcnt,
     &     lstrop,lstrtroptmp,
     &     dmap_troptmpc,dmap_troptmpa,
     &     dmap_opc,dmap_opa
     &     )

      ifree = mem_flushmark()

      if (ntest.ge.100) then
        if (type_xret.ne.0)
     &       write(lulog,*) 'xret on exit = ',xret(1)
      end if

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
! threading, single thread version
module thread
implicit none
contains

subroutine int_thread(nthread)
integer, intent(in) :: nthread
if (nthread > 1) stop 'Multithread not enabled'
end subroutine

end module

max_stars_repo_name: gely/coseis

max_issues_repo_name: gely/coseis

max_forks_repo_name: gely/coseis


text:
      SUBROUTINE DSBMV( UPLO, N, K, ALPHA, A, LDA, X, INCX,
     $                  BETA, Y, INCY )
*
*  -- Automatically Tuned Linear Algebra Software (ATLAS)
*     (C) Copyright 2000 All Rights Reserved
*
*  -- ATLAS routine -- F77 Interface -- Version 3.9.24 -- December 25, 2000
*
*  Author         : Antoine P. Petitet
*  Originally developed at the University of Tennessee,
*  Innovative Computing Laboratory,  Knoxville TN, 37996-1301, USA.
*
*  ---------------------------------------------------------------------
*
*  -- Copyright notice and Licensing terms:
*
*  Redistribution  and  use in  source and binary forms, with or without
*  modification, are  permitted provided  that the following  conditions
*  are met:
*
*  1. Redistributions  of  source  code  must retain the above copyright
*     notice, this list of conditions and the following disclaimer.
*  2. Redistributions in binary form must reproduce  the above copyright
*     notice,  this list of conditions, and the  following disclaimer in
*     the documentation and/or other materials provided with the distri-
*     bution.
*  3. The name of the University,  the ATLAS group,  or the names of its
*     contributors  may not be used to endorse or promote products deri-
*     ved from this software without specific written permission.
*
*  -- Disclaimer:
*
*  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 UNIVERSITY
*  OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,  INDIRECT, INCIDENTAL, SPE-
*  CIAL, 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 THEO-
*  RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT  (IN-
*  CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
*  THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*  ---------------------------------------------------------------------
*
*     .. Scalar Arguments ..
      CHARACTER*1        UPLO
      INTEGER            INCX, INCY, K, LDA, N
      DOUBLE PRECISION   ALPHA, BETA
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
*     ..
*
*  Purpose
*  =======
*
*  DSBMV performs the matrix-vector operation
*
*     y := alpha*A*x + beta*y,
*
*  where alpha and beta are scalars, x and y are n-element vectors and A
*  is an n by n symmetric band matrix, with k super-diagonals.
*
*  Arguments
*  =========
*
*  UPLO    (input)                       CHARACTER*1
*          On entry, UPLO  specifies whether the upper or lower triangu-
*          lar part of the band matrix A is being supplied as follows:
*
*              UPLO = 'U' or 'u'   The upper triangular part of A is
*                                  being supplied.
*
*              UPLO = 'L' or 'l'   The lower triangular part of A is
*                                  being supplied.
*
*          Unchanged on exit.
*
*  N       (input)                       INTEGER
*          On entry, N specifies the order of the matrix A. N must be at
*          least zero. Unchanged on exit.
*
*  K       (input)                       INTEGER
*          On entry, K  specifies  the number of  super-diagonals of the
*          matrix A. K must satisfy 0 .le. K. Unchanged on exit.
*
*  ALPHA   (input)                       DOUBLE PRECISION
*          On entry, ALPHA specifies the scalar alpha.   When  ALPHA  is
*          supplied as zero then  A and X  need not be set on input. Un-
*          changed on exit.
*
*  A       (input)                       DOUBLE PRECISION array
*          On entry, A is an array of DIMENSION ( LDA, n ). Before entry
*          with UPLO = 'U' or 'u',  the leading (k + 1) by n part of the
*          array  A  must  contain the upper triangular band part of the
*          symmetric matrix, supplied column by column, with the leading
*          diagonal of the matrix in row (k + 1) of the array, the first
*          super-diagonal  starting at position 2 in row k,  and  so on.
*          The  top left k by k  triangle of the array A is not referen-
*          ced.  The following program  segment  will transfer the upper
*          triangular part of a symmetric band  matrix from conventional
*          full matrix storage to band storage:
*
*                DO 20, J = 1, N
*                   M = K + 1 - J
*                   DO 10, I = MAX( 1, J - K ), J
*                      A( M + I, J ) = matrix( I, J )
*             10    CONTINUE
*             20 CONTINUE
*
*          Before entry with  UPLO = 'L' or 'l', the leading (k+1)  by n
*          part of  the array A must contain the lower  triangular  band
*          part of the symmetric matrix, supplied column by column, with
*          the leading diagonal of the matrix in row 1 of the array, the
*          first sub-diagonal starting at position  1  in row  2, and so
*          on.  The bottom right  k by k  triangle of the array A is not
*          referenced.  The  following program segment will transfer the
*          lower triangular part of a symmetric band matrix from conven-
*          tional full matrix storage to band storage:
*
*                DO 20, J = 1, N
*                   M = 1 - J
*                   DO 10, I = J, MIN( N, J + K )
*                      A( M + I, J ) = matrix( I, J )
*             10    CONTINUE
*             20 CONTINUE
*
*          Unchanged on exit.
*
*  LDA     (input)                       INTEGER
*          On entry, LDA  specifies the first dimension of A as declared
*          in the calling (sub) program. LDA must be at least k + 1. Un-
*          changed on exit.
*
*  X       (input)                       DOUBLE PRECISION array
*          On entry,  X  is an incremented array of dimension  at  least
*          ( 1 + ( n - 1 ) * abs( INCX ) ).  Before entry, the incremen-
*          ted array X must contain the vector x. Unchanged on exit.
*
*  INCX    (input)                       INTEGER
*          On entry, INCX specifies the increment for the elements of X.
*          INCX must not be zero. Unchanged on exit.
*
*  BETA    (input)                       DOUBLE PRECISION
*          On entry,  BETA  specifies the scalar  beta.   When  BETA  is
*          supplied as zero then Y  need not be set on input.  Unchanged
*          on exit.
*
*  Y       (input/output)                DOUBLE PRECISION array
*          On entry,  Y  is an incremented array of dimension  at  least
*          ( 1 + ( n - 1 ) * abs( INCY ) ).  Before entry with BETA non-
*          zero, the incremented array  Y  must contain the vector y. On
*          exit, Y is overwritten by the updated vector y.
*
*  INCY    (input)                       INTEGER
*          On entry, INCY specifies the increment for the elements of Y.
*          INCY must not be zero. Unchanged on exit.
*
*  Further Details
*  ===============
*
*  For further information on the Level 1 BLAS specification, see:
*
*  ``A Proposal for Standard Linear Algebra Subprograms''  by R. Hanson,
*  F. Krogh and C. Lawson, ACM SIGNUM Newsl., 8(16), 1973,
*
*  ``Basic Linear Algebra Subprograms for Fortran Usage''  by C. Lawson,
*  R. Hanson, D. Kincaid and F. Krogh,  ACM Transactions on Mathematical
*  Software, 5(3) pp 308-323, 1979.
*
*  For further information on the Level 2 BLAS specification, see:
*
*  ``An  Extended Set of  FORTRAN  Basic Linear Algebra Subprograms'' by
*  J. Dongarra,  J. Du Croz,  S. Hammarling and R. Hanson,  ACM Transac-
*  tions on Mathematical Software, 14(1) pp 1-17, 1988.
*
*  ``Algorithm 656: An extended Set of Basic Linear Algebra Subprograms:
*  Model Implementation and Test Programs''  by J. Dongarra, J. Du Croz,
*  S. Hammarling and R. Hanson,  ACM  Transactions on Mathematical Soft-
*  ware, 14(1) pp 18-32, 1988.
*
*  For further information on the Level 3 BLAS specification, see:
*
*  ``A Set of Level 3 Basic Linear Algebra Subprograms'' by J. Dongarra,
*  J. Du Croz, I. Duff and S. Hammarling, ACM Transactions on Mathemati-
*  cal Software, 16(1), pp 1-17, 1990.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            ILOWER, IUPPER
      PARAMETER          ( IUPPER = 121, ILOWER = 122 )
*     ..
*     .. Local Scalars ..
      INTEGER            INFO, IUPLO
*     ..
*     .. External Subroutines ..
      EXTERNAL           ATL_F77WRAP_DSBMV, XERBLA
*     ..
*     .. External Functions ..
      EXTERNAL           LSAME
      LOGICAL            LSAME
*     ..
*     .. Executable Statements ..
*
      INFO = 0
*
      IF(      LSAME( UPLO  , 'U' ) ) THEN
         IUPLO = IUPPER
      ELSE IF( LSAME( UPLO  , 'L' ) ) THEN
         IUPLO = ILOWER
      ELSE IF( INFO.EQ.0 ) THEN
         INFO  = 1
      END IF
*
      IF( INFO.EQ.0 ) THEN
         IF( N.LT.0 ) THEN
            INFO = 2
         ELSE IF( K.LT.0 ) THEN
            INFO = 3
         ELSE IF( LDA.LT.( K + 1 ) ) THEN
            INFO = 6
         ELSE IF( INCX.EQ.0 ) THEN
            INFO = 8
         ELSE IF( INCY.EQ.0 ) THEN
            INFO = 11
         END IF
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DSBMV ', INFO )
         RETURN
      END IF
*
      CALL ATL_F77WRAP_DSBMV( IUPLO, N, K, ALPHA, A, LDA, X, INCX,
     $                        BETA, Y, INCY )
*
      RETURN
*
*     End of DSBMV
*
      END

ext: f

max_stars_repo_name: kevleyski/math-atlas

max_issues_repo_name: kevleyski/math-atlas

max_forks_repo_name: kevleyski/math-atlas


text:
      subroutine clst2dCF(ns,iflop,clst)
c ---------------------------------------------------------------------
c PROGRAM WRITTEN BY: Donald J. Jacobs                    Sept  7, 1996
c         This program takes the connection information of a 2d generic
c network as defined by the CF-bond linkage table (ie. linkcf() ) and
c the pebble index  pebble(,)  containing all rigidity information, and
c decomposes the network into its unique set of rigid clusters. It also
c calculates the number of floppy modes using the cluster statistics
c sum rule as an additional self-consistency check. Various moments of
c the cluster statistics is calculated as well.
c ---------------------------------------------------------------------
      include 'par_2dcfrp'
      integer s,so,s1,sf,stest,btag,rigid,smax,point(nsmax)
      integer linkcf(nb2max),clst(nb2max),tag,tag1
      integer multcf(nsmax),state(nb2max),stressed
      integer pivots(nsmax),mark,covered
      integer   stat(nsmax),gradus
      integer pebble(id,0:nmax),shell(0:nmax),block(0:nmax)
      integer btrack(0:nmax),len,lth

      data rigid /2147483644/  mark /32/
c          rigid = 2**31 - 4
c =====================================================================
      common/rigidity/pebble,block,tag,shell,nsfail
      common/network/n,linkcf,point,multcf,state,covered,stressed
      common/search/btrack,btag
c =====================================================================
      smax = 1
      nsvoid = 0
      tag1 = 0
      do s=1,ns
       stat(s) = 0
       pivots(s) = 0
      enddo
c --------------------------------------------- identify isolated sites
         do so=1,ns
          if( multcf(so) .lt. 1 ) stat(1) = stat(1) + 1
         enddo
c ----------------------------------------------------- void correction
      stat(1) = stat(1) - nsvoid
      nscorr = ns - nsvoid
c ----------------------------------------- decompose all clusters with
c                                                    more than one bond
         do so=1,ns
          ipo = point(so)
            do j=1,multcf(so)
               if( state(ipo+j) .lt. mark ) then
               sf = linkcf(ipo+j)
c --------------------------------------------------------- map cluster
               call free2d2(so)
               call free2d1(sf,so,np1)

c +++++++++++++++++++++++++++++++++++++++++++++++ REMOVE
                  if( np1 .eq. 2 ) then
                  write(86,6543)
                  write(*,*)'in decomp',so,sf
 6543             format(5x,'More than three pebbles collected ',
     &                  'across a single bond in mapping cluster')
                  stop
                  endif
                  if( np1 .ne. 1 ) then
                  write(86,6544)
 6544             format(5x,'Three pebbles can not be collected ',
     &                  'across a single bond in mapping cluster')
                  stop
                  endif
c +++++++++++++++++++++++++++++++++++++++++++++++

c ------------------------------------------------- initialize the burn
c                                      btag marks sites already checked
               btag = btag + 1
               btrack(so) = btag
               btrack(sf) = btag
c -------------------------------------------------- mark rigid cluster
               block(so) = rigid
               block(sf) = rigid
               shell(1) = so
               shell(2) = sf
               kmax = 2
c ----------------------------------------------- burn out rigid region
               k = 0
  100          k = k + 1
                  if( k .gt. kmax ) then
c ------------------------------------------------------ record cluster
                  tag1 = tag1+1
                  stat(kmax) = stat(kmax) + 1
                  if( kmax .gt. smax ) smax = kmax
                     do k=1,kmax
                     s = shell(k)
                     pivots(s) = pivots(s) + 1
                     ips = point(s)
                        do jj=1,multcf(s)
                        stest = linkcf(ips+jj)
                           if( block(stest) .eq. rigid ) then
                           state(ips+jj) = state(ips+jj) + mark
                           clst(ips+jj) = tag1
                           endif
                        enddo
                     enddo
c ------------------------------------------------------- re-initialize
                     do k=1,kmax
                     block( shell(k) ) = 0
                     enddo
                  else
c -------------------------------------------------------- grow cluster
                  s = shell(k)
                  ips = point(s)
c ------------------------------------------- search and expand about s
                     do jtest=1,multcf(s)
c ------------------------------------------------------------ unmarked
                        if( state(ips+jtest) .lt. mark ) then
                        stest = linkcf(ips+jtest)
c ----------------------------------------------------------- unchecked
                           if( btrack(stest) .lt. btag ) then
                           btrack(stest) = btag
                           call free2d0(stest,kmax)
                           endif
                        endif
                     enddo
                  goto 100
                  endif
               endif
            enddo
         enddo
c ------------------------------------ calculate number of floppy modes
       ic = 0
          do s=2,smax
          ic = ic + stat(s)* (2*s - 3)
          enddo
       iflop2 = id* nscorr - ic
       npivot = 0
         do s=1,ns
         if( pivots(s) .gt. 1 ) npivot = npivot+1
         enddo
c --------------------------------------------- self consistency checks
       iflag = 1
       itemp = 0
          do s=1,ns
             if( pivots(s) .gt. mcoord ) then
                if( iflag .eq. 1 ) then
                write(87,6723)
 6723           format(5x,'PIVOT # exceeds maximum coordination')
                endif
             iflag = 0
             endif
          if( pivots(s) .eq. 0 ) itemp = itemp +  1
          enddo
       itemp2 = stat(1) + nsvoid
          if( itemp .ne. itemp2 ) then
          iflag = 0
          write(87,6724)
 6724     format(5x,'PIVOT(s) = 0 is inconsistent with ',
     &              ' stat(1) + nsvoid')
          write(87,*) itemp,stat(1),nsvoid
          endif
       if( iflop2 .ne. iflop ) then
       write(87,6725) iflop,iflop2
 6725  format(5x,'ERROR: F_direct = ',i8,' .NE. F_sumrule = ',i8)
       iflag = 0
       endif
!       if( iflag .ne. 1 ) stop

c ------------------------------------------------------ restore state
         do s=1,ns
          ips = point(s)
            do j=1,multcf(s)
            if( state(ips+j) .ge. mark ) state(ips+j) =
     &         state(ips+j) - mark
            enddo
         enddo

      return
      end

max_stars_repo_name: VardaHagh/pebbleGame

max_issues_repo_name: VardaHagh/pebbleGame

max_forks_repo_name: VardaHagh/pebbleGame


text:
c
c     ===================================================================
      function alpha2(rs,dz,sp,iexch,exchg)
c     ===================================================================
c
      implicit none
c
      integer  iexch
      integer  incof
      integer  i
c
      real*8   alpha2
      real*8   rs
      real*8   dz
      real*8   sp
      real*8   exchg
      real*8   ccp
      real*8   rp
      real*8   ccf
      real*8   rf
      real*8   a(3)
      real*8   b(3)
      real*8   c(3)
      real*8   x0(3)
      real*8   cst
      real*8   aip
      real*8   fnot
      real*8   bip
      real*8   for3
      real*8   thrd
      real*8   vxx(2)
      real*8   vxcs(2)
      real*8   g(3)
      real*8   dg(3)
      real*8   tbq(3)
      real*8   tbxq(3)
      real*8   bxx(3)
      real*8   q(3)
      real*8   bb(3)
      real*8   cx(3)
      real*8   fm
      real*8   fdz
      real*8   ex
      real*8   exf
      real*8   xp
      real*8   xf
      real*8   gp
      real*8   gf
      real*8   exc
      real*8   excf
      real*8   dedz
      real*8   gpp
      real*8   gfp
      real*8   depd
      real*8   defd
      real*8   decd
      real*8   bfc
      real*8   zp1
      real*8   zm1
      real*8   xr
      real*8   pex
      real*8   xrsq
      real*8   qi
      real*8   txb
      real*8   fx
      real*8   arct
      real*8   dxs
      real*8   vcc
      real*8   facc
      real*8   ecp
      real*8   zp3
      real*8   zm3
      real*8   zp3m3
      real*8   fx1
      real*8   z4
      real*8   fz
      real*8   beta
      real*8   ec
      real*8   f3ex
c
c     data for von barth-hedin
c
      data     ccp,rp,ccf,rf/0.0450d+00,21.0d+00,
     >                       0.02250d+00,52.9166820d+00/
c
c     data for vosko-wilk-nusair
c
      data     incof/0/
      data     a/-0.0337740d+00,0.06218140d+00,0.03109070d+00/
      data     b/1.131070d+00,3.727440d+00,7.060420d+00/
      data     c/13.00450d+00,12.93520d+00,18.05780d+00/
      data     x0/-0.00475840d+00,-0.104980d+00,-0.32500d+00/
      data     cst,aip/1.923661050d+00,0.916330590d+00/
      data     fnot,bip,for3,thrd/1.709920950d+00,0.259921050d+00,
     >                            1.33333333330d+00,0.33333333330d+00/
c
c     the following are constants needed to obtain potential
c     which are in g.s. painter's paper
c     =====given here for check=====(generated below)
c
c     data q/.7123108918d+01,.6151990820d+01,.4730926910d+01/
c     data bxx/-.4140337943d-03,-.3116760868d-01,-.1446006102/
c     data tbq/.3175776232d+00,.1211783343d+01,.2984793524d+01/
c     data tbxq/.3149055315d+00,.1143525764d+01,.2710005934d+01/
c     data bb/.4526137444d+01,.1534828576d+02,.3194948257d+02/
c
c     write(6,'('' iexch,rs,dz,sp '',i5,3f10.5)') iexch,rs,dz,sp
c
c*********** for positron potential construction iexch.ge.100 *********
      if(iexch.ge.100)then
        alpha2=0.d0
        exchg=0.d0
        return
      endif
c**********************************************************************
      go to (10,20) iexch
c
c     von barth-hedin  exch-corr potential
c     j. phys. c5,1629(1972)
c
c
  10  continue
      fm=2.0d+00**(4.0d+00/3.0d+00)-2.0d+00
      fdz = ((1.0d+00+dz)**(4.0d+00/3.0d+00)
     >     +(1.0d+00-dz)**(4.0d+00/3.0d+00)-2.0d+00)/fm
      ex=-0.916330d+00/rs
      exf=ex*2.0d+00**0.333333330d+00
      xp=rs/rp
      xf=rs/rf
      gp = (1.0d+00+xp**3)*log(1.0d+00+1.0d+00/xp)
     >    -xp*xp +xp/2.0d+00 - 0.333333330d+00
      gf = (1.0d+00+xf**3)*log(1.0d+00+1.0d+00/xf)
     >    -xf*xf +xf/2.0d+00 - 0.333333330d+00
      exc = ex-ccp*gp
      excf=exf-ccf*gf
      dedz= (4.0d+00/3.0d+00)*(excf-exc)
     >     *((1.0d+00+dz)**(1.0d+00/3.0d+00)
     >     -(1.0d+00-dz)**(1.0d+00/3.0d+00))/fm
      gpp = 3.0d+00*xp*xp*log(1.0d+00+1.0d+00/xp)-1.0d+00/xp
     >     +1.50d+00-3.0d+00*xp
      gfp = 3.0d+00*xf*xf*log(1.0d+00+1.0d+00/xf)-1.0d+00/xf
     >     +1.50d+00-3.0d+00*xf
      depd=-ex/rs-ccp/rp*gpp
      defd=-exf/rs-ccf/rf*gfp
      decd=depd+(defd-depd)*fdz
c     exchange-correlation energy
      exchg= exc + (excf-exc)*fdz
c     exchange-correlation potential
      alpha2 = exc+(excf-exc)*fdz-rs*decd/3.0d+00
     >        +sp*(1.0d+00-sp*dz)*dedz
      return
c
c
  20  continue
      if(incof.ne.0) go to 30
      incof=2
c
c     vosko-wilk-nusair exch-corr potential
c     taken from g.s. painter
c     phys. rev. b24 4264,1981
c
c     generate constant coefficients for the parameterization (v-w-n)
c
      do i=1,3
         cx(i)= x0(i)**2 + b(i)*x0(i) + c(i)
         bfc= 4.0d+00*c(i) - b(i)**2.0d+00
         q(i)= sqrt(bfc)
         bxx(i)= b(i)*x0(i)/cx(i)
         tbq(i)= 2.0d+00*b(i)/q(i)
         tbxq(i)= tbq(i) + 4.0d+00*x0(i)/q(i)
         bb(i)= 4.0d+00*b(i)*( 1 - x0(i)*(b(i) + 2.0d+00*x0(i))/cx(i) )
      enddo
c
  30  continue
      zp1= 1.0d+00 + dz
      zm1= 1.0d+00 - dz
      xr=sqrt(rs)
      pex= -aip/rs
      xrsq= rs
c
c     generate g(i)=alpha,epsilon fct.s
c     and their derivatives dg(i)
c     1=alpha(spin stiffness)  2=ecp  3=ecf
c
      do i=1,3
         qi=q(i)
         txb= 2.0d+00*xr + b(i)
         fx= xrsq + xr*b(i) + c(i)
         arct= atan2(qi,txb)
         dxs= (xr-x0(i))**2/fx
         g(i)=a(i)*( log(xrsq/fx) + tbq(i)*arct-bxx(i)*(log(dxs)
     >              +tbxq(i)*arct) )
         dg(i)=a(i)*( 2.0d+00/xr - txb/fx
     >               -bxx(i)*(2.0d+00/(xr-x0(i))-txb/fx)
     >               -bb(i)/(qi**2 + txb**2) )
      enddo
c
      ecp=g(2)
      zp3=zp1**thrd
      zm3=zm1**thrd
      zp3m3=zp3-zm3
c     part of last term in vx   eq(13)
      fx1=.50d+00*for3*pex*zp3m3
      z4= dz**4
      fz= cst*(zp1**for3 + zm1**for3 - 2.0d+00)
      beta= fnot*( g(3)-g(2) )/g(1) -1.0d+00
      ec= ecp + fz*g(1)*( 1.0d+00 + z4*beta )/fnot
      ex= pex*( 1.0d+00 + fz*bip )
      f3ex= for3*ex
c     echange-correlation energy
      exchg= ec + ex
c     exchange potential
      vxx(1)= f3ex + fx1*zm1
      vxx(2)= f3ex - fx1*zp1
c     correlation potential
      vcc= ec - xr*( (1.0d+00-z4*fz)*dg(2) + z4*fz*dg(3)
     >              +(1.0d+00 - z4)*fz*dg(1)/fnot )/6.0d+00
c
      facc= 4.0d+00*g(1)*( dz**3*fz*beta
     >                    +( 1.0d+00 + beta*z4 )*zp3m3/(6.0d+00*bip) )
     >                  /fnot
c
c     exch-corr. potential for each spin as called in newpot
c
      vxcs(1)= vcc + zm1*facc + vxx(1)
      vxcs(2)= vcc - zp1*facc + vxx(2)
c
      if( sp.ge.0 ) alpha2= vxcs(1)
      if( sp.lt.0 ) alpha2= vxcs(2)
      return
      end

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: hkershaw-brown/MuST