text:       subroutine invsqrt(mel_inp,mel_inv,nocc_cls,half,
       &     get_u,mel_u,pass_spc,mel_spc,lmodspc,
       &     op_info,orb_info,str_info,strmap_info)
  *----------------------------------------------------------------------*
  *     Routine to calculate S^(-0.5) of density matrices.
  *     The ME list is split into matrices that can either contain
  *     elements from a single distribution or from several distributions.
  *     half = true: only makes a half transform U*s^(-0.5)
  *     half = false: also does half transform and in addition
  *                   returns projector matrix U*1s*U^+ on input list
  *     get_u= true: return unitary matrix U on mel_u
  *     lmodspc=true: special mode: no sequential projections; symmetrize;
  *                   just diagonal blocks; no spin splitting for ms=0;
  *                   return diagonalized list on inp and unitary on inv
  *
  *     works also for sums of density matrices of the structure:
  *       /0 0 0 0\
  *       \0 0 x 0/  (not necessarily igastp=3)
  *       /0 0 y 0\
  *       \0 0 y 0/
  *       /0 0 x 0\
  *       \0 0 0 0/
  *
  *     IMPORTANT NOTE: At present the sequence matters (for project!=1):
  *        e.g.  (,V;V,V;V,) then (,V;V,;,) then (,;,V;V,) and then (,;,;,)
  *
  *     matthias, dec 2009 (adopted from invert.f)
  *----------------------------------------------------------------------*
  
        implicit none
  
        include 'stdunit.h'
        include 'opdim.h'
        include 'def_orbinf.h'
        include 'mdef_operator_info.h'
        include 'def_graph.h'
        include 'def_strinf.h'
        include 'def_strmapinf.h'
        include 'ifc_memman.h'
        include 'hpvxseq.h'
        include 'multd2h.h'
        include 'ifc_input.h'
        include 'routes.h'
        include 'molpro_out.h'
  
        type gam_generator_t
            integer::last_gam
            integer::ngam
            integer::total_gam
        end type
  
        type ms_generator_t
          integer::last_ms
          integer::nel
          integer::total_ms
        end type
  
        integer, parameter ::
       &     ntest = 00
        character(len=*),parameter::
       &     i_am = "invsqrt"
  
        type(orbinf), intent(in) ::
       &     orb_info
        type(operator_info), intent(in) ::
       &     op_info
        type(strinf), intent(in) ::
       &     str_info
        type(strmapinf), intent(in) ::
       &     strmap_info
        type(me_list), intent(in) ::
       &     mel_inp, mel_inv, mel_u, mel_spc
        integer, intent(in) ::
       &     nocc_cls
        logical, intent(in) ::
       &     half, get_u, pass_spc, lmodspc
  
        integer, parameter ::
       &     maxrank = 5
        type(graph), pointer ::
       &     graphs(:)
  
        logical ::
       &     bufin, bufout, first, ms_fix, fix_success, onedis, transp,
       &     logdum, sgrm, bufu, svdonly, reg_tik, first2,lcycle
        logical, pointer ::
       &     blk_used(:), blk_redundant(:)
  c      logical ::
  c     &     loop(nocc_cls)
        integer ::
       &     ifree, nbuff, idxmsa, iocc_cls, iexc_cls,
       &     msmax, msa, igama, idx, jdx, ngam,
       &     ioff, njoined,
       &     idxdis, lenca, iblkoff, ncblk, nablk,
       &     msc, igamc, idxmsc, mscmax,
       &     ndis, ndim, off_col, off_line,
       &     idxa1, idxa2, idxc1, idxc2, icol, iline,
       &     msmax_sub, ms1, ms2, msc1, msc2, msa1, msa2,
       &     gamc1, gamc2, gama1, gama2, igam, na1, na2, nc1, nc2, ij,
       &     maxbuf, ngraph, ioff2, idxmsa2, ndis2, idxdis2, idx2, nel,
       &     nsing, ising, itrip, ntrip, icnt_sv, icnt_sv0,
       &     nalph, nbeta, nc1mx, na1mx, jocc_cls, jblkoff,
       &     ncblk2, nablk2, len2(4), off_line2, off_col2, off_colmax,
  c dbg
  c     &     ipass,
  c dbgend
       &     off_linmax, maxbuf_tmp,
       &     rankdim(maxrank), rankoff(maxrank), nrank,
       &     irank, jrank, idxst, idxnd, rdim, idxst2, idxnd2, rdim2,
       &     icnt_cur, ih, ip, iexc, tocc_cls, gno, project,
       &     krank, idxst3, idxnd3, rdim3, curr_rec, len_rec,
       &     iprint
        real(8) ::
       &     fac, xmax, xmin, xdum, omega2
        real(8), pointer ::
       &     buffer_in(:), buffer_out(:), scratch(:,:), scratch2(:,:),
       &     sing(:,:), trip(:,:), sing2(:,:), trip2(:,:),
       &     palph(:), pbeta(:), proj(:,:), norm(:), scratch3(:,:),
       &     scratch4(:,:), sing3(:,:), trip3(:,:), buffer_u(:),
       &     svs(:),
       &     scratch_tmp1(:,:),scratch_tmp2(:,:),scratch_tmp3(:,:),
       &     scratch_tmp4(:,:)
        real(8), target ::
       &     xdummy(1,1)
  c dbg
  c      integer, pointer ::
  c     &     matrix(:,:)
  c dbgend
  
        integer, pointer ::
       &     hpvx_csub(:),hpvx_asub(:),
       &     occ_csub(:), occ_asub(:),
       &     graph_csub(:), graph_asub(:),
       &     msdis_c(:),  msdis_a(:),
       &     idxmsdis_c(:),  idxmsdis_a(:),
       &     msdis_c2(:),  msdis_a2(:),
       &     idxmsdis_c2(:),  idxmsdis_a2(:),
       &     gamdis_c(:), gamdis_a(:),
       &     len_str(:),
       &     hpvx_occ(:,:,:), ca_occ(:,:), idx_graph(:,:,:),
       &     ldim_opin_c(:), ldim_opin_a(:),
       &     istr_csub(:), istr_asub(:),
       &     istr_csub_flip(:), istr_asub_flip(:),
       &     flipmap_c(:), flipmap_a(:),
       &     flmap(:,:),
       &     iocc(:,:,:), idx_g(:,:,:), msdst(:,:,:), igamdst(:,:,:),
       &     idorb(:), idspn(:), idspc(:), lexlscr(:,:),
       &     hpvx_csub2(:),hpvx_asub2(:),
       &     occ_csub2(:), occ_asub2(:),
       &     graph_csub2(:), graph_asub2(:),
       &     iocc2(:,:,:),idx_g2(:,:,:),
       &     bins(:,:), ex2occ_cls(:)
  
        type(filinf), pointer ::
       &     ffinp, ffinv, ffu
        type(operator), pointer ::
       &     op_inv, op_inp, op_u, op_t
  
        integer, external ::
       &     ielprd, idx_msgmdst2, idx_str_blk3, msa2idxms4op, idxcount,
       &     idxlist, idx_oplist2
        logical, external ::
       &     next_tupel_ca, occ_is_diag_blk
        type(gam_generator_t)::
       &     gam_gen1,gam_gen2
        type(ms_generator_t)::
       &     ms_gen1,ms_gen2
  
        iprint=max(iprlvl,ntest)
  
        if (iprint.ge.100) write(lulog,*) '==== entered invsqrt ===='
  
  c dbg
  c      ipass = 0
  c dbgend
        call get_argument_value('method.MR','project',ival=project)
        call get_argument_value('method.MR','GNO',ival=gno)
        sgrm = project.eq.1.and.gno.eq.0.or.lmodspc ! only diagonal blocks
        if (lmodspc.and.(get_u.or..not.half))
       &   call quit(1,'invsqrt','lmodspc only with get_u=F, half=T')
  
        reg_tik = tikhonov.ne.0d0
        omega2 = tikhonov**2
  
        ffinp => mel_inp%fhand
        ffinv => mel_inv%fhand
        ffu => mel_u%fhand
        op_inp => mel_inp%op
        op_inv => mel_inv%op
        op_u => mel_u%op
        njoined = op_inp%njoined
  
        ngraph = str_info%ngraph
        hpvx_occ => op_inp%ihpvca_occ
        idx_graph => mel_inp%idx_graph
        ca_occ => op_inp%ica_occ
        graphs => str_info%g
  
  
        ms_fix = fix_ms_h(mel_inp, mel_inv, mel_u)
  
        ! Check whether files are buffered.
        bufin = .false.
        if(ffinp%buffered) bufin = .true.
        bufout = .false.
        if(ffinv%buffered) bufout = .true.
        bufu = .false.
        if (get_u.and.ffu%buffered) bufu = .true.
        if (ntest.ge.100) then
          write(lulog,*) 'in:  ',trim(mel_inp%label),' ',bufin
          write(lulog,*) 'out: ',trim(mel_inv%label),' ',bufout
          if (get_u)
       &    write(lulog,*) 'U:   ',trim(mel_u%label),' ',bufu
        end if
  
        ifree = mem_setmark('invsqrt')
  
        ! Number of irreps in symmetry group.
        ngam = orb_info%nsym
  
        ! Allocations made to maximum block length to save time.
        if(.not.bufin)then
          curr_rec=ffinp%current_record
          len_rec =ffinp%length_of_record
          nbuff = 0
          do iocc_cls = 1, nocc_cls
            if(op_inp%formal_blk(iocc_cls))
       &         cycle
  
            nbuff = nbuff+mel_inp%len_op_occ(iocc_cls)
          enddo
          ifree = mem_alloc_real(buffer_in,nbuff,'buffer_in')
          call get_vec(ffinp,buffer_in,(curr_rec-1)*len_rec+1,
       &       (curr_rec-1)*len_rec+nbuff)
        else
          if(iprint.ge.100)
       &       write(lulog,*)'Invert: input incore'
          buffer_in => ffinp%buffer(1:)
        endif
  
        if(.not.bufout)then
          nbuff = 0
          do iocc_cls = 1, nocc_cls
            nbuff = nbuff + mel_inv%len_op_occ(iocc_cls)
          enddo
          ifree= mem_alloc_real(buffer_out,nbuff,'buffer_out')
          buffer_out(1:nbuff) = 0d0
        else
          if(iprint.ge.100)
       &       write(lulog,*)'Invert: output incore'
          buffer_out => ffinv%buffer(1:)
        endif
  
        if(get_u.and..not.bufu)then
          nbuff = 0
          do iocc_cls = 1, nocc_cls
            nbuff = nbuff + mel_u%len_op_occ(iocc_cls)
          enddo
          ifree= mem_alloc_real(buffer_u,nbuff,'buffer_u')
          buffer_u(1:nbuff) = 0d0
        else if (get_u) then
          if(iprint.ge.100)
       &       write(lulog,*)'Invert: output (2) incore'
          buffer_u => ffu%buffer(1:)
        endif
  
  
        icnt_sv  = 0 ! we will count the
        icnt_sv0 = 0 ! number of singular values below threshold
        xmax = 0d0   ! largest excluded singular value
        xmin = 1234567890d0   ! smallest included singular value
        allocate(blk_used(nocc_cls),blk_redundant(nocc_cls),
       &         bins(17,nocc_cls),ex2occ_cls(nocc_cls))
        bins = 0 ! binning for singular values:
                 ! >10E0,>10E-1,...,>10E-15,0
        blk_used(1:nocc_cls) = .false.
        blk_redundant(1:nocc_cls) = .true.
        iexc_cls = 0
        tocc_cls = 0
  
        op_t => get_cluster_op_h(op_info) ! 'T' 'C' or null
  
        if (.not.half.and. (iprint.ge.3) ) write(lulog,*)
       &     'Input list will be overwritten by projector.'
  
        ! Loop over occupation class.
        if (ntest.ge.100) write(lulog,*) 'entering iocc_loop:'
        iocc_loop: do iocc_cls = 1, nocc_cls !loops over occupations of input operator
          iblkoff = (iocc_cls-1)*njoined
          if (ntest.ge.100) then
             write(lulog,'(1x,72("="))')
             write(lulog,*) 'iocc_cls = ',iocc_cls
             write(lulog,*) 'formal:   ',op_inp%formal_blk(iocc_cls)
             write(lulog,*) 'blk_used: ',blk_used(iocc_cls)
             write(lulog,*) 'len_op_occ: ', mel_inp%len_op_occ(iocc_cls)
            call wrt_occ_n(lulog,op_inp%ihpvca_occ(1,1,iblkoff+1),njoined)
          end if
  
          if (occ_is_diag_blk(hpvx_occ(1,1,iblkoff+1),njoined))
       &     tocc_cls = tocc_cls + 1 ! increment occ.cls of corresp. Op.
          if(op_inp%formal_blk(iocc_cls)) cycle iocc_loop
          if (blk_used(iocc_cls)) cycle iocc_loop
  
          blk_used(iocc_cls) = .true.
          iexc_cls = iexc_cls + 1
          ex2occ_cls(iexc_cls) = tocc_cls
  
          if (iprint.ge.10) write(lulog,*) 'current occ_cls: ',iocc_cls
  
          ! only one element? easy!
          ! (also regularization is never needed in this case)
          if (mel_inp%len_op_occ(iocc_cls).eq.1) then
            icnt_cur = icnt_sv - icnt_sv0
            ioff = mel_inp%off_op_gmo(iocc_cls)%gam_ms(1,1)
            buffer_out(ioff+1) = buffer_in(ioff+1)
            if (lmodspc) then ! only special mode
              call mat_svd_traf(1,buffer_out(ioff+1),buffer_in(ioff+1),
       &                       icnt_sv,icnt_sv0,xmax,xmin,
       &                       bins(1,iexc_cls))
            else if (get_u) then
              call invsqrt_mat(1,buffer_out(ioff+1),buffer_in(ioff+1),
       &                         half,buffer_u(ioff+1),get_u,
       &                         xdum,icnt_sv,icnt_sv0,xmax,xmin,
       &                         bins(1,iexc_cls))
  
            else
              call invsqrt_mat(1,buffer_out(ioff+1),buffer_in(ioff+1),
       &                       half,xdummy,get_u, !buffer_u: dummy
       &                       xdum,icnt_sv,icnt_sv0,xmax,xmin,
       &                       bins(1,iexc_cls))
            end if
            if (sgrm.and.icnt_cur.lt.icnt_sv-icnt_sv0)
       &       blk_redundant(iocc_cls) = .false.
            cycle iocc_loop
          end if
  
  
          ifree = mem_setmark('invsqrt_blk')
  
  !-----------
          ! for simplicity we distinguish the following two types:
          !  /0 0 0 0\       /0 0 0 0\  i.e. occupation classes  with
          !  \0 0 x 0/       \0 0 x 0/  several distributions per
          !  /0 0 y 0\  and  /0 0 0 0\  per Ms(A)/GAMMA(A) block
          !  \0 0 y 0/       \0 0 0 0/  and those with only single
          !  /0 0 x 0\       /0 0 x 0\  distributions
          !  \0 0 0 0/       \0 0 0 0/
          call get_num_subblk(ncblk,nablk,
       &       hpvx_occ(1,1,iblkoff+1),njoined)
          ! for simplicity we distinguish the following two types:
          !  /0 0 0 0\       /0 0 0 0\  i.e. occupation classes  with
          !  \0 0 x 0/       \0 0 x 0/  several distributions per
          !  /0 0 y 0\  and  /0 0 0 0\  per Ms(A)/GAMMA(A) block
          !  \0 0 y 0/       \0 0 0 0/  and those with only single
          !  /0 0 x 0\       /0 0 x 0\  distributions
          !  \0 0 0 0/       \0 0 0 0/
          if (ncblk.eq.1.and.nablk.eq.1) then
            onedis = .true.
          else if (ncblk.eq.2.and.nablk.eq.2) then
            onedis = .false.
          else
            call quit(1,'invsqrt','not adapted for this case')
          end if
  
  !----------
          allocate(hpvx_csub(ncblk),hpvx_asub(nablk),
       &           occ_csub(ncblk), occ_asub(nablk),
       &           graph_csub(ncblk), graph_asub(nablk),
       &           msdis_c(ncblk),  msdis_a(nablk),
       &           idxmsdis_c(ncblk),  idxmsdis_a(nablk),
       &           msdis_c2(ncblk),  msdis_a2(nablk),
       &           idxmsdis_c2(ncblk),  idxmsdis_a2(nablk),
       &           gamdis_c(ncblk), gamdis_a(nablk),
       &           len_str(ncblk+nablk),
       &           istr_csub(ncblk),istr_asub(nablk),
       &           istr_csub_flip(ncblk),istr_asub_flip(nablk),
       &           ldim_opin_c(ncblk),ldim_opin_a(nablk))
  
          ! set HPVX and OCC info
          call condense_occ(occ_csub, occ_asub,
       &                    hpvx_csub,hpvx_asub,
       &                    hpvx_occ(1,1,iblkoff+1),njoined,hpvxblkseq)
          ! do the same for the graph info
          call condense_occ(graph_csub, graph_asub,
       &                    hpvx_csub,hpvx_asub,
       &                    idx_graph(1,1,iblkoff+1),njoined,hpvxblkseq)
  
          ! set flip maps
          call strmap_man_flip(
       &       maxbuf,
       &       graph_csub,ncblk,
       &       str_info,strmap_info,orb_info)
  c dbg
  c        print *,'maxbuf flipmap_c:',maxbuf
          maxbuf = maxbuf + 5
  c dbgend
          ifree = mem_alloc_int(flipmap_c,maxbuf,'flipmap_c')
          call strmap_man_flip(
       &       maxbuf,
       &       graph_asub,nablk,
       &       str_info,strmap_info,orb_info)
  c dbg
  c        print *,'maxbuf flipmap_a:',maxbuf
          maxbuf = maxbuf + 5
  c dbgend
          ifree = mem_alloc_int(flipmap_a,maxbuf,'flipmap_a')
  
          ! simple case: only single distributions:
          if (onedis) then
            icnt_cur = icnt_sv - icnt_sv0
  
            ! we also need to distinguish:
            ! /0 0 0 0\     /0 0 0 0\
            ! \0 0 x 0/     \0 0 0 0/ i.e. if creators come first,
            ! /0 0 0 0\ and /0 0 x 0\ we need to transpose the
            ! \0 0 0 0/     \0 0 x 0/ scratch matrix
            ! /0 0 x 0\     /0 0 0 0\ (important if result is non-symmetric!)
            ! \0 0 0 0/     \0 0 0 0/
            transp = .false.
            do ij = 1, njoined
              if (sum(hpvx_occ(1:ngastp,1,iblkoff+ij)).gt.0) then
                transp = .true.
                exit
              else if (sum(hpvx_occ(1:ngastp,2,iblkoff+ij)).gt.0) then
                exit
              end if
            end do
            if (transp.and.iprint.ge.100) then
              write(lulog,*) 'Using transposed scratch matrix!'
            end if
  
            ! Loop over Ms of annihilator string.
            idxmsa = 0
            msmax = op_inp%ica_occ(1,iocc_cls)
            mscmax = op_inp%ica_occ(2,iocc_cls)
            if (msmax.ne.mscmax) call quit(1,'invsqrt',
       &         'need particle conserving operator')
            msa_loop : do msa = msmax, -msmax, -2
  
              idxmsa = idxmsa+1
              msc = msa + mel_inp%mst
              idxmsc = (msmax-msc)/2 + 1
  
              ! Loop over Irrep of annihilator string.
              igama_loop: do igama =1, ngam
  
                igamc = multd2h(igama,mel_inp%gamt)
                ndis = mel_get_ndis_h(mel_inp,iocc_cls,igamc,idxmsc)
  
                ndim = isqrt_h(
       &             mel_get_lenblock_h(mel_inv,iocc_cls,igama,idxmsa))
  
                if (ndim.gt.0.and.ndis.ne.1)
       &                call quit(1,i_am,'cannot handle this')
                if (ndim.eq.0) cycle igama_loop
  
                if (iprint.ge.15)
       &           write(lulog,'(a,3i8)') 'msa, gama, ndim:',msa,igama,
       &             ndim
  c     &           int(sqrt(dble(mel_inp%len_op_gmo(iocc_cls)%
  c     &                         gam_ms(igama,idxmsa))))
  
                if (iprint.ge.100)
       &           write(lulog,*) ' len = ',
       &             mel_inp%len_op_gmo(iocc_cls)%gam_ms(igama,idxmsa),
       &             ' ndis = ',ndis
  
                ioff = mel_inv%off_op_gmo(iocc_cls)%gam_ms(igama,idxmsa)
  
                call get_flipmap_blk(flipmap_c,
       &            ncblk,occ_csub,ndim,graph_csub,idxmsc,igamc,
       &            strmap_info,ngam,ngraph)
  
                ! single distribution can simply be read in as simple matrix
                allocate(scratch(ndim,ndim),svs(ndim))
                if (.not.half.or.lmodspc) then
                  allocate(scratch2(ndim,ndim))
                else
                  scratch2 => xdummy
                end if
                if (get_u) then
                  allocate(scratch3(ndim,ndim))
                else
                  scratch3 => xdummy
               end if
               call copy_buffer_1_h(buffer_in, scratch, !buffer_in -> scatch
       &            ndim,ndim,ioff,transp)
  
  
               if (msc.eq.0.and..not.lmodspc) then
                  ! here a splitting into "singlet" and "triplet" blocks is needed:
  
  c dbg
  c                write(lulog,*) 'flmap:'
  c                do icol = 1, ndim
  c                  write(lulog,'(i4,3i6)') icol,flipmap_c(icol)
  c                end do
  c dbgend
                  nsing = ndim
                  do icol = 1, ndim
                    if (abs(flipmap_c(icol)).eq.icol) nsing = nsing
       &                                  + sign(1,flipmap_c(icol))
                  end do
                  nsing = nsing/2
                  ntrip = ndim - nsing
  c dbg
  c                print *,'nsing: ',nsing
  c dbgend
  
                  ! do the pre-diagonalization
                  allocate(sing(nsing,nsing),trip(ntrip,ntrip))
                  if (.not.half) then
                    allocate(sing2(nsing,nsing),trip2(ntrip,ntrip))
                  else
                    sing2 => xdummy
                    trip2 => xdummy
                  end if
                  if (get_u) then
                    allocate(sing3(nsing,nsing),trip3(ntrip,ntrip))
                  else
                    sing3 => xdummy
                    trip3 => xdummy
                  end if
                  call spinsym_traf(1,ndim,scratch,flipmap_c,nsing,
       &                            sing,trip,.false.)
  
                  ! calculate T^(-0.5) for both blocks
                  if (iprint.ge.100) write(lulog,*) '+ case, nsing=',nsing
                  call invsqrt_mat(nsing,sing,sing2,half,sing3,get_u,
       &                           svs,icnt_sv,icnt_sv0,
       &                           xmax,xmin,bins(1,iexc_cls))
                  if (iprint.ge.100) write(lulog,*) '- case, ntrip=',ntrip
                  call invsqrt_mat(ntrip,trip,trip2,half,trip3,get_u,
       &                           svs(nsing+min(1,ntrip)),!avoid segfault
       &                           icnt_sv,icnt_sv0,
       &                           xmax,xmin,bins(1,iexc_cls))
  
                  ! partial undo of pre-diagonalization: Upre*T^(-0.5)
                  call spinsym_traf(2,ndim,scratch,flipmap_c,nsing,
       &                            sing,trip,.true.)
                  if (.not.half) then
                    ! full undo of pre-diagonalization for projector
                    call spinsym_traf(2,ndim,scratch2,flipmap_c,nsing,
       &                              sing2,trip2,.false.)
                    deallocate(sing2,trip2)
                  end if
                  if (get_u) then
                    ! partial undo of pre-diagonalization: Upre*U
                    call spinsym_traf(2,ndim,scratch3,flipmap_c,nsing,
       &                              sing3,trip3,.true.)
                    deallocate(sing3,trip3)
                  end if
                  deallocate(sing,trip)
               else ! msc.eq.0 .and. .not.lmodspc
  
                 if (lmodspc) then
                    call mat_svd_traf(ndim,scratch,scratch2,
       &                 icnt_sv,icnt_sv0,xmax,xmin,
       &                 bins(1,iexc_cls))
                 else
                    ! calculate S^(-0.5)
                    call invsqrt_mat(ndim,scratch,scratch2,
       &                             half,scratch3,get_u,svs,
       &                             icnt_sv,icnt_sv0,xmax,xmin,
       &                             bins(1,iexc_cls))
                 end if
               end if ! msc.eq.0 .and. .not.lmodspc
  
                ! Tikhonov regularization?
                if (reg_tik.and..not.lmodspc)
       &           call regular_tikhonov(ndim,ndim,scratch,svs,omega2)
  
                ! write to output buffer (contains transformation matrix X)
                call copy_buffer_2_h(scratch, buffer_out,
       &             ndim,ndim,ioff,transp)
  
                deallocate(scratch,svs)
  
                if (.not.half.or.lmodspc) then
  ! write projector to input buffer
  c dbg
                  if (ntest.ge.150) write(lulog,*) 
       &                'wrote prj to buffer_in'
  c dbg
                   ! write to ouptut buffer (contains projector)
                   call copy_buffer_2_h(scratch2, buffer_in, 
       &                ndim,ndim,ioff,transp)
                   deallocate(scratch2)
                end if
  
                if (get_u) then
  ! write unitary matrix to buffer
                   call copy_buffer_2_h(scratch3, buffer_u,
       &                ndim,ndim,ioff,transp)
  
                   deallocate(scratch3)
                end if
             enddo igama_loop
          enddo msa_loop
  
            deallocate(hpvx_csub,hpvx_asub,
       &             occ_csub, occ_asub,
       &             graph_csub, graph_asub,
       &             msdis_c,  msdis_a,
       &             idxmsdis_c,  idxmsdis_a,
       &             msdis_c2,  msdis_a2,
       &             idxmsdis_c2,  idxmsdis_a2,
       &             gamdis_c, gamdis_a,
       &             len_str,
       &             istr_csub,istr_asub,
       &             istr_csub_flip,istr_asub_flip,
       &             ldim_opin_c,ldim_opin_a)
            ifree = mem_flushmark('invsqrt_blk')
            if (sgrm.and.icnt_cur.lt.icnt_sv-icnt_sv0)
       &         blk_redundant(iocc_cls) = .false.
  
            cycle iocc_loop ! skip rest of the loop
         end if !onedis, i.e. the simple cases
  
  
          ! Here comes the complicated part for densities with 3 vertices
          ! and multiple distributions per Ms(A)/GAMMA(A) block
          ! Here also elements from different Ms(A)/GAMMA(A) blocks couple
          ! since we need to sort the whole occupation class into
          ! vectors of certain A1/C1 tuples
          ! i.e. the matrix must be resorted from A | C to A1/C1 | A2/C2
  
          ! determine number of c/a-operators of A1/C1 tuple
          na1mx = sum(hpvx_occ(1:ngastp,2,iblkoff+1))
          nc1mx = sum(hpvx_occ(1:ngastp,1,iblkoff+2))
          na2 = sum(hpvx_occ(1:ngastp,2,iblkoff+2))
          nc2 = sum(hpvx_occ(1:ngastp,1,iblkoff+3))
          msmax_sub = na1mx + nc1mx
          ! must be symmetric and so on and so forth
          if (njoined.ne.3.or.msmax_sub.ne.na2+nc2.or.na1mx.ne.nc2.or.
       &      na2.ne.nc1mx)
       &     call quit(1,'invsqrt','too difficult operator!')
          msmax = op_inp%ica_occ(1,iocc_cls)
          mscmax = op_inp%ica_occ(2,iocc_cls)
          if (msmax.ne.mscmax.or.msmax.ne.na1mx+na2
       &      .or.mscmax.ne.nc1mx+nc2)
       &            call quit(1,'invsqrt','need part.cons.op')
  
          ! occupation and so on of A2/C2 vector
          nel = na2 + nc2
          allocate(idorb(nel),idspn(nel),idspc(nel),lexlscr(nel,3),
       &           iocc(ngastp,2,2),idx_g(ngastp,2,2),msdst(ngastp,2,2),
       &           igamdst(ngastp,2,2))
          iocc = 0
          idx_g = 0
          iocc(1:ngastp,2,1) = hpvx_occ(1:ngastp,2,iblkoff+2)
          iocc(1:ngastp,1,2) = hpvx_occ(1:ngastp,1,iblkoff+3)
          idx_g(1:ngastp,2,1) = idx_graph(1:ngastp,2,iblkoff+2)
          idx_g(1:ngastp,1,2) = idx_graph(1:ngastp,1,iblkoff+3)
  
  
          ! loop over Ms/Gamma combinations of A1/C1 tuple
          ! which are the decoupled blocks of the A1/C1 | A2/C2 matrix
          ! Ms/Gamma of C2/A2 is then already defined
  ! For now we consider densities: Ms(total)=0, Gamma(total) = 1
  
          do ms1 = msmax_sub, -msmax_sub, -2
             ms2 = -ms1           ! particle conserving operator
  
           do igam = 1, ngam ! irrep of the second dimension must be the same
              call calculate_rank_dimensions_h(
       &           mel_inp,op_inp,graphs,iocc_cls, njoined,
       &           ms1, igam,
       &           sgrm,project,
       &           na1mx,nc1mx,
       &           ndim,nrank,rankoff,rankdim)
  
            if (ndim.eq.0) cycle
  
            if (project.eq.0) then
              nrank = 1
              rankdim(1) = ndim
              rankoff(1) = 0
            else if (ndim.ne.rankoff(1)+rankdim(1)) then
              write(lulog,'(1x,a,3i12)')
       &        'ndim,rankoff(1),rankdim(1): ',ndim,rankoff(1),rankdim(1)
              call quit(1,i_am,'dimensions don''t add up!')
            end if
  
            if (iprint.ge.15)
       &       write(lulog,'(a,3i8)') 'ms1, igam, ndim:',ms1,igam,ndim
            if (iprint.ge.15)
       &       write(lulog,'(a,5i8)') 'dim. per rank:',rankdim(1:nrank)
  
            allocate(scratch(ndim,ndim),flmap(ndim,3),svs(ndim))
  
            scratch = 0d0
            svs = 0d0
            if (.not.half.or.lmodspc) then
              allocate(scratch2(ndim,ndim))
              ! initialize: important for off-dia. blks for GNO/seq.orth.
              scratch2(1:ndim,1:ndim) = 0d0
            end if
            if (get_u) allocate(scratch3(ndim,ndim))
            if (get_u) scratch3(1:ndim,1:ndim) = 0d0
  c dbg
  c          allocate(matrix(ndim,ndim))
  c dbgend
  
            flmap(1:ndim,1:3) = 1
            ! loops over coupling blocks. Must be in correct order!
            jocc_cls = iocc_cls
  
            if (ntest.ge.100)
       &         write(lulog,'(" looking for coupling blocks")')
  
            jblkoff = (jocc_cls-1)*njoined
            na1 = na1mx
            nc1 = nc1mx
            off_line2 = 0
            do while(min(na1,nc1).ge.0.and.na1+nc1.ge.abs(ms1))
             na2 = nc1mx
             nc2 = na1mx
             off_col2 = 0
             off_linmax = 0
             do while(min(na2,nc2).ge.0.and.na2+nc2.ge.abs(ms1))
              ! no off-diagonal blocks in SepO (separate orthog.)
              if (sgrm.and.na1.ne.nc2) then
                na2 = na2 - 1
                nc2 = nc2 - 1
                off_col2 = off_colmax
                cycle
              end if
              ! exit if there is no block like this
              if (na1.ne.sum(hpvx_occ(1:ngastp,2,jblkoff+1)).or.
       &          nc1.ne.sum(hpvx_occ(1:ngastp,1,jblkoff+2)).or.
       &          na2.ne.sum(hpvx_occ(1:ngastp,2,jblkoff+2)).or.
       &          nc2.ne.sum(hpvx_occ(1:ngastp,1,jblkoff+3))) exit
              blk_used(jocc_cls) = .true.
              ! exception for pure inactive block:
              if (na1+nc1+na2+nc2.eq.0) then
                if (igam.ne.1) exit
                ioff = mel_inp%off_op_gmox(jocc_cls)%
       &                 d_gam_ms(1,1,1)
                scratch(off_line2+1,off_col2+1) = buffer_in(ioff+1)
                flmap(off_col2+1,1:2) = ioff+1
                flmap(off_col2+1,3) = 1
                exit
              end if
  c dbg
  c            print *,'na1,nc1,na2,nc2,jocc_cls:',na1,nc1,na2,nc2,
  c     &          jocc_cls
  c            print *,'off_line2, off_col2: ',off_line2,off_col2
  c dbgend
              msmax = op_inp%ica_occ(1,jocc_cls)
              call get_num_subblk(ncblk2,nablk2,
       &           hpvx_occ(1,1,jblkoff+1),njoined)
              allocate(hpvx_csub2(ncblk2),hpvx_asub2(nablk2),
       &               occ_csub2(ncblk2), occ_asub2(nablk2),
       &               graph_csub2(ncblk2), graph_asub2(nablk2),
       &               iocc2(ngastp,2,2),idx_g2(ngastp,2,2))
              iocc2 = 0
              idx_g2 = 0
              iocc2(1:ngastp,2,1) = hpvx_occ(1:ngastp,2,jblkoff+2)
              iocc2(1:ngastp,1,2) = hpvx_occ(1:ngastp,1,jblkoff+3)
              idx_g2(1:ngastp,2,1) = idx_graph(1:ngastp,2,jblkoff+2)
              idx_g2(1:ngastp,1,2) = idx_graph(1:ngastp,1,jblkoff+3)
              ! set HPVX and OCC info
              call condense_occ(occ_csub2, occ_asub2,
       &                      hpvx_csub2,hpvx_asub2,
       &                      hpvx_occ(1,1,jblkoff+1),njoined,hpvxblkseq)
              ! do the same for the graph info
              call condense_occ(graph_csub2, graph_asub2,
       &                      hpvx_csub2,hpvx_asub2,
       &                      idx_graph(1,1,jblkoff+1),njoined,hpvxblkseq)
  
              ! read in current Ms1/GAMMA block
              ! loop over Ms(A1)/GAMMA(A1) --> Ms(C1)/GAMMA(C1) already defined
              off_line = off_line2
              off_colmax = 0
  
              call init_ms_generator_h(ms_gen1,na1,ms1)
              do while(next_ms_h(ms_gen1,msa1,msc1))
                 if (.not. is_possible_h(nc1,msc1,1) ) cycle
  
               call init_gam_generator_h(gam_gen1,ngam,igam)
               do while(next_gam_h(gam_gen1, gama1,gamc1))
                if (.not. is_possible_h(na1,0,gama1) ) cycle
                if (.not. is_possible_h(nc1,0,gamc1) ) cycle
  
                off_col = off_col2
  
  ! loop over Ms(C2)/GAMMA(C2) --> Ms(A2)/GAMMA(A2) already defined
                call init_ms_generator_h(ms_gen2,nc2,-ms2)
                do while(next_ms_h(ms_gen2,msc2,msa2) )
                   if (.not. is_possible_h(na2, msa2,1)) cycle
  
                 msa = msa1 + msa2
                 idxmsa = (msmax-msa)/2 + 1
                 idxmsa2 = msa2idxms4op(-msa,ms1+ms2,msmax,msmax)
  
                 call init_gam_generator_h(gam_gen2,ngam,igam)
                 do while(next_gam_h(gam_gen2, gamc2,gama2))
                  if (.not. is_possible_h(na2,0,gama2) ) cycle
                  if (.not. is_possible_h(nc2,0,gamc2) ) cycle
  
                  igama = multd2h(gama1,gama2)
  
                  call set_dist_msgam_h(
       &               gamdis_c,gamdis_a,
       &               msdis_c,msdis_a,
       &               msc1,msa1,
       &               msc2,msa2,
       &               gamc1,gama1,
       &               gamc2,gama2,
       &               nc1,na1,
       &               nc2,na2,
       &               ncblk2,nablk2)
  
  
                  ! determine the distribution in question
                  call ms2idxms(idxmsdis_c,msdis_c,occ_csub2,ncblk2)
                  call ms2idxms(idxmsdis_a,msdis_a,occ_asub2,nablk2)
  
                  call set_len_str(len_str,ncblk2,nablk2,
       &                       graphs,
       &                       graph_csub2,idxmsdis_c,gamdis_c,hpvx_csub2,
       &                       graph_asub2,idxmsdis_a,gamdis_a,hpvx_asub2,
       &                       hpvxseq,.false.)
                  lenca = ielprd(len_str,ncblk2+nablk2) ! yeah, a product
                  call expand_len_str_h(len2,len_str,nablk2,ncblk2,nc1,nc2
       &               ,na1,na2)
                  if (lenca.eq.0) cycle
  
                  ndis = mel_get_ndis_h(mel_inp,jocc_cls,igama,idxmsa)
                  idxdis = 1
                  if (ndis.gt.1)
       &             idxdis =
       &                 idx_msgmdst2(.true.,
       &                  jocc_cls,idxmsa,igama,
       &                  occ_csub2,idxmsdis_c,gamdis_c,ncblk2,
       &                  occ_asub2,idxmsdis_a,gamdis_a,nablk2,
       &                  .false.,-1,-1,mel_inp,ngam)
  c dbg
  c                print *,'igama,idxmsa: ',igama,idxmsa
  c                print *,'lenca,ndis,idxdis,len:',lenca,ndis,idxdis,
  c     &           mel_inp%len_op_gmox(jocc_cls)%
  c     &               d_gam_ms(idxdis,igama,idxmsa)
  c dbgend
                  if (lenca.ne.mel_inp%len_op_gmox(jocc_cls)%
       &               d_gam_ms(idxdis,igama,idxmsa))
       &             call quit(1,i_am,'inconsistency!')
  
                  call set_op_ldim_c(ldim_opin_c,ldim_opin_a,
       &               hpvx_csub2,hpvx_asub2,
       &               len_str,ncblk2,nablk2,.false.)
  
                  ioff = mel_inp%off_op_gmox(jocc_cls)%
       &                 d_gam_ms(idxdis,igama,idxmsa)
  
                  ! now for spin-flipped counterpart:
                  msdis_c2(1:ncblk2) = -msdis_c(1:ncblk2)
                  msdis_a2(1:nablk2) = -msdis_a(1:nablk2)
                  call ms2idxms(idxmsdis_c2,msdis_c2,occ_csub2,ncblk2)
                  call ms2idxms(idxmsdis_a2,msdis_a2,occ_asub2,nablk2)
  
                  ndis2 =mel_inp%off_op_gmox(jocc_cls)%ndis(igama,idxmsa2)
                  idxdis2 = 1
                  if (ndis2.gt.1)
       &               idxdis2 =
       &                   idx_msgmdst2(.true.,
       &                    jocc_cls,idxmsa2,igama,
       &                    occ_csub2,idxmsdis_c2,gamdis_c,ncblk2,
       &                    occ_asub2,idxmsdis_a2,gamdis_a,nablk2,
       &                    .false.,-1,-1,mel_inp,ngam)
  
                  ioff2 = mel_inp%off_op_gmox(jocc_cls)%
       &                 d_gam_ms(idxdis2,igama,idxmsa2)
  
                  ! set flip maps (if not existing yet)
                  call strmap_man_flip(
       &               maxbuf_tmp,graph_csub2,ncblk2,
       &               str_info,strmap_info,orb_info)
                  if (maxbuf_tmp.gt.maxbuf) call quit(1,i_am,
       &             'add more to maxbuf')
                  call strmap_man_flip(
       &               maxbuf_tmp,graph_asub2,nablk2,
       &               str_info,strmap_info,orb_info)
                  if (maxbuf_tmp.gt.maxbuf) call quit(1,i_am,
       &             'add more to maxbuf')
  
                  call get_flipmap_blk(flipmap_c,
       &              ncblk2,occ_csub2,len_str,
       &              graph_csub2,idxmsdis_c,gamdis_c,
       &              strmap_info,ngam,ngraph)
                  call get_flipmap_blk(flipmap_a,
       &              nablk2,occ_asub2,len_str(ncblk2+1),
       &              graph_asub2,idxmsdis_a,gamdis_a,
       &              strmap_info,ngam,ngraph)
  
                  ! assemble distribution of A2/C2 tuple
                  ! inelegant: we are currently using next_tupel_ca for setting up
                  ! final flipmap. But there should be a way to do this using
                  ! only the flipmaps from above.
                  msdst = 0
                  igamdst = 1
                  do idx = 1, ngastp
                    if (iocc2(idx,2,1).ne.0) then
                      msdst(idx,2,1) = msdis_a(nablk2)
                      igamdst(idx,2,1) = gamdis_a(nablk2)
                    end if
                    if (iocc2(idx,1,2).ne.0) then
                      msdst(idx,1,2) = msdis_c(ncblk2)
                      igamdst(idx,1,2) = gamdis_c(ncblk2)
                    end if
                  end do
  
  c dbg
  c                write(lulog,'(a,4i4)') 'ms  : ',msa1,msc1,msa2,msc2
  c                write(lulog,'(a,4i4)') 'gam : ',gama1,gamc1,gama2,gamc2
  c                write(lulog,'(a,2i4)') 'msa, igama: ',msa,igama
  c                write(lulog,'(a,2i4)') 'dist, len: ',idxdis, lenca
  c                write(lulog,'(a,1i4)') 'dist2: ',idxdis2
  c                write(lulog,'(a,2i4)') 'off_line/col: ',off_line,off_col
  c                print *,'len1: ',len_str(1)*len_str(3)
  c                print *,'flipmap_c: len=',len_str(1:ncblk2)
  c                print '(10i6)',flipmap_c(1:sum(len_str(1:ncblk2)))
  c                print *,'flipmap_a:'
  c                print '(10i6)',
  c     &               flipmap_a(1:sum(len_str(ncblk2+1:ncblk2+nablk2)))
  c dbgend
  
                  ! copy all required elements of this distribution
                  ! to their block in the A1/C1 | A2/C2 matrix
                  first = .true.
                  iline = off_line
                  do idxc1 = 1, len2(1)
                    if (nc1.ne.0) then
                     istr_csub(1) = idxc1-1
                     istr_csub_flip(1) = abs(flipmap_c(idxc1))-1
                    end if
                    do idxa1 = 1, len2(3)
                      if (na1.ne.0) then
                       istr_asub(1) = idxa1-1
                       istr_asub_flip(1) = abs(flipmap_a(idxa1))-1
                      end if
                      iline = iline + 1
                      icol = off_col
                      do idxa2 = 1, len2(4)
                        if (na1.eq.0.and.na2.ne.0) then
                         istr_asub(1) = idxa2-1
                         istr_asub_flip(1) =
       &                          abs(flipmap_a(idxa2))-1
                        else if (na2.ne.0) then
                         istr_asub(2) = idxa2-1
                         istr_asub_flip(2) =
       &                          abs(flipmap_a(len2(3)+idxa2))-1
                        end if
                        do idxc2 = 1, len2(2)
                          if (nc1.eq.0.and.nc2.ne.0) then
                           istr_csub(1) = idxc2-1
                           istr_csub_flip(1) =
       &                          abs(flipmap_c(idxc2))-1
                          else if (nc2.ne.0) then
                           istr_csub(2) = idxc2-1
                           istr_csub_flip(2) =
       &                          abs(flipmap_c(len2(1)+idxc2))-1
                          end if
                          icol = icol + 1
                          idx = ioff + idx_str_blk3(istr_csub,istr_asub,
       &                         ldim_opin_c,ldim_opin_a,ncblk2,nablk2)
                          scratch(iline,icol) = buffer_in(idx)
  c dbg
  c                        matrix(iline,icol) = idx
  c dbgend
                          if (ms1.eq.0.and.iline.eq.icol) then
                            idx2 = ioff2 + idx_str_blk3(istr_csub_flip,
       &                           istr_asub_flip,
       &                           ldim_opin_c,ldim_opin_a,ncblk2,nablk2)
                            flmap(icol,1) = idx
                            flmap(icol,2) = idx2
                            logdum = next_tupel_ca(idorb,idspn,idspc,
       &                     na2+nc2,2,iocc2,idx_g2,
       &                     msdst,igamdst,first,
       &                     str_info%igas_restr,
       &                     orb_info%mostnd,orb_info%igamorb,
       &                     orb_info%nsym,orb_info%ngas,
       &                     orb_info%ngas_hpv,orb_info%idx_gas,
       &                     hpvxseq,lexlscr)
                            first = .false.
                            if (.not.logdum) call quit(1,'invsqrt',
       &                         'no next tuple found!')
                            if (mod(na2+nc2
       &                            -idxcount(2,idspn,na2+nc2,1),4).ne.0)
       &                          flmap(icol,3) = -1
  c dbg
  c                          write(lulog,'(i8,x,4i4,x,4i4)')idx,idorb,idspn
  c dbgend
                          end if
                        end do
                      end do
                    end do
                  end do
  
                  off_col = off_col + len2(2)*len2(4)
                  off_colmax = max(off_colmax,off_col)
                 end do
                end do
                off_line = off_line + len2(1)*len2(3)
                off_linmax = max(off_linmax,off_line)
               end do
              end do
  
              deallocate(hpvx_csub2,hpvx_asub2,occ_csub2,
       &              occ_asub2,graph_csub2,graph_asub2,
       &           iocc2,idx_g2)
              call advance_block_n2_h(
       &           nc2,na2,
       &           ms1,
       &           nc1mx,na1mx,
       &           jocc_cls)
              jblkoff = (jocc_cls-1)*njoined
              off_col2 = off_colmax
             end do
             na1 = na1 - 1
             nc1 = nc1 - 1
             off_line2 = off_linmax
            end do
  c dbg
            if (ntest.ge.100) then
                write(lulog,*) 'initial overlap matrix: jocc_cls = ',
       &                        jocc_cls
                call wrtmat3(scratch,ndim,ndim,ndim,ndim)
            end if
  c dbgend
  c dbg
  c          print *,'index matrix:'
  c          write(*,'(5x,18i4)') (icol, icol=1,ndim)
  c          do iline = 1, ndim
  c            write(*,'(i4,x,18i4)') iline, matrix(iline,1:ndim)
  c          end do
  c dbgend
            if (ntest.ge.150) write(lulog,*) 'entering irank loop'
  
            ! loop over blocks that should be orthogonalized separately
            do irank = 1, nrank
             icnt_cur = icnt_sv - icnt_sv0
             rdim = rankdim(irank)
             idxst = rankoff(irank) + 1
             idxnd = rankoff(irank) + rdim
             if (ntest.ge.100) write(lulog,*) '---- irank = ',irank
  
             ! build projector and apply to current block
             if (irank.ge.2) then
              if (iprint.ge.15)
       &        write(lulog,'(x,a,i8)') 'next rank:',irank
              if (ntest.ge.100) write(lulog,*) 'building projector'
  
              select case(project)
              case (1,2)
               if (project.eq.1) then
                ! write tensors "p" on off-diagonal blocks of scratch
                do jrank = 1, irank-1
                  rdim2 = rankdim(jrank)
                  idxst2 = rankoff(jrank) + 1
                  idxnd2 = rankoff(jrank) + rdim2
                 call create_submatrix_h(
       &               idxst,idxnd,
       &               idxst2,idxnd2,
       &               scratch_tmp1)
                 scratch_tmp1=0d0
                 call extract_submatrix_h(scratch,
       &              ndim, idxst2, idxnd2,
       &              ndim, idxst2, idxnd2,
       &              scratch_tmp2)
                  call get_vec2remove(scratch_tmp1,
       &                             scratch_tmp2,
       &                             rdim,rdim2,ms1,igam,irank-jrank,
       &                             na1mx-nrank+irank,nc1mx-nrank+irank,
       &                             gno,pass_spc,mel_spc,na1mx,nc1mx,
       &                             orb_info,str_info,strmap_info)
                 call delete_submatrix_h(
       &              scratch_tmp2
       &              )
                 call insert_submatrix_h(scratch,
       &              ndim, idxst, idxnd,
       &              ndim, idxst2, idxnd2,
       &              scratch_tmp1)
  
                end do
               else ! if project.eq.2
                ! multiply off-diagonal blocks with lower trafo matrices
                do jrank = 1, irank-1
                 rdim2 = rankdim(jrank)
                 idxst2 = rankoff(jrank) + 1
                 idxnd2 = rankoff(jrank) + rdim2
  ! Evil Hacking ahead:
  ! this actually is scratch(idxst:idxnd,idxst2:idxnd2)=scratch(idxst2:idxnd2,idxst:idxnd)*scratch(idxst2:idxnd2,idxst2:idxnd2)
  ! But defining it with ranges creates temorary matrices on the stack (at least with ifort) leading stackoverflow
  ! we could create those matrices on the heap and our elements. (we do this for other routines)
  ! to avoid copying : we give the first element of the submatrix in question and tell dgemm, that
  !                  scratch(idxst2,idxst+1) is ndim elements behind scratch(idxst2,idxst)
  !  this is how dgemm is intended to be used.
  ! I miss C++
  
                 call dgemm('t','n',rdim,rdim2,rdim2,
       &                   1d0,scratch(idxst2,idxst),ndim,
       &                   scratch(idxst2,idxst2),ndim,
       &                   0d0,scratch(idxst,idxst2),ndim)
  
  c dbg
                 ! transpose (only needed for double-check below)
                 do idx = 0, rdim2-1
                   scratch(idxst2+idx,idxst:idxnd)
       &           = scratch(idxst:idxnd,idxst2+idx)
                 end do
  c dbgend
                end do
                ! singular value decompose off-diagonal block
                ! and overwrite by nonredundant part of W^+ (left-hand)
                call svd_get_left(rdim,ndim-idxnd,
       &                          scratch(idxst:idxnd,idxnd+1:ndim)) ! TODO avoid stack overflow
               end if
               allocate(proj(rdim,rdim))
               if (project.eq.1) then
                allocate(scratch4(rdim,rdim))
  ! form outer product of these vectors
                !print *, "to be tested or not to be tested"
  
  ! Evil Hacking ahead:
  ! see above
                call dgemm('n','t',rdim,rdim,ndim-idxnd,
       &                   1d0,scratch(idxst,idxnd+1),ndim,
       &                   scratch(idxst,idxnd+1),ndim,
       &                   0d0,scratch4,rdim)
                ! assemble projector: 1 - sum(p*p^T)*S
                call dgemm('n','n',rdim,rdim,rdim,
       &                   -1d0,scratch4,rdim,
       &                   scratch(idxst,idxst),ndim,
       &                   0d0,proj,rdim)
               else ! if project.eq.2
  ! form outer product -W^+ * W
  ! Evil Hacking ahead:
  ! see above
                call dgemm('n','t',rdim,rdim,ndim-idxnd,
       &                   -1d0,scratch(idxst,idxnd+1),ndim,
       &                   scratch(idxst,idxnd+1),ndim,
       &                   0d0,proj,rdim)
               end if
               do idx = 1, rdim
                 proj(idx,idx) = proj(idx,idx) + 1d0
               end do
               if (iprint.ge.100) then
                 write(lulog,*) 'Projector for removing lower-rank exc.:'
                 call wrtmat3(proj,rdim,rdim,rdim,rdim)
               end if
               ! Apply projector to current metric block
               ! (a) Q^+*S
               !     This step is not strictly necessary, but numerically
               !     it helps to keep the matrix S*Q symmetric
               if (project.eq.2) allocate(scratch4(rdim,rdim))
               call dgemm('t','n',rdim,rdim,rdim,
       &                  1d0,proj,rdim,
       &                  scratch(idxst,idxst),ndim,
       &                  0d0,scratch4,rdim)
               ! (b) (Q^+*S)*Q
               call dgemm('n','n',rdim,rdim,rdim,
       &            1d0,scratch4,rdim,proj,rdim,
       &            0d0,scratch(idxst,idxst),ndim)
               deallocate(scratch4)
  c dbg
               if (project.eq.2) then
  ! Evil Hacking ahead:
  ! see above
                 ! check that off-diagonal blocks are projected out
                 call dgemm('t','t',rdim,ndim-idxnd,rdim,
       &                    1d0,proj,rdim,
       &                    scratch(idxnd+1,idxst),ndim,
       &                    0d0,scratch(idxst,idxnd+1),ndim)
                 if (iprint.ge.100) then
                   write(lulog,*) 'Projected off-diagonal block:'
                   call wrtmat3(scratch(idxst:idxnd,idxnd+1:ndim),
       &                        rdim,ndim-idxnd,rdim,ndim-idxnd)
                 end if
                 if (any(abs(scratch(idxst:idxnd,idxnd+1:ndim)).gt.1d-12))
       &           call warn('invsqrt','off-dia block not projected out?')
               end if
  c dbgend
               if (gno.eq.1.or.project.eq.2) then
                 ! zero off-diagonal blocks (assuming the projector works)
                 scratch(idxst:idxnd,idxnd+1:ndim) = 0d0
                 scratch(idxnd+1:ndim,idxst:idxnd) = 0d0
               end if
  
  c ==========================================
  c Here the usual Gram-Schmidt stuff happens:
  c ==========================================
              case (3,4)
               ! Gram-Schmidt step: substract lower-rank components from
               ! the diagonal block and the blocks above
  
               ! the way through all this is a bit involved; follow the marks
               ! (1) (2) (3) (4)
               ! (1) is at the end (look for "fourth mm")
  
               if (ntest.ge.100) write(lulog,*) 'GS projection follows'
  
               do jrank = irank, nrank ! target block: row j, column i
                rdim2 = rankdim(jrank)
                idxst2 = rankoff(jrank) + 1
                idxnd2 = rankoff(jrank) + rdim2
                do krank = 1, irank-1 ! input blocks: row k, columns j and i
                 rdim3 = rankdim(krank)
                 idxst3 = rankoff(krank) + 1
                 idxnd3 = rankoff(krank) + rdim3
  ! Evil Hacking ahead:
  ! see above
  c dbg
  c               write(lulog,'(" (",i4,",",i4,")x(",i4,",",i4,")'//
  c     &                    ' -> (",i4,",",i4,")")') 
  c     &                    idxst3,idxst2,idxst3,idxst,idxst2,idxst 
  c dbg
                 ! (2)
                 ! we reach this point for irank=2
                 ! scatch() contains up to this point: illustrated for doubles:
                 !/   S2     S21 \
                 !\ X1+.S12   X1 /
                 !
                 if (ntest.ge.150) write(lulog,'(" first mm for",2i4)') 
       &                       jrank, krank
                 if (ntest.ge.150) then 
                    write(lulog,'(" offsets and dimensions: ")') 
                    write(lulog,'(" matrix 1: ",4i6)') 
       &                     idxst3,idxst2,rdim3,rdim2 
                    write(lulog,'(" matrix 2: ",4i6)') 
       &                     idxst3,idxst,rdim3,rdim 
                    write(lulog,'(" matrix 3: ",4i6)') 
       &                     idxst2,idxst,rdim2,rdim 
                 end if
                 call dgemm('t','n',rdim2,rdim,rdim3,
       &                  -1d0,scratch(idxst3,idxst2),ndim,  ! (X1+.S12)+ = S21.X1
       &                  scratch(idxst3,idxst),ndim,        ! X1.S12
       &                  1d0,scratch(idxst2,idxst),ndim)    ! -> subtract from S2
                 ! state afterwards (S2' is projected S2):
                 !/   S2'    S21 \
                 !\ X1+.S12   X1 /
                end do
               end do
               ! compute the blocks below (intermediately stored on the right)
               ! for the formation of the trafo matrix
               do jrank = 1, irank-1 ! target block: row i, column j
                rdim2 = rankdim(jrank)
                idxst2 = rankoff(jrank) + 1
                idxnd2 = rankoff(jrank) + rdim2
                scratch(idxst:idxnd,idxst2:idxnd2) = 0d0 !initialize
                do krank = jrank, irank-1 ! input blocks: (k,i) and (j,k)
                 rdim3 = rankdim(krank)
                 idxst3 = rankoff(krank) + 1
                 idxnd3 = rankoff(krank) + rdim3
  ! Evil Hacking ahead:
  ! see above
  c dbg
  c               write(lulog,'(" (",i4,",",i4,")x(",i4,",",i4,")'//
  c     &                    ' -> (",i4,",",i4,")")') 
  c     &                    idxst3,idxst,idxst2,idxst3,idxst,idxst2 
  c dbg
                 if (ntest.ge.150) write(lulog,'(" second mm for",2i4)') 
       &                       jrank, krank
                 if (ntest.ge.150) then 
                    write(lulog,'(" offsets and dimensions: ")') 
                    write(lulog,'(" matrix 1: ",4i6)') 
       &                     idxst3,idxst,rdim3,rdim 
                    write(lulog,'(" matrix 2: ",4i6)') 
       &                     idxst2,idxst3,rdim3,rdim2 
                    write(lulog,'(" matrix 3: ",4i6)') 
       &                     idxst,idxst2,rdim,rdim2 
                 end if
                 ! (3)
                 call dgemm('t','t',rdim,rdim2,rdim3,
       &                  -1d0,scratch(idxst3,idxst),ndim, ! (X1+.S12)+ = S21.X1
       &                  scratch(idxst2,idxst3),ndim,     ! X1^+
       &                  1d0,scratch(idxst,idxst2),ndim)  ! --> -S21.X1.X1+
                 ! state afterwards (before now sending S2' to invsqrt_mat)
                 !/   S2'   -S21.X1.X1+ \
                 !\ X1+.S12     X1      /
                end do
               end do
  
              case default
                call quit(1,'invsqrt','unknown case for keyword project')
              end select
             end if
  c dbg
              if (iprint.ge.100) then
                write(lulog,*) 'overlap matrix after projections ',
       &                        jocc_cls
                call wrtmat3(scratch,ndim,ndim,ndim,ndim)
              end if
  c dbgend
            ! core step: singular value decomposition
            if (ms1.eq.0.and.rdim.gt.1.and..not.lmodspc) then
              ! here a splitting into "singlet" and "triplet" blocks is needed:
  c dbg
              if (ntest.ge.100) write(lulog,*) 
       &          'have to split matrix into "singlet" and "triplet"'
  c dbg
  
  c dbg
  c            write(lulog,*) 'flmap:'
  c            do icol = idxst, idxnd
  c              write(lulog,'(i4,2i6)') icol,flmap(icol,1:2)
  c            end do
  c dbgend
              do icol = idxst, idxnd
                idx = idxlist(flmap(icol,1),flmap(idxst:idxnd,2),rdim,1)
                if (idx.eq.-1) call quit(1,'invsqrt','idx not found!')
                flmap(icol,3) = flmap(icol,3)*idx
              end do
  c dbg
  c            write(lulog,*) 'flmap:'
  c            do icol = idxst, idxnd
  c              write(lulog,'(i4,3i6)') icol,flmap(icol,1:3)
  c            end do
  c dbgend
              nsing = rdim
              do icol = idxst, idxnd
                if (abs(flmap(icol,3)).eq.icol-idxst+1) nsing = nsing
       &                              + sign(1,flmap(icol,3))
              end do
              nsing = nsing/2
              ntrip = rdim - nsing
  
              ! do the pre-diagonalization
              allocate(sing(nsing,nsing),trip(ntrip,ntrip))
              if (.not.half) then
                allocate(sing2(nsing,nsing),trip2(ntrip,ntrip))
              else
                sing2 => xdummy
                trip2 => xdummy
              end if
              if (get_u) then
                allocate(sing3(nsing,nsing),trip3(ntrip,ntrip))
              else
                sing3 => xdummy
                trip3 => xdummy
              end if
              call extract_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
              call spinsym_traf(1,rdim,
       &                        scratch_tmp1,
       &                        flmap(idxst:idxnd,3),nsing,
       &                        sing,trip,.false.)
              call insert_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
              ! calculate T^(-0.5) for both blocks
              call invsqrt_mat(nsing,sing,sing2,half,sing3,get_u,
       &                       svs(idxst),icnt_sv,icnt_sv0,
       &                       xmax,xmin,bins(1,iexc_cls))
              call invsqrt_mat(ntrip,trip,trip2,half,trip3,get_u,
       &                       svs(idxst-1+nsing+min(1,ntrip)),!avoid segfault
       &                       icnt_sv,icnt_sv0,
       &                       xmax,xmin,bins(1,iexc_cls))
  
              ! partial undo of pre-diagonalization: Upre*T^(-0.5)
              call extract_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
              call spinsym_traf(2,rdim,
       &                        scratch_tmp1,
       &                        flmap(idxst:idxnd,3),nsing,
       &                        sing,trip,.true.)
              call insert_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
  
              if (.not.half) then
                ! full undo of pre-diagonalization for projector
                call extract_submatrix_h(scratch2,
       &             ndim, idxst, idxnd,
       &             ndim, idxst, idxnd,
       &             scratch_tmp1)
                call spinsym_traf(2,rdim,
       &                          scratch_tmp1,
       &                          flmap(idxst:idxnd,3),nsing,
       &             sing2,trip2,.false.)
              call insert_submatrix_h(scratch2,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
  
                deallocate(sing2,trip2)
              end if
              if (get_u) then
                ! partial undo of pre-diagonalization: Upre*U
              call extract_submatrix_h(scratch3,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
              call spinsym_traf(2,rdim,
       &                          scratch_tmp1,
       &                          flmap(idxst:idxnd,3),nsing,
       &                          sing3,trip3,.true.)
              call insert_submatrix_h(scratch3,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
              deallocate(sing3,trip3)
              end if
              deallocate(sing,trip)
  
            else if (lmodspc) then
  c dbg
              if (ntest.ge.100) write(lulog,*) 'lmodspc case'
  c dbg
  
              call extract_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
              call extract_submatrix_h(scratch2,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp2)
              call mat_svd_traf(rdim,scratch_tmp1,
       &                       scratch_tmp2,
       &                       icnt_sv,icnt_sv0,xmax,xmin,
       &                       bins(1,iexc_cls))
              call insert_submatrix_h(scratch2,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp2)
              call insert_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
  
            else if (.not.half.and.get_u) then
  ! calculate S^(-0.5)
  c dbg
              if (ntest.ge.100) write(lulog,*) '.not.half.and.get_u case'
  c dbg
              call extract_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
              call extract_submatrix_h(scratch2,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp2)
              call extract_submatrix_h(scratch3,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp3)
              call invsqrt_mat(rdim,scratch_tmp1,
       &                       scratch_tmp2,
       &                       half,
       &                       scratch_tmp3,get_u,
       &                       svs(idxst),icnt_sv,icnt_sv0,xmax,xmin,
       &                       bins(1,iexc_cls))
              call insert_submatrix_h(scratch3,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp3)
              call insert_submatrix_h(scratch2,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp2)
              call insert_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
            else if (.not.half) then
  c dbg
              if (ntest.ge.100) write(lulog,*) '.not.half case'
  c dbg
              call extract_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
              call extract_submatrix_h(scratch2,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp2)
              call invsqrt_mat(rdim,scratch_tmp1,! out: X
       &                       scratch_tmp2,     ! out: Projector X X-1
       &                       half,xdummy,get_u, !scratch3: dummy
       &                       svs(idxst),icnt_sv,icnt_sv0,xmax,xmin,
       &                       bins(1,iexc_cls))
              call insert_submatrix_h(scratch2, ! contains now Proj.
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp2)
              call insert_submatrix_h(scratch,  ! contains now X
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
            else if (get_u) then
  c dbg
              if (ntest.eq.100) write(lulog,*) 'get_u case'
  c dbg
              ! calculate S^(-0.5)
              call extract_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
              call extract_submatrix_h(scratch3,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp3)
              call invsqrt_mat(rdim,scratch_tmp1,
       &                       xdummy,half, !scratch2: dummy
       &                       scratch_tmp3,get_u,
       &                       svs(idxst),icnt_sv,icnt_sv0,xmax,xmin,
       &                       bins(1,iexc_cls))
              call insert_submatrix_h(scratch3,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp3)
              call insert_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
            else
  c dbg
              if (ntest.ge.100) write(lulog,*) 'else case'
  c dbg
              ! calculate S^(-0.5)
              call extract_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
              call invsqrt_mat(rdim,scratch_tmp1,
       &                       xdummy,half, !scratch2: dummy
       &                       xdummy,get_u, !scratch3: dummy
       &                       svs(idxst),icnt_sv,icnt_sv0,xmax,xmin,
       &                       bins(1,iexc_cls))
              call insert_submatrix_h(scratch,
       &           ndim, idxst, idxnd,
       &           ndim, idxst, idxnd,
       &           scratch_tmp1)
            end if
  
             select case(project)
             case (1,2)
              ! apply projector again: X = Q*U*s^(-0.5)
              if (irank.ge.2) then
                allocate(scratch4(rdim,rdim))
  ! Evil Hacking ahead:
  ! see above
                call dgemm('n','n',rdim,rdim,rdim,
       &                   1d0,proj,rdim,
       &                   scratch(idxst,idxst),ndim,
       &                   0d0,scratch4,rdim)
                scratch(idxst:idxnd,idxst:idxnd) = scratch4
                if (iprint.ge.100) then
                  write(lulog,*) 'Trafo matrix:'
                  call wrtmat3(scratch(idxst:idxnd,idxst:idxnd),
       &                       rdim,rdim,rdim,rdim)
                end if
                if (.not.half) then
  ! Evil Hacking ahead:
  ! see above
                  call dgemm('n','n',rdim,rdim,rdim,
       &                     1d0,proj,rdim,
       &                     scratch2(idxst,idxst),ndim,
       &                     0d0,scratch4,rdim)
                  scratch2(idxst:idxnd,idxst:idxnd) = scratch4
                  if (iprint.ge.100) then
                    write(lulog,*) 'Projector matrix:'
                    call wrtmat3(scratch2(idxst:idxnd,idxst:idxnd),
       &                         rdim,rdim,rdim,rdim)
                  end if
                end if
                deallocate(scratch4,proj)
              end if
  
             case (3,4)
              ! multiply off-diagonal blocks of trafo matrix with dia.blks
              ! and set upper off-diagonal blocks to zero
              do jrank = 1, irank-1 ! target: (j,i) input: (i,j) and (i,i)
               rdim2 = rankdim(jrank)
               idxst2 = rankoff(jrank) + 1
               idxnd2 = rankoff(jrank) + rdim2
  ! Evil Hacking ahead:
  ! see above
               ! for the doubles, this is our last stop
               ! state so far is:
               !/   X2    -S21.X1.X1+ \
               !\ X1+.S12     X1      /
               if (ntest.ge.150) write(lulog,'(" third mm for",2i4)') 
       &                       irank, jrank
               if (ntest.ge.150) then 
                  write(lulog,'(" offsets and dimensions: ")') 
                  write(lulog,'(" matrix 1: ",4i6)') 
       &                     idxst,idxst2,rdim,rdim2 
                  write(lulog,'(" matrix 2: ",4i6)') 
       &                     idxst,idxst,rdim,rdim 
                  write(lulog,'(" matrix 3: ",4i6)') 
       &                     idxst2,idxst,rdim2,rdim 
               end if
               call dgemm('t','n',rdim2,rdim,rdim,
       &            1d0,scratch(idxst,idxst2),ndim, ! -(S21.X1.X1+)+ = -X1.X1+.S12
       &            scratch(idxst,idxst),ndim,      ! X2
       &            0d0,scratch(idxst2,idxst),ndim) ! --> -X1.X1+.S12.X2
               if(project.eq.4)then
                 call dgemm('t','n',rdim2,rdim,rdim,
       &            1d0,scratch(idxst,idxst2),ndim,  ! -(S21.X1.X1+)+ = -X1.X1+.S12 
       &            scratch2(idxst,idxst),ndim,      ! P2
       &            0d0,scratch2(idxst2,idxst),ndim) ! --> -X1.X1+.S12.X2
               end if
               scratch(idxst:idxnd,idxst2:idxnd2) = 0d0
               ! so we have finally:
               !/   X2              0 \
               !\ -X1.X1+.S12.X2   X1 /
               ! and for the projector (with term in [] for project==4)
               !/    P2             0 \
               !\ [-X1.X1+.S12.P2] P1 /
              end do
              do jrank = irank+1, nrank ! target: (i,j) input: (i,i) and (j,i)
               rdim2 = rankdim(jrank)
               idxst2 = rankoff(jrank) + 1
               idxnd2 = rankoff(jrank) + rdim2
  
               ! (1) actually the first place that the program will reach
               ! for the doubles case, this will be reached for irank=1, jrank=2
               if (ntest.ge.150) write(lulog,'(" fourth mm for",i4)') 
       &                       irank, jrank
               if (ntest.ge.150) then 
                  write(lulog,'(" offsets and dimensions: ")') 
                  write(lulog,'(" matrix 1: ",4i6)') 
       &                     idxst,idxst,rdim,rdim 
                  write(lulog,'(" matrix 2: ",4i6)') 
       &                     idxst2,idxst,rdim,rdim2 
                  write(lulog,'(" matrix 3: ",4i6)') 
       &                     idxst,idxst2,rdim,rdim2 
               end if
  
               ! comments refer to doubles case, we get here in the first round (irank=1)
               call dgemm('t','t',rdim,rdim2,rdim,
       &                  1d0,scratch(idxst:idxnd,idxst:idxnd),rdim,  ! X1
       &                  scratch(idxst2:idxnd2,idxst:idxnd),rdim2,   ! S21
       &                  0d0,scratch(idxst:idxnd,idxst2:idxnd2),rdim)! -> X1+.S12
              ! state after this multiply is then
              !/   S2     S21 \
              !\ X1+.S12   X1 /
              end do ! jrank
  
             end select
   
             if (iprint.ge.100) then
                write(lulog,*) 
       &           'current state of overlap matrix, after irank = ',irank
                call wrtmat3(scratch,ndim,ndim,ndim,ndim)
             end if
  
             if (sgrm.and.icnt_cur.lt.icnt_sv-icnt_sv0)
  c     &        blk_redundant(iocc_cls+min(na1mx,nc1mx)+irank-nrank)
  c     &        blk_redundant(iocc_cls+irank-1)
       &        blk_redundant(iocc_cls+nrank-irank)
       &        = .false.
            end do ! irank
  
            ! Tikhonov regularization?
            if (reg_tik.and..not.lmodspc)
       &       call regular_tikhonov(ndim,ndim,scratch,svs,omega2)
  
  c dbg
            if (iprint.ge.100) then
                write(lulog,*) 'final transformation matrix:'
                call wrtmat3(scratch,ndim,ndim,ndim,ndim)
                write(lulog,*) 'final projection matrix:'
                call wrtmat3(scratch2,ndim,ndim,ndim,ndim)
            end if
  c dbgend
  
            ! write to output buffer
            ! loops over coupling blocks. Must be in correct order!
            jocc_cls = iocc_cls
            jblkoff = (jocc_cls-1)*njoined
            na1 = na1mx
            nc1 = nc1mx
            off_line2 = 0
            do while(min(na1,nc1).ge.0.and.na1+nc1.ge.abs(ms1))
             na2 = nc1mx
             nc2 = na1mx
             off_col2 = 0
             off_linmax = 0
             do while(min(na2,nc2).ge.0.and.na2+nc2.ge.abs(ms1))
              ! no off-diagonal blocks in SepO (separate orthog.)
              if (sgrm.and.na1.ne.nc2) then
                na2 = na2 - 1
                nc2 = nc2 - 1
                off_col2 = off_colmax
                cycle
              end if
              ! exit if there is no block like this
              if (na1.ne.sum(hpvx_occ(1:ngastp,2,jblkoff+1)).or.
       &          nc1.ne.sum(hpvx_occ(1:ngastp,1,jblkoff+2)).or.
       &          na2.ne.sum(hpvx_occ(1:ngastp,2,jblkoff+2)).or.
       &          nc2.ne.sum(hpvx_occ(1:ngastp,1,jblkoff+3))) exit
              ! exception for pure inactive block:
              if (na1+nc1+na2+nc2.eq.0) then
                if (igam.ne.1) exit
                ioff = mel_inp%off_op_gmox(jocc_cls)%
       &                 d_gam_ms(1,1,1)
                buffer_out(ioff+1) = scratch(off_line2+1,off_col2+1)
                if (.not.half.or.lmodspc) ! copy projector to input buffer
       &           buffer_in(ioff+1) = scratch2(off_line2+1,off_col2+1)
                if (get_u)then
                   buffer_u(ioff+1) = scratch3(off_line2+1,off_col2+1)
                end if
                exit
              end if
              msmax = op_inp%ica_occ(1,jocc_cls)
              call get_num_subblk(ncblk2,nablk2,
       &           hpvx_occ(1,1,jblkoff+1),njoined)
              allocate(hpvx_csub2(ncblk2),hpvx_asub2(nablk2),
       &               occ_csub2(ncblk2), occ_asub2(nablk2),
       &               graph_csub2(ncblk2), graph_asub2(nablk2))
              ! set HPVX and OCC info
              call condense_occ(occ_csub2, occ_asub2,
       &                      hpvx_csub2,hpvx_asub2,
       &                      hpvx_occ(1,1,jblkoff+1),njoined,hpvxblkseq)
              ! do the same for the graph info
              call condense_occ(graph_csub2, graph_asub2,
       &                      hpvx_csub2,hpvx_asub2,
       &                      idx_graph(1,1,jblkoff+1),njoined,hpvxblkseq)
  
              ! loop over Ms(A1)/GAMMA(A1) --> Ms(C1)/GAMMA(C1) already defined
              off_line = off_line2
              off_colmax = 0
              call init_ms_generator_h(ms_gen1,na1,ms1)
              do while(next_ms_h(ms_gen1,msa1,msc1))
                 if (.not. is_possible_h(nc1,msc1,1) ) cycle
  
               call init_gam_generator_h(gam_gen1,ngam,igam)
               do while(next_gam_h(gam_gen1, gama1,gamc1))
                if (.not. is_possible_h(na1,0,gama1) ) cycle
                if (.not. is_possible_h(nc1,0,gamc1) ) cycle
  
                ! loop over Ms(C2)/GAMMA(C2) --> Ms(A2)/GAMMA(A2) already defined
                off_col = off_col2
                call init_ms_generator_h(ms_gen2,nc2,-ms2)
                do while(next_ms_h(ms_gen2,msc2,msa2) )
                   if (.not. is_possible_h(na2, msa2,1)) cycle
                   msa = msa1 + msa2
                   idxmsa = (msmax-msa)/2 + 1
  
                 call init_gam_generator_h(gam_gen2,ngam,igam)
                 do while(next_gam_h(gam_gen2, gamc2,gama2))
                    if (.not. is_possible_h(na2,0,gama2) ) cycle
                    if (.not. is_possible_h(nc2,0,gamc2) ) cycle
  
                  igama = multd2h(gama1,gama2)
                  call set_dist_msgam_h(
       &               gamdis_c,gamdis_a,
       &               msdis_c,msdis_a,
       &               msc1,msa1,
       &               msc2,msa2,
       &               gamc1,gama1,
       &               gamc2,gama2,
       &               nc1,na1,
       &               nc2,na2,
       &               ncblk2,nablk2)
  
                  ! determine the distribution in question
                  call ms2idxms(idxmsdis_c,msdis_c,occ_csub2,ncblk2)
                  call ms2idxms(idxmsdis_a,msdis_a,occ_asub2,nablk2)
  
                  call set_len_str(len_str,ncblk2,nablk2,
       &                       graphs,
       &                       graph_csub2,idxmsdis_c,gamdis_c,hpvx_csub2,
       &                       graph_asub2,idxmsdis_a,gamdis_a,hpvx_asub2,
       &               hpvxseq,.false.)
  
                  call expand_len_str_h(len2,len_str,nablk2,ncblk2,nc1,nc2
       &               ,na1,na2)
                  if (ielprd(len_str,ncblk2+nablk2).eq.0) cycle
  
                  ndis = mel_inp%off_op_gmox(jocc_cls)%ndis(igama,idxmsa)
                  idxdis = 1
                  if (ndis.gt.1)
       &             idxdis =
       &                 idx_msgmdst2(.true.,
       &                  jocc_cls,idxmsa,igama,
       &                  occ_csub2,idxmsdis_c,gamdis_c,ncblk2,
       &                  occ_asub2,idxmsdis_a,gamdis_a,nablk2,
       &                  .false.,-1,-1,mel_inp,ngam)
  
                  call set_op_ldim_c(ldim_opin_c,ldim_opin_a,
       &               hpvx_csub2,hpvx_asub2,
       &               len_str,ncblk2,nablk2,.false.)
  
                  ioff = mel_inp%off_op_gmox(jocc_cls)%
       &                 d_gam_ms(idxdis,igama,idxmsa)
  
                  ! copy all required elements of this distribution
                  ! to their block in the A1/C1 | A2/C2 matrix1
                  iline = off_line
                  do idxc1 = 1, len2(1)
                    if (nc1.ne.0) istr_csub(1) = idxc1-1
                    do idxa1 = 1, len2(3)
                      if (na1.ne.0) istr_asub(1) = idxa1-1
                      iline = iline + 1
                      icol = off_col
                      do idxa2 = 1, len2(4)
                        if (na1.eq.0.and.na2.ne.0) then
                         istr_asub(1) = idxa2-1
                        else if (na2.ne.0) then
                         istr_asub(2) = idxa2-1
                        end if
                        do idxc2 = 1, len2(2)
                          if (nc1.eq.0.and.nc2.ne.0) then
                           istr_csub(1) = idxc2-1
                          else if (nc2.ne.0) then
                           istr_csub(2) = idxc2-1
                          end if
                          icol = icol + 1
                          idx = ioff + idx_str_blk3(istr_csub,istr_asub,
       &                         ldim_opin_c,ldim_opin_a,ncblk2,nablk2)
                          buffer_out(idx) = scratch(iline,icol)
                          if (.not.half.or.lmodspc) ! copy projector to input buffer
       &                     buffer_in(idx) = scratch2(iline,icol)
                          if (get_u)then
                             buffer_u(idx) = scratch3(iline,icol)
                          end if
                        end do
                      end do
                    end do
                  end do
  
                  off_col = off_col + len2(2)*len2(4)
                  off_colmax = max(off_colmax,off_col)
                 end do
                end do
                off_line = off_line + len2(1)*len2(3)
                off_linmax = max(off_linmax,off_line)
               end do
              end do
  
              deallocate(hpvx_csub2,hpvx_asub2,occ_csub2,
       &              occ_asub2,graph_csub2,graph_asub2)
              na2 = na2 - 1
              nc2 = nc2 - 1
              jocc_cls = jocc_cls + 1
              if (na2+nc2.lt.abs(ms1)) !jump to next line
       &                jocc_cls = jocc_cls + min(na1mx,nc1mx)
       &                -(na1mx+nc1mx-abs(ms1))/2
              jblkoff = (jocc_cls-1)*njoined
              off_col2 = off_colmax
             end do
             na1 = na1 - 1
             nc1 = nc1 - 1
             off_line2 = off_linmax
            end do
  
            deallocate(scratch,flmap,svs)
            if (.not.half.or.lmodspc) deallocate(scratch2)
            if (get_u) deallocate(scratch3)
  c dbg
  c          deallocate(matrix)
  c dbgend
  
           end do
          end do
  
          deallocate(hpvx_csub,hpvx_asub,
       &           occ_csub, occ_asub,
       &           graph_csub, graph_asub,
       &           msdis_c,  msdis_a,
       &           idxmsdis_c,  idxmsdis_a,
       &           msdis_c2,  msdis_a2,
       &           idxmsdis_c2,  idxmsdis_a2,
       &           gamdis_c, gamdis_a,
       &           len_str,
       &           istr_csub,istr_asub,
       &           istr_csub_flip,istr_asub_flip,
       &           ldim_opin_c,ldim_opin_a,
       &           idorb,idspn,idspc,lexlscr,
       &           iocc,idx_g,msdst,
       &           igamdst)
  
          ifree = mem_flushmark('invsqrt_blk')
  
        enddo iocc_loop
        deallocate(blk_used)
  
        ! Set n_sv, sv_max, sv_min from molpro_out.h for use in the molpro
        ! output.
        ! Assuming that half=false, when calling invsqrt for the
        ! orthogonalisation (and not preconditioner) step.
        if (.not. half) then
           n_sv = icnt_sv-icnt_sv0
           sv_min = xmin
           sv_max = xmax
        end if
  
        if (iprint.ge.5) then
          write(lulog,'(x,77("="))')
          write(lulog,'(x,a)')
       &       'Singular value histogram (by excitation classes)'
          if (lmodspc) then
            write(lulog,'(x,a,x,14i10)') 'class:', (idx,idx=1,iexc_cls)
          else if (associated(op_t)) then
            write(lulog,'(x,a,x,14i10)') 'n_h = ',
       &         op_t%ihpvca_occ(IHOLE,2,ex2occ_cls(1:iexc_cls))
            write(lulog,'(x,a,x,14i10)') 'n_p = ',
       &         op_t%ihpvca_occ(IPART,1,ex2occ_cls(1:iexc_cls))
          end if
          write(lulog,'(x,77("-"))')
          write(lulog,'(x,a,x,14i10)') ' 1E+00',bins(1,1:iexc_cls)
  c        if (lmodspc) then
  c          do idx = 2, 8
  c            write(lulog,'(x,a,i2.2,x,14i10)') ' 1E-',idx-1,
  c     &                                        bins(idx,1:iexc_cls)
  c          end do
  c          write(lulog,'(x,a,x,14i10)') '     0',bins(9,1:iexc_cls)
  c          do idx = 10, 16
  c            write(lulog,'(x,a,i2.2,x,14i10)') '-1E-',17-idx,
  c     &                                        bins(idx,1:iexc_cls)
  c          end do
  c          write(lulog,'(x,a,x,14i10)') '-1E+00',bins(17,1:iexc_cls)
  c          write(lulog,'(x,77("="))')
  c          write(lulog,'(x,i8,a,i8,a)') icnt_sv-icnt_sv0,
  c     &          ' out of ',icnt_sv,' eigenvalues were included'
  c          write(lulog,'(x,a,E19.10)')
  c     &          'The  largest excluded eigenvalue is ',xmax
  c          write(lulog,'(x,a,E19.10)')
  c     &          'The smallest included eigenvalue is ',xmin
  c        else
          do idx = 2, 16
            write(lulog,'(x,a,i2.2,x,14i10)') ' 1E-',idx-1,
       &                                      bins(idx,1:iexc_cls)
          end do
          write(lulog,'(x,a,x,14i10)') '     0',bins(17,1:iexc_cls)
          write(lulog,'(x,77("="))')
          write(lulog,'(x,i8,a,i8,a)') icnt_sv-icnt_sv0,
       &        ' out of ',icnt_sv,' singular values were included'
          write(lulog,'(x,a,E19.10)')
       &        'The  largest excluded singular value is ',xmax
          write(lulog,'(x,a,E19.10)')
       &        'The smallest included singular value is ',xmin
  c        end if
        end if
        deallocate(bins,ex2occ_cls)
  
        if (sgrm.and.any(blk_redundant(1:nocc_cls))) then
          ! Print out which blocks are redundant
          write(lulog,'(x,a)') 'There are redundant blocks:'
  c        write(lulog,'(x,a)') 'There are redundant blocks in T:'
          write(lulog,'(x,a,26i3)') 'Block #  :',(idx,idx=1,nocc_cls)
          write(lulog,'(x,a,26(2x,L1))') 'Redundant?',
       &                               (blk_redundant(idx),idx=1,nocc_cls)
          call get_argument_value('method.MR','svdonly',lval=svdonly)
          if (svdonly.and..not.lmodspc) then
            ! Excplicitly print restrictions for input file
            if (.not.associated(op_t)) call quit(0,'invsqrt',
       &         'inconsistency in program: fix reference to op_t')
            write(lulog,*)
            write(lulog,'(x,a)') 'Copy the following into the input file:'
            ih = -1
            ip = -1
            iexc = 0
            first = .false.
            do iocc_cls = 1, nocc_cls
              first2 = .false.
              if (ih.ne.op_t%ihpvca_occ(IHOLE,2,iocc_cls)) then
                ih = op_t%ihpvca_occ(IHOLE,2,iocc_cls)
                first = .true.
                first2 = .true.
              end if
              if (ip.ne.op_t%ihpvca_occ(IPART,1,iocc_cls)) then
                ip = op_t%ihpvca_occ(IPART,1,iocc_cls)
                first = .true.
                first2 = .true.
              end if
              if (first2) then
                ! how many blocks belong to this excitation class?
                jocc_cls = nocc_cls
                do while(op_t%ihpvca_occ(IHOLE,2,jocc_cls).ne.ih.or.
       &                 op_t%ihpvca_occ(IPART,1,jocc_cls).ne.ip)
                  jocc_cls = jocc_cls - 1
                end do
                nrank = jocc_cls - iocc_cls
                irank = 0
              else
                irank = irank + 1
              end if
              iexc = op_t%ica_occ(1,iocc_cls)
              if (blk_redundant(iocc_cls+nrank-2*irank).and.first) then
                first = .false.
                write(lulog,'(x,a,i1,a,i1,a,i1,a,i1,a,i1,a)')
       &          'MR excrestr=(',ih,',',ih,',',ip,',',ip,',1,',iexc-1,')'
              else if (.not.blk_redundant(iocc_cls+nrank-2*irank).and.
       &               .not.first) then
                write(lulog,'(x,a,i4)') 'Watch out: Non-redundant block:',
       &                                iocc_cls
                call warn('invsqrt',
       &                  'non-redundant block beyond redundant one?')
              end if
            end do
            write(lulog,*)
          end if
        end if
        deallocate(blk_redundant)
  
        if(.not.bufout)then
          curr_rec=ffinv%current_record
          len_rec =ffinv%length_of_record
          call put_vec(ffinv,buffer_out,(curr_rec-1)*len_rec+1,
       &       (curr_rec-1)*len_rec+nbuff)
        endif
        if(.not.bufin.and.(.not.half.or.lmodspc))then
          curr_rec=ffinp%current_record
          len_rec =ffinp%length_of_record
          ! return projector matrix on input list
          call put_vec(ffinp,buffer_in,(curr_rec-1)*len_rec+1,
       &       (curr_rec-1)*len_rec+nbuff)
        endif
        curr_rec=ffu%current_record
        len_rec =ffu%length_of_record
        if (.not.bufu.and.get_u)then
           call put_vec(ffu,buffer_u,
       &     (curr_rec-1)*len_rec+1,(curr_rec-1)*len_rec+nbuff)
        end if
        ifree = mem_flushmark('invsqrt')
  
        return
        contains
  
  !-----------------------------------------------------------------------!
  !!    determines if the ms symmetry is fixed.
  !!
  !!
  !-----------------------------------------------------------------------!
        function compare_iocc_cls(occ_descr, occ_def, njoined )
        integer, parameter ::
       &     maxdef = 64
        logical ::
       &     compare_iocc_cls
        character(len=*), intent(in) ::
       &     occ_descr
        integer,intent(in)::
       &     njoined
        integer::
       &     occ_def(ngastp,2,njoined)
        integer ::
       &     occ_def_cmp(ngastp,2,maxdef)
        logical, external ::
       &     iocc_equal_n
        integer::
       &     nocc_cmp,iocc_cmp
  
        call process_occ_descr(occ_def_cmp,
       &     nocc_cmp,occ_descr,njoined,maxdef)
        do iocc_cmp = 1,nocc_cmp
           if(iocc_equal_n(
       &        occ_def_cmp(1,1,(iocc_cmp-1)*njoined+1),.false.,
       &        occ_def(1,1,1),.false.,njoined))then
              compare_iocc_cls=.true.
           end if
        end do
        compare_iocc_cls=.false.
  
        end function
  !----------------------------------------------------------------------- !
  !!    determines if the ms symmetry is fixed.
  !!
  !!
  !-----------------------------------------------------------------------!
        function fix_ms_h(mel_inp, mel_inv, mel_u)
        implicit none
        logical ::
       &    fix_ms_h
        type(me_list),intent(in)::
       &     mel_inp, mel_inv, mel_u
  
        fix_ms_h = .false.
        if(mel_inp%fix_vertex_ms.or.mel_inv%fix_vertex_ms
       &   .or.mel_u%fix_vertex_ms)then
          fix_ms_h = mel_inp%fix_vertex_ms.and.mel_inv%fix_vertex_ms
       &           .and.mel_u%fix_vertex_ms
          if(.not.fix_ms_h) call quit(1,'invsqrt',
       &                              'fix ms or not?')
        endif
        end function
  !-----------------------------------------------------------------------!
  !!    tries to guess the cluster operator for the singular value histogramm
  !!
  !!    returns a pointer to the cluster operator or a null pointer if no guess could be generated
  !-----------------------------------------------------------------------!
        function get_cluster_op_h(op_info)
  !-----------------------------------------------------------------------!
        implicit none
        type(operator),pointer::
       &     get_cluster_op_h
        type(operator_info),intent(in)::
       &     op_info
        integer::
       &     idxop
  
        idxop = idx_oplist2('T',op_info)
  
        if (is_keyword_set('method.MRCI').gt.0)
       &      idxop = idx_oplist2('C',op_info)
  
        if (idxop.gt.0)then
           get_cluster_op_h => op_info%op_arr(idxop)%op
        else
           get_cluster_op_h => null()
        end if
  
        end function get_cluster_op_h
  !-----------------------------------------------------------------------!
  !!   copies a one dimensional buffer onto a two dimensional buffer
  !-----------------------------------------------------------------------!
        subroutine copy_buffer_1_h(buffer_in,buffer_out,
       &     dim1,dim2,ioff,transpose )
  !-----------------------------------------------------------------------!
        implicit none
        real(8),intent(in)::
       &     buffer_in(:)
        real(8),intent(inout)::
       &     buffer_out(:,:)
        integer,intent(in)::
       &     dim1,dim2,
       &     ioff
        logical,intent(in)::
       &     transpose
  
        integer::
       &     idx,jdx
  
        if(transpose)then
           do idx = 1,dim1
              do jdx = 1,dim2
                 buffer_out(jdx,idx) =  buffer_in(ioff+(idx-1)*dim2+jdx)
              end do
           end do
        else
           do idx = 1,dim1
              do jdx = 1,dim2
                 buffer_out(idx,jdx) =  buffer_in(ioff+(idx-1)*dim2+jdx)
              end do
           end do
        end if
        end subroutine
  !-----------------------------------------------------------------------!
  !! copies a two dimensional buffer onto a two dimensional buffer
  !-----------------------------------------------------------------------!
        subroutine copy_buffer_2_h(buffer_in,buffer_out,
       &     dim1,dim2,ioff,transpose )
  !-----------------------------------------------------------------------!
        implicit none
        real(8),intent(in)::
       &     buffer_in(:,:)
        real(8),intent(inout)::
       &     buffer_out(:)
        integer,intent(in)::
       &     dim1,dim2,
       &     ioff
        logical,intent(in)::
       &     transpose
  
        integer::
       &     idx,jdx
  
        if(transpose)then
           do idx = 1,dim1
              do jdx = 1,dim2
                 buffer_out(ioff+(idx-1)*dim2+jdx)  =  buffer_in(jdx,idx)
              end do
           end do
        else
           do idx = 1,dim1
              do jdx = 1,dim2
                 buffer_out(ioff+(idx-1)*dim2+jdx)  =  buffer_in(idx,jdx)
              end do
           end do
        end if
        end subroutine
  !-----------------------------------------------------------------------!
  !>    returns the number of distributions within a symmetry /MS block
  !!
  !!    @param mel the me-list in question
  !!    @param occupation block number
  !!    @param igam symetry number (of anihilation symmetry)
  !!    @param ms index      (of anihilation MS
  !-----------------------------------------------------------------------!
        pure function mel_get_ndis_h(mel,iocc,igam,idxms)
  !-----------------------------------------------------------------------!
        implicit none
  
        integer::
       &     mel_get_ndis_h
  
        type(me_list),intent(in)::
       &     mel
        integer,intent(in)::
       &     iocc,                !occupation block number
       &     igam,                !symetry number (of anihilation symmetry
       &     idxms                ! ms index      (of anihilation MS
  
        mel_get_ndis_h=
       &     mel%off_op_gmox(iocc)%ndis(igam,idxms) !yep, it is saved there.
        end function
  
  !-----------------------------------------------------------------------!
  !>    returns length of a  symmetry /MS block
  !!
  !!    @param mel the me-list in question
  !!    @param occupation block number
  !!    @param igam symetry number (of anihilation symmetry)
  !!    @param ms index  of anihilation MS
  !-----------------------------------------------------------------------!
        pure function  mel_get_lenblock_h(mel,iocc,igam,idxms)
  !-----------------------------------------------------------------------!
        implicit none
  
        integer::
       &     mel_get_lenblock_h
  
        type(me_list),intent(in)::
       &     mel
        integer,intent(in)::
       &     iocc,
       &     igam,
       &     idxms
        mel_get_lenblock_h =
       &     mel%len_op_gmo(iocc_cls)%gam_ms(igama,idxmsa)
        end function mel_get_lenblock_h
  
  !-----------------------------------------------------------------------!
  !>    finds the square root of an integer (delivered as integer)
  !!
  !!     undefined if a negative integer is passed
  !!    @param in integer
  !-----------------------------------------------------------------------!
        function isqrt_h(in)
  !-----------------------------------------------------------------------!
        implicit none
  
        integer::
       &     isqrt_h
        integer,intent(in)::
       &     in
        integer ::
       &     last_step
  
        isqrt_h = int(sqrt(dble(in)))
  
        if(isqrt_h*isqrt_h.ne.in) call quit(1,i_am,"square root error")
        end function isqrt_h
  !----------------------------------------------------------------------!
        subroutine expand_len_str_h(len_str_out, len_str,
       &     nablk,ncblk,
       &     nc1,nc2,
       &     na1,na2)
  !-----------------------------------------------------------------------!
        implicit none
        integer, intent(inout)::
       &     len_str_out(4)
        integer, intent(in)::
       &     len_str(*),
       &     nablk,ncblk,
       &     nc1,nc2,
       &     na1,na2
        len_str_out(1:4) = 1
        if (ncblk.eq.1.and.nc2.eq.0) len_str_out(1) = len_str(1)
        if (ncblk.eq.1.and.nc1.eq.0) len_str_out(2) = len_str(1)
        if (ncblk.eq.2) then
           len_str_out(1) = len_str(1)
           len_str_out(2) = len_str(2)
        end if
        if (nablk.eq.1.and.na2.eq.0) len_str_out(3) =len_str(ncblk+1)
        if (nablk.eq.1.and.na1.eq.0) len_str_out(4) =len_str(ncblk+1)
        if (nablk.eq.2) then
           len_str_out(3) = len_str(ncblk+1)
           len_str_out(4) = len_str(ncblk+2)
        end if
  
        end subroutine
  !-----------------------------------------------------------------------!
        subroutine calculate_rank_dimensions_h(
       &     mel,op,graphs,iocc_cls, njoined,
       &     ms1, igam,
       &     sgrm,project,
       &     na1mx,nc1mx,
       &     ndim,nrank,rankoff,rankdim)
  !-----------------------------------------------------------------------!
        implicit none
        type(me_list),intent(in)::
       &     mel
        type(operator),intent(in)::
       &     op
        type(graph), pointer,intent(in)::
       &     graphs(:)
        integer,intent(in)::
       &     na1mx,nc1mx,         !maximum number of anihilation/creation operators in vertex pair1
       &     iocc_cls,njoined,
       &     ms1, igam,           !ms and symmetry filter
       &     project
        logical, intent(in)::
       &     sgrm                 ! if separate orthogonalization is requested
        integer, intent(inout)::
       &     ndim, nrank
        integer, dimension(maxrank),intent(out)::
       &     rankoff, rankdim
        integer::
       &     jocc_cls, jblkoff,   ! running indices for occupation class number and offset
       &     na1,nc1,
       &     na2,nc2,
       &     rdim
  
        ndim = 0
        nrank = 0
  
        jocc_cls = iocc_cls
        jblkoff = occoff_h(jocc_cls,njoined)
  
        na1 = na1mx
        nc1 = nc1mx
        rankoff = 0
        rankdim = 0
        blk_loop: do while(min(na1,nc1).ge.0.and.na1+nc1.ge.abs(ms1))
             na2 = nc1mx
             nc2 = na1mx
             rdim = 0
  
             do while(min(na2,nc2).ge.0.and.na2+nc2.ge.abs(ms1))
  
             ! no off-diagonal blocks in SepO (separate orthog.)
             if (sgrm.and.na1.ne.nc2) then
                na2 = na2 - 1
                nc2 = nc2 - 1
                cycle
             end if
  
              ! exit if there is no block like this
             if ( na1.ne.sum( op_inp%ihpvca_occ(1:ngastp,2,jblkoff+1)).or.
       &          nc1.ne.sum( op_inp%ihpvca_occ(1:ngastp,1,jblkoff+2)).or.
       &          na2.ne.sum( op_inp%ihpvca_occ(1:ngastp,2,jblkoff+2)).or.
       &          nc2.ne.sum( op_inp%ihpvca_occ(1:ngastp,1,jblkoff+3)))
       &          exit blk_loop
  
  ! skip off-diagonal blocks
              if (na1.ne.nc2) then
                 call advance_block_n2_h(nc2,na2,
       &              ms1,nc1mx,na1mx,jocc_cls)
                 jblkoff = occoff_h(jocc_cls,njoined)
                cycle
              end if
              ! fasttrack pure inactive block:
              if (na1.eq.0.and.nc1.eq.0) then
                 if (igam.eq.1)then
                    ndim = ndim + 1
                    call update_rankarrays_h(nrank,rankdim,rankoff,1)
                 end if
                exit blk_loop
             end if
  
  ! ASSUMPTION rdim = 0 at this point
  ! currently ensured by:
  ! because only diagonal blocks reach this point.
  ! so only 1 block per na1,nc1 touple
  ! a new na1 nc1 touple is generated at the end of an blk_loop iteration
  ! at the begin of next iteration rdim is set to zero
             call get_combineddistlen_h(
       &          rdim,
       &          op,mel,
       &          jocc_cls,
       &          ms1,igam,
       &          graphs)
             call advance_block_n2_h(nc2,na2,ms1,nc1mx,na1mx,jocc_cls)
             jblkoff = occoff_h(jocc_cls,njoined)
             end do
  
             na1 = na1 - 1
             nc1 = nc1 - 1
  
             ndim = ndim + rdim
             if (project.gt.0.and.rdim.gt.0) then
                call update_rankarrays_h(nrank,rankdim,rankoff,rdim)
             end if
          end do blk_loop
  
        end subroutine calculate_rank_dimensions_h
  !-----------------------------------------------------------------------!
  !>  adds the number of all strings in the ioff block to len
  !
  !!   @param[inout] len cumulative number of all strings with this ms restrictions
  !!   @param hpvx_occ (occupational numbers for all spaces for all vertices)
  !!   @param idx_graph graph indices for all spaces of all vertices
  !!   @param iocc, index of the occupation block
  !!   @param ms1,gam1 combined ms and symmetry of C1/A1 tuple
  !
  !    Note: only works for derived density matrices ms_total =0, totalsymmetrisch , form:
  !
  !!    / 0 0 0 0 \
  !!    \ 0 0 x 0 /  msc1 gamc1 nc1
  !!    / 0 0 y 0 \  msa1 gama1 na1
  !!    \ 0 0 y 0 /  msc2 gamc2 nc2
  !!    / 0 0 x 0 \  msa2 gama2 na2
  !!    \ 0 0 0 0 /
  !-----------------------------------------------------------------------!
         subroutine get_combineddistlen_h(
       &     len,
       &     op,mel,
       &     iocc,
       &     ms1,gam1,
       &     graphs)
  !-----------------------------------------------------------------------!
        implicit none
        type(me_list),intent(in)::
       &     mel
        type(operator),intent(in)::
       &     op
        integer,intent(in)::
       &     iocc
        integer,intent(in)::
       &     ms1,gam1
        type(graph), pointer,intent(in)::
       &     graphs(:)
  
        integer,intent(inout)::
       &     len
        integer::
       &     ncblk,nablk,ioff,njoined
  
        njoined = op%njoined
  
        ioff = (iocc-1)*njoined
  
        call get_num_subblk(ncblk,nablk,
       &     op%ihpvca_occ(1,1,ioff+1),njoined)
        call get_combineddistlen_core_h(
       &     len,
       &     op%ihpvca_occ,mel%idx_graph,
       &     ioff, njoined,
       &     ncblk,nablk,
       &     ms1,igam,
       &     graphs)
  
        end subroutine
  
  !-----------------------------------------------------------------------!
        subroutine get_combineddistlen_core_h(
       &     len,
       &     hpvx_occ,idx_graph,
       &     ioff,njoined,
       &     ncblk,nablk,
       &     ms1,gam1,
       &     graphs)
  !-----------------------------------------------------------------------!
        implicit none
  
        integer::
       &     hpvx_occ(ngastp,2,*),         !ngastp from hostinclude
       &     idx_graph(ngastp,2,*)
        integer,intent(in)::
       &     ioff,                !offset of block
       &     njoined,
       &     nablk,ncblk
  
        integer,intent(in)::
       &     ms1,gam1
        type(graph), pointer,intent(in)::
       &     graphs(:)
        integer,intent(inout)::
       &     len
        integer,pointer::
       &     len_str(:)
        integer::
       &     hpvx_csub(ncblk),hpvx_asub(nablk),
       &     occ_csub(ncblk),occ_asub(nablk),
       &     graph_csub(ncblk), graph_asub(nablk)
        integer ::
       &     gamdis_c(ncblk),gamdis_a(nablk),
       &     msdis_c(ncblk),msdis_a(nablk)
        integer ::
       &     idxmsdis_c(ncblk),idxmsdis_a(nablk) !index of
        integer::
       &     msc1,gama1,gamc1
        type(gam_generator_t)::
       &     gama_gen
        type(ms_generator_t)::
       &     msa_gen
  
  
        allocate(len_str(nablk+ncblk))
  
        na1=sum(hpvx_occ(1:ngastp,2,ioff+1))
        nc1=sum(hpvx_occ(1:ngastp,1,ioff+2))
        call condense_occ(occ_csub, occ_asub,
       &     hpvx_csub,hpvx_asub,
       &     hpvx_occ(1,1,ioff+1),njoined,hpvxblkseq) !hpvxblkseq from hostinclude
  
        call condense_occ(graph_csub, graph_asub,
       &     hpvx_csub,hpvx_asub,
       &     idx_graph(1,1,ioff+1),njoined,hpvxblkseq)
  
        call init_ms_generator_h(msa_gen,na1,ms1)
        do while(next_ms_h(msa_gen,msa1,msc1) )
           if (.not. is_possible_h(nc1, msc1,1)) cycle
  
           call init_gam_generator_h(gama_gen,ngam,igam)
           do while(next_gam_h(gama_gen,gama1,gamc1) )
              if (.not. is_possible_h(na1,0,gama1) ) cycle
              if (.not. is_possible_h(nc1,0,gamc1) ) cycle
  
              call set_dist_msgam_h(
       &           gamdis_c,gamdis_a,
       &           msdis_c,msdis_a,
       &           msc1,msa1,
       &           msa1,msc1,
       &           gamc1,gama1,
       &           gama1,gamc1,
       &           nc1,na1,
       &           na1,nc1,
       &           ncblk,nablk)
  
  
              call ms2idxms(idxmsdis_c,msdis_c,occ_csub,ncblk) !translates ms numbers on msdis into indexes
              call ms2idxms(idxmsdis_a,msdis_a,occ_asub,nablk)
  
              call set_len_str(len_str,ncblk,nablk, ! finds the number of strings for for the given distribution
       &           graphs,
       &           graph_csub,idxmsdis_c,gamdis_c,hpvx_csub,
       &           graph_asub,idxmsdis_a,gamdis_a,hpvx_asub,
       &           hpvxseq,.false.)
              if (nablk.eq.1) then
                 len = len + len_str(1)
              else
                 len = len + len_str(1)*len_str(3)
              end if
           end do
        end do
        deallocate(len_str)
        end subroutine
  
  
  !-----------------------------------------------------------------------!
  !!
  !!  sets the ms and gam values for the distributions
  !!
  !!   @param[out] gamdis_c,gamdis_a symmetry of the creator/anihilator strings of a specific distribution (see also condense_occ.f)
  !!   @param[out] msmdis_c,msmdis_a symmetry of the creator/anihilator strings of a specific distribution
  !!    / 0 0 0 0 \
  !!    \ 0 0 x 0 /  msc1 gamc1 nc1
  !!    / 0 0 y 0 \  msa1 gama1 na1
  !!    \ 0 0 y 0 /  msc2 gamc2 nc2
  !!    / 0 0 x 0 \  msa2 gama2 na2
  !!    \ 0 0 0 0 /
  !-----------------------------------------------------------------------!
        subroutine set_dist_msgam_h(
       &     gamdis_c,gamdis_a,
       &     msdis_c,msdis_a,
       &     msc1,msa1,
       &     msc2,msa2,
       &     gamc1,gama1,
       &     gamc2,gama2,
       &     nc1,na1,
       &     nc2,na2,
       &     ncblk,nablk)
  !-----------------------------------------------------------------------!
        implicit none
        integer,intent(inout)::
       &     gamdis_c(*),gamdis_a(*),
       &     msdis_c(*),msdis_a(*)
        integer,intent(in)::
       &     msc1,msa1,
       &     msc2,msa2,
       &     gamc1,gama1,
       &     gamc2,gama2,
       &     nc1,na1,
       &     nc2,na2,
       &     ncblk,nablk
  
        msdis_c(1) = msc1
        msdis_a(1) = msa1
        gamdis_c(1) = gamc1
        gamdis_a(1) = gama1
  
         if(nablk.eq.1 .and. na1.eq.0)then
            msdis_a(1) = msa2
            gamdis_a(1) = gama2
         else if (nablk.eq.2)then
            msdis_a(2) = msa2
            gamdis_a(2) = gama2
         end if
  
        if(ncblk.eq.1 .and. nc1.eq.0)then
           msdis_c(1) = msc2
           gamdis_c(1) = gamc2
        else if (ncblk.eq.2)then
           msdis_c(2) = msc2
           gamdis_c(2) = gamc2
        end if
  
        end subroutine
  !------------------------------------------------------------------!
  !!    calculates the occupation offset of an occupation on the hpvx_occ array
  !!
  !!   @param iocc index of the occupation
  !!   @param njoined numbers of vertexes the operator has
  !------------------------------------------------------------------!
        pure function occoff_h(iocc,njoined)
  !------------------------------------------------------------------!
        implicit none
        integer::
       &     occoff_h
        integer,intent(in)::
       &     iocc,njoined
        occoff_h = (iocc-1)*njoined
        end function
  !------------------------------------------------------------------!
  !!
  !!    generates the next block by reducing nc2 and na2:
  !!    ,v;v,vv;vv, -> ,v;v,v;v,
  !------------------------------------------------------------------!
        subroutine advance_block_n2_h(
       &     nc2,na2,
       &     ms1,
       &     nc1mx,na1mx,
       &     jocc_cls)
  !------------------------------------------------------------------!
        implicit none
        integer,intent(in)::
       &     ms1,
       &     nc1mx,na1mx
  
        integer,intent(inout)::
       &     nc2,na2,
       &     jocc_cls
  
        na2 = na2 - 1
        nc2 = nc2 - 1
        jocc_cls = jocc_cls + 1
        if (na2+nc2.lt.abs(ms1))  !jump to next line
       &     jocc_cls = jocc_cls + min(na1mx,nc1mx)
       &     -(na1mx+nc1mx-abs(ms1))/2
        end subroutine
  !------------------------------------------------------------------!
  !>
  !!    so the dimension for an excitation rank lies on rankdim(nrank+1 )
  !!    because then the information for ,;,;, lies on rankdim(1)
  !!    on rankoff lies the offset for the respective rank.
  !------------------------------------------------------------------!
        subroutine update_rankarrays_h(nrank,rankdim,rankoff,rdim)
  !------------------------------------------------------------------!
        implicit none
        integer,intent(inout)::
       &     nrank,
       &     rankdim(maxrank),
       &     rankoff(maxrank)
        integer,intent(in)::
       &     rdim
        integer::
       &     idx
  
        nrank = nrank+1
        if(nrank.gt.maxrank)call quit(1,i_am,
       &     'increase maxrank')
  
        do idx=nrank,2,-1
           rankdim(idx)=rankdim(idx-1)
           rankoff(idx)=rankoff(idx-1)
        end do
        rankdim(1)=rdim
        rankoff(1)=0
        if(nrank.ge.2) rankoff(1)=rankoff(2)+rankdim(2)
  !      print *, nrank,rdim
  !      print *, rankoff
  !      print *, rankdim
        end subroutine
  
  !###################################################################!
  !!  Generators: They iterate through a number of possible states:
  !!
  !!  usage:
  !!  call init_generator(mygenerator, ....   )
  !!  do  while (next(mygenerator, ...values) ) ! values are intent(out) or intent(inout)
  !!     ! code, that uses values
  !!  end do
  !!  call del_generator(mygenerator) ! mygenerator may now be reused (new call to init_generator...)
  !
  
  !##################################################################!
  !! ms Generator: generates all combinations  of creator ms values and anihilator ms values possible for a given number of electrons
  !!  Attention: ms values correspond to 2*physical ms (here and everywhere else) to avoid using floating point units
  
  !------------------------------------------------------------------!
  !!   initializes an ms Generator
  !!
  !!
  !!   @param ms_gen the ms_generator
  !!   @param nel number of electrons
  !!   @param total_ms total ms
  !------------------------------------------------------------------!
        subroutine init_ms_generator_h(ms_gen,nel,total_ms)
  !------------------------------------------------------------------!
        implicit none
        type(ms_generator_t),intent(inout)::
       &     ms_gen
        integer,intent(in)::
       &     nel,total_ms
  
        ms_gen%last_ms = nel+2 !msa starts at last_ms -2 so first value is msa=nel
        ms_gen%nel = nel
        ms_gen%total_ms = total_ms
  
        end subroutine
  
  
  !------------------------------------------------------------------!
        subroutine del_ms_generator_h(ms_gen)
  !------------------------------------------------------------------!
        implicit none
        type(ms_generator_t),intent(inout)::
       &     ms_gen
        continue !noting to do
        end subroutine
  !------------------------------------------------------------------!
  !!  generates actual ms values
  !!
  !!  note that msa and msc should not be used if the function returns false
  !!  @param msa the ms value of the anihilator string
  !!  @param msc the ms value of the creator string
  !------------------------------------------------------------------!
        function next_ms_h(ms_gen,msa,msc)
  !------------------------------------------------------------------!
        implicit none
        logical::
       &     next_ms_h
        type(ms_generator_t),intent(inout)::
       &     ms_gen
        integer,intent(inout)::
       &     msa,msc
  
        if (ms_gen%last_ms .le. -ms_gen%nel)then
           ! setting some defined values
           msa=0
           ms_gen%last_ms = msa
           msc = msa + ms_gen%total_ms
           next_ms_h =.false.
        else
           msa= ms_gen%last_ms-2
           ms_gen%last_ms = msa
           msc = msa + ms_gen%total_ms
           next_ms_h =.true.
        end if
  
        end function
  
  
  
  
  !##################################################################!
  !! gam Generator: generates all combinations  of creator symmetry indexes and anihilator indexes possible for a given number of electrons
  
  !------------------------------------------------------------------!
        subroutine init_gam_generator_h(gam_gen,ngam,total_gam)
  !------------------------------------------------------------------!
        implicit none
        type(gam_generator_t),intent(inout)::
       &     gam_gen
        integer,intent(in)::
       &     ngam,total_gam
  
        gam_gen%last_gam=0 !gama starts at last_gam+1
        gam_gen%ngam = ngam
        gam_gen%total_gam = total_gam
        end subroutine
  
  
  !------------------------------------------------------------------!
        subroutine del_gam_generator_h(gam_gen)
  !------------------------------------------------------------------!
        implicit none
        type(gam_generator_t),intent(inout)::
       &     gam_gen
        continue !noting to do
        end subroutine
  
  !------------------------------------------------------------------!
  !!
  !!
  !!  note that gama and gamc should not be used if the function returns false
  !------------------------------------------------------------------!
  !!  generates actual ms values
  !!
  !!  note that gama and gamc should not be used if the function returns false
  !!  @param gama the symmetry number of the anihilator string
  !!  @param gamc the symmetry number of the creator string
  !------------------------------------------------------------------!
        function next_gam_h(gam_gen,gama,gamc)
  !------------------------------------------------------------------!
        implicit none
        logical::
       &     next_gam_h
        integer,intent(inout)::
       &     gama,gamc
        type(gam_generator_t),intent(inout)::
       &     gam_gen
  
        if (gam_gen%last_gam .ge. gam_gen%ngam)then
           ! setting some defined values
           gama=0
           gam_gen%last_gam = gama
           gamc = multd2h(gama,gam_gen%total_gam)
           next_gam_h =.false.
        else
           gama= gam_gen%last_gam+1
           gam_gen%last_gam = gama
           gamc = multd2h(gama,gam_gen%total_gam)
           next_gam_h =.true.
        end if
        end function
  !##################################################################!
  
  
  !------------------------------------------------------------------!
  !!     this could be easier done by just transferring mat(a_start:a_end,b_start:b_end)
  !!     where mat is used, but ifort then allocates the matrix on the stack => possible stackoverflow
  !!     instead we explicitly allocate them on the heap
  !!
  !!     I miss C++
  !------------------------------------------------------------------!
        subroutine extract_submatrix_h(mat,
       &     dim1, a_start, a_end,
       &     dim2, b_start, b_end,
       &     outmat)
  !------------------------------------------------------------------!
        implicit none
        integer,intent(in)::
       &     dim1, a_start, a_end,
       &     dim2, b_start, b_end
  
        double precision, intent(in)::
       &     mat(dim1,dim2)
        double precision, pointer::
       &     outmat(:,:)
        allocate(outmat(a_end-a_start+1,b_end-b_start+1))
        outmat = mat(a_start:a_end,b_start:b_end)
        end subroutine
  
  !------------------------------------------------------------------!
        subroutine insert_submatrix_h(
       &     mat,
       &     dim1, a_start, a_end,
       &     dim2, b_start, b_end,
       &     inmat)
  !------------------------------------------------------------------!
        implicit none
        integer,intent(in)::
       &     dim1, a_start, a_end,
       &     dim2, b_start, b_end
  
        double precision, pointer::
       &     inmat(:,:)
        double precision::
       &     mat(dim1,dim2)
        mat(a_start:a_end,b_start:b_end) = inmat
        deallocate(inmat)
        inmat=> null()
        end subroutine
  !------------------------------------------------------------------!
        subroutine create_submatrix_h(
       &     a_start, a_end,
       &     b_start, b_end,
       &     outmat)
  !------------------------------------------------------------------!
        implicit none
        integer,intent(in)::
       &     a_start, a_end,
       &     b_start, b_end
        double precision, pointer::
       &     outmat(:,:)
        allocate(outmat(a_end-a_start+1,b_end-b_start+1))
        end subroutine
  !------------------------------------------------------------------!
  !------------------------------------------------------------------!
        subroutine delete_submatrix_h(
       &     inmat)
  !------------------------------------------------------------------!
        implicit none
        double precision, pointer::
       &     inmat(:,:)
        deallocate(inmat)
        inmat=>null()
        end subroutine
  !------------------------------------------------------------------!
  
  !------------------------------------------------------------------!
  !!    checks if a given combination of electron number, ms and symmetry is possible for a gam/ms block
  !!    ms=0 and gam =1 are qlways possible and check only the other condition
  !------------------------------------------------------------------!
        pure function is_possible_h(nel, ms,gam)
  !------------------------------------------------------------------!
        implicit none
        logical::
       &     is_possible_h
  
        integer,intent(in)::
       &     nel,ms,gam
  
        is_possible_h = .true.
        if (abs(ms).gt. nel)is_possible_h = .false.
        if (nel .eq.0 .and. gam.ne.1)is_possible_h = .false.
  
        end function
  !------------------------------------------------------------------!
  ! we need to distinguish:
  ! /0 0 0 0\     /0 0 0 0\
  ! \0 0 x 0/     \0 0 0 0/ i.e. if creators come first,
  ! /0 0 0 0\ and /0 0 x 0\ we need to transpose the
  ! \0 0 0 0/     \0 0 x 0/ scratch matrix
  ! /0 0 x 0\     /0 0 0 0\ (important if result is non-symmetric!)
  ! \0 0 0 0/     \0 0 0 0/
  !------------------------------------------------------------------!
        function use_transposed_scratch_h(hpvx_occ, njoined)
        character(len=*),parameter ::
       &     i_am = 'invsqrt:use_transposed_scratch_h'
        logical::
       &     use_transposed_scratch_h
  
        integer, intent(in)::
       &     njoined,
       &     hpvx_occ(ngastp,2,njoined)
        logical ::
       &     error
        integer ::
       &     ij
  
        if ( njoined .ne. 3) then
          call quit(2, i_am, "not prepared for njonied !=3")
        end if
        use_transposed_scratch_h = .false. ! necessary to give an result if there are no creators or anihilators
        do ij = 1, njoined
          if (sum(hpvx_occ(1:ngastp,1,iblkoff+ij)).gt.0) then
            use_transposed_scratch_h = .true.
            exit
          else if (sum(hpvx_occ(1:ngastp,2,iblkoff+ij)).gt.0) then
            use_transposed_scratch_h = .false.
            exit
          end if
        end do
        return
        end function
  !------------------------------------------------------------------!
        end
  