text:
! -*- mode: F90; mode: font-lock -*-
! ------------------------------------------------------------------------------
! $Id$
! ------------------------------------------------------------------------------
! Module cover_module
! ------------------------------------------------------------------------------
! Code area 8: indexing
! -------------------------------------------------------------------------------

!!****h* Conquest/cover_module
!!  NAME
!!   cover_module
!!  PURPOSE
!!   This module holds variables and routines related to covering
!!   sets.  A covering set is (at the moment) the smallest orthorhombic
!!   collection of small groups that will fully enclose a primary set
!!   out to a specified radius.  It's used for searching for neighbours
!!   and communicating details about small groups between processors.
!!
!!   More details (exhaustive details !) can be found in the Conquest
!!   notes "Matrix multiplication in Conquest: Practical details"
!!   which are in the Conquest documentation repository (stored
!!   in the directory TechnicalMatMult/)
!!  USES
!!   GenComms, basic_types, datatypes, global_module
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99
!!  MODIFICATION HISTORY
!!   19/02/00 by D.R.Bowler
!!    Implementing changes for transposes (adding ncover_rem and
!!    routines to find iprim_part)
!!   03/03/00 by DRB
!!    Finishing up transposes
!!   29/03/00 by DRB
!!    Creating a generalised cover type and altering make_gcs to use it
!!   11/04/00 by DRB
!!    Correcting cover type
!!   19/04/00 by DRB
!!    Adding general group and primary set types to make GCS totally
!!    general
!!   20/04/00 by DRB
!!    Removed type declaration and placed in basic_types module
!!    Continued generalisation of make_cs
!!    Removed GOTOs from indexx (replaced with do...while)
!!   21/06/2001 dave
!!    Added ROBODoc headers and removed stop statements and MPI_Abort
!!    statements throughout
!!   19/04/2002 drb
!!    Corrected reference to notes and added USE field for whole module
!!   17/06/2002 dave
!!    Added RCS Id tag and improved headers a little
!!   2008/02/04 17:12 dave
!!    Changes for output to file not stdoutx
!!   2008/05/16 ast
!!    Added some timers
!!   2014/09/15 18:30 lat
!!    fixed call start/stop_timer to timer_module (not timer_stdlocks_module !)
!!   2016/01/28 16:42 dave
!!    Changed ewald_CS to ion_ion_CS
!!   2019/11/04 15:12 dave
!!    Replace indexx from Numerical Recipes with generic heapsort
!!  SOURCE
!!
module cover_module

  ! Module use
  use datatypes
  use basic_types
  use global_module,          only: io_lun
  use timer_module,           only: start_timer, stop_timer
  use timer_stdclocks_module, only: tmr_std_indexing, tmr_std_allocation

  implicit none

  save

  ! DCS = domain CS; BCS = bundle CS
  type(cover_set) :: DCS_blocks
  type(cover_set) :: DCS_parts
  type(cover_set) :: BCS_blocks
  type(cover_set) :: BCS_parts
  type(cover_set) :: ion_ion_CS
  type(cover_set) :: D2_CS ! for DFT-D2

!!***
contains

!!****f* cover_module/make_cs *
!!
!!  NAME
!!   make_cs - makes a covering set
!!  USAGE
!!
!!  PURPOSE
!!   sbrt make_cs: create a covering set out of groups given by the
!!   group_set variable, for a cutoff gr_rcut around the primary set
!!   defined by the widths (nw_x,y,z) and left-spans (nlx,y,z).
!!   N.B. The widths and left spans must have been converted into units
!!   of the CS groups !
!!   The origin is also passed; this is a problem is we're making a
!!   CS for one group type out of the other.  It MUST be the group
!!   nearest the origin partition of the primary set.
!!   The results are put into the derived type set.
!!   The logical variable members decides whether or not the
!!   x,y,zcover variables are created (NOT for the cross-group CS)
!!
!!   The notes above are a little out of date - the widths and left-spans
!!   are now part of the primary set type, and conversion is done by
!!   convert_primary (below)
!!  INPUTS
!!
!!
!!  USES
!!
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99
!!  MODIFICATION HISTORY
!!   21/06/2001 dave
!!    Added ROBODoc header and cq_abort
!!   2008/05/16 ast
!!    Added some timers
!!   2013/07/01 M.Arita
!!    Added allocation of ig_cover
!!   2013/08/21 M.Arita
!!    Initialised and calculated ig_cover for loading matrices
!!  SOURCE
!!
  subroutine make_cs(myid, gr_rcut, set, groups, prim, mx_mcell, &
                     x_mem_cell, y_mem_cell, z_mem_cell)

    ! Module usage
    use global_module
    use basic_types
    use GenComms, ONLY: cq_abort, inode, ionode
    use memory_module, ONLY: reg_alloc_mem, reg_dealloc_mem, type_int,&
         type_dbl
    use functions, ONLY: heapsort_integer_index

    implicit none

    ! Passed variables
    type(cover_set) :: set
    type(group_set), target   :: groups
    type(primary_set), target :: prim
    integer, intent(in) :: myid
    integer, OPTIONAL :: mx_mcell
    real(double), dimension(:), intent(IN), OPTIONAL :: x_mem_cell
    real(double), dimension(:), intent(IN), OPTIONAL :: y_mem_cell
    real(double), dimension(:), intent(IN), OPTIONAL :: z_mem_cell
    real(double) :: gr_rcut ! Cutoff radius for CS

    ! Local variables
    ! Number of periodic images of a group
    integer :: nrepx(groups%mx_gedge)
    integer :: nrepy(groups%mx_gedge)
    integer :: nrepz(groups%mx_gedge)
    ! x,y,z numbering of CS groups (for CC labels)
    integer, allocatable, dimension(:)  :: nx_in_cover
    integer, allocatable, dimension(:)  :: ny_in_cover
    integer, allocatable, dimension(:)  :: nz_in_cover
    ! Variables for irreducible CS
    integer :: ind_min(groups%mx_gcell)
    integer :: ngcx_min(groups%mx_gcell)
    integer :: ngcy_min(groups%mx_gcell)
    integer :: ngcz_min(groups%mx_gcell)
    integer :: min_sort(groups%mx_gcell)
    ! Converted units of primary cell size
    integer :: nx_o,ny_o,nz_o
    ! General scalars
    integer :: ng_in_cell,nmodx,nmody,nmodz
    integer :: noccx,nremx,minx,ngcx,noccy,nremy,miny,ngcy,noccz,nremz
    integer :: minz,ngcz,ng_in_min,ind,nqx,nqy,nqz,ind_qart,ino,ind_cover
    integer :: nrx,nry,nrz,nsx,nsy,nsz,ni,nnd,stat
    integer :: nm_in_cover
    real(double) :: dcellx,dcelly,dcellz,xadd,yadd,zadd
    logical :: members

    call start_timer(tmr_std_indexing)
    ! This determines whether or not we find information about all
    ! the members of the groups in the covering set
    members = (PRESENT(mx_mcell).AND.PRESENT(x_mem_cell).AND.&
         PRESENT(y_mem_cell).AND.PRESENT(z_mem_cell))
    nnd = myid+1
    ! Start the derived type - assign useful variables
    set%rcut = gr_rcut     ! Cutoff around primary set
    set%groups => groups   ! Groups out of which to build CS
    set%prim => prim       ! Primary set around which to build CS
    ! Find the origin, width and left span of the CS
    call convert_primary(groups,prim,set,nx_o,ny_o,nz_o)
    ! This gives the total number of groups
    set%ng_cover=set%ncoverx*set%ncovery*set%ncoverz
    set%mx_gcover = set%ng_cover
    call allocate_cs(set,numprocs,members)
    call start_timer(tmr_std_allocation)
    allocate(nx_in_cover(set%ng_cover),ny_in_cover(set%ng_cover),nz_in_cover(set%ng_cover),STAT=stat)
    if(stat/=0) call cq_abort("Error allocating nx_in_cover: ",set%ng_cover,stat)
    call reg_alloc_mem(area_index, 3*set%ng_cover,type_int)
    call stop_timer(tmr_std_allocation)
    !if(set%ng_cover>set%mx_gcover) then
    !   call cq_abort('make_cs: too many groups in GCS',&
    !        set%ng_cover,set%mx_gcover)
    !endif
    ! Conversion factors from unit cell lengths->groups
    dcellx=rcellx/real(groups%ngcellx,double)
    dcelly=rcelly/real(groups%ngcelly,double)
    dcellz=rcellz/real(groups%ngcellz,double)
    ! Fully explained in notes mentioned above
    nmodx=((groups%ngcellx+set%nspanlx-1)/groups%ngcellx)*groups%ngcellx
    nmody=((groups%ngcelly+set%nspanly-1)/groups%ngcelly)*groups%ngcelly
    nmodz=((groups%ngcellz+set%nspanlz-1)/groups%ngcellz)*groups%ngcellz
    ! create apparatus for periodic boundary conditions ---------------
    ! ... x-direction
    noccx=set%ncoverx/groups%ngcellx
    nremx=set%ncoverx-noccx*groups%ngcellx
    minx=min(set%ncoverx,groups%ngcellx)
    if(minx>groups%mx_gedge) then
       call cq_abort('make_cs: too many groups in x-edge',minx)
    endif
    do ngcx=1,minx
       if(ngcx<=nremx) then
          nrepx(ngcx)=noccx+1
       else
          nrepx(ngcx)=noccx
       endif
    enddo
    ! ... y-direction
    noccy=set%ncovery/groups%ngcelly
    nremy=set%ncovery-noccy*groups%ngcelly
    miny=min(set%ncovery,groups%ngcelly)
    if(miny>groups%mx_gedge) then
       call cq_abort('make_cs: too many groups in y-edge',miny)
    endif
    do ngcy=1,miny
       if(ngcy<=nremy) then
          nrepy(ngcy)=noccy+1
       else
          nrepy(ngcy)=noccy
       endif
    enddo
    ! ... z-direction
    noccz=set%ncoverz/groups%ngcellz
    nremz=set%ncoverz-noccz*groups%ngcellz
    minz=min(set%ncoverz,groups%ngcellz)
    if(minz>groups%mx_gedge) then
       call cq_abort('make_cs: too many groups in z-edge',minz)
    endif
    do ngcz=1,minz
       if(ngcz<=nremz) then
          nrepz(ngcz)=noccz+1
       else
          nrepz(ngcz)=noccz
       endif
    enddo
    ! go over groups in GCS periodic-irreducible set, calculating
    ! simulation-cell (node-order, home-start) label of each
    ng_in_cell = groups%ngcellx*groups%ngcelly*groups%ngcellz
    ng_in_min = minx*miny*minz
    ind=0
    do ngcx=1,minx
       do ngcy=1,miny
          do ngcz=1,minz
             ind=ind+1
             nqx=1+mod(nx_o+ngcx-set%nspanlx-2+nmodx,groups%ngcellx)
             nqy=1+mod(ny_o+ngcy-set%nspanly-2+nmody,groups%ngcelly)
             nqz=1+mod(nz_o+ngcz-set%nspanlz-2+nmodz,groups%ngcellz)
             ind_qart=(nqx-1)*groups%ngcelly*groups%ngcellz+(nqy-1)*groups%ngcellz+nqz
             ino=groups%inv_ngnode(ind_qart)
             ind_min(ind)=1+mod(ino-groups%inode_beg(nnd)+ ng_in_cell,ng_in_cell)
             ngcx_min(ind)=ngcx
             ngcy_min(ind)=ngcy
             ngcz_min(ind)=ngcz
          enddo
       enddo
    enddo
    ! sort minimum CS by nodes
    call heapsort_integer_index(ng_in_min,ind_min,min_sort)

    ! go over all GCS groups in node-periodic-grouped order
    ind_cover=0
    do ind=1,ng_in_min
       ngcx=ngcx_min(min_sort(ind))
       ngcy=ngcy_min(min_sort(ind))
       ngcz=ngcz_min(min_sort(ind))
       do nrx=1,nrepx(ngcx)
          do nry=1,nrepy(ngcy)
             do nrz=1,nrepz(ngcz)
                ind_cover=ind_cover+1
                nx_in_cover(ind_cover)=ngcx-1-set%nspanlx+(nrx-1)*groups%ngcellx
                ny_in_cover(ind_cover)=ngcy-1-set%nspanly+(nry-1)*groups%ngcelly
                nz_in_cover(ind_cover)=ngcz-1-set%nspanlz+(nrz-1)*groups%ngcellz
             enddo
          enddo
       enddo
    enddo
    ! make covering set for primary set
    !if(members) set%icover_ibeg(1)=1
    if(members) nm_in_cover = 0
    do ind_cover=1,set%ng_cover
       nsx=nx_in_cover(ind_cover)
       nsy=ny_in_cover(ind_cover)
       nsz=nz_in_cover(ind_cover)
       nqx=1+mod(nx_o+nsx+nmodx-1,groups%ngcellx)
       nqy=1+mod(ny_o+nsy+nmody-1,groups%ngcelly)
       nqz=1+mod(nz_o+nsz+nmodz-1,groups%ngcellz)
       xadd=(nx_o+nsx-nqx)*dcellx
       yadd=(ny_o+nsy-nqy)*dcelly
       zadd=(nz_o+nsz-nqz)*dcellz
       ind_qart=(nqx-1)*groups%ngcelly*groups%ngcellz+&
            (nqy-1)*groups%ngcellz+nqz
       set%lab_cell(ind_cover)=ind_qart
       set%lab_cover(ind_cover)=&
            (nsx+set%nspanlx)*set%ncovery*set%ncoverz+&
            (nsy+set%nspanly)*set%ncoverz+(nsz+set%nspanlz)+1
       !!  TM  temporary ????
       set%inv_lab_cover(set%lab_cover(ind_cover))=ind_cover
       !!  TM  temporary ????
       if(members) nm_in_cover = nm_in_cover + groups%nm_group(ind_qart)
       !if(members.AND.groups%nm_group(ind_qart)>0) then
       !   set%n_ing_cover(ind_cover)=groups%nm_group(ind_qart)
       !   if(set%icover_ibeg(ind_cover)+groups%nm_group(ind_qart)-1> &
       !        set%mx_mcover) then
       !      call cq_abort('make_cs: xcover dimension exceeded',&
       !           set%icover_ibeg(ind_cover)+groups%nm_group(ind_qart)-1, &
       !           set%mx_mcover)
       !   endif
       !   if(groups%icell_beg(ind_qart)+groups%nm_group(ind_qart)-1> &
       !        mx_mcell) then
       !      call cq_abort('make_cs: x_atom_cover dim. exceeded', &
       !           groups%icell_beg(ind_qart)+groups%nm_group(ind_qart)-1)
       !   endif
       !   do ni=1,groups%nm_group(ind_qart)
       !      set%xcover(set%icover_ibeg(ind_cover)+ni-1)= &
       !           x_mem_cell(groups%icell_beg(ind_qart)+ni-1)+xadd
       !      set%ycover(set%icover_ibeg(ind_cover)+ni-1)= &
       !           y_mem_cell(groups%icell_beg(ind_qart)+ni-1)+yadd
       !      set%zcover(set%icover_ibeg(ind_cover)+ni-1)= &
       !           z_mem_cell(groups%icell_beg(ind_qart)+ni-1)+zadd
       !   enddo
       !   if(ind_cover.lt.set%ng_cover) then
       !      set%icover_ibeg(ind_cover+1)=set%icover_ibeg(ind_cover)+ &
       !           groups%nm_group(ind_qart)
       !   endif
       !else if(members) then
       !   if(ind_cover < set%ng_cover) then   !TM 26/Jun/2003
       !      set%icover_ibeg(ind_cover+1)=set%icover_ibeg(ind_cover)+ &
       !           groups%nm_group(ind_qart)
       !   endif  ! (ind_cover < set%ng_cover)  !TM 26/Jun/2003
       !   set%n_ing_cover(ind_cover)=0
       !endif ! members.AND.nm_group>0
    enddo
    if(members) then ! Now generate member information
       if(inode==ionode.AND.iprint_index>1) write(io_lun,*) 'Members in covering set: ',nm_in_cover
       if(members) set%icover_ibeg(1)=1
       set%mx_mcover = nm_in_cover
       call start_timer(tmr_std_allocation)
       allocate(set%xcover(set%mx_mcover),STAT=stat)
       if(stat/=0) call cq_abort('Error allocating cover set members x: ',set%mx_mcover,stat)
       allocate(set%ycover(set%mx_mcover),STAT=stat)
       if(stat/=0) call cq_abort('Error allocating cover set members y: ',set%mx_mcover,stat)
       allocate(set%zcover(set%mx_mcover),STAT=stat)
       if(stat/=0) call cq_abort('Error allocating cover set members z: ',set%mx_mcover,stat)
       allocate(set%ig_cover(set%mx_mcover),STAT=stat)
       if(stat/=0) call cq_abort('Error allocating ig_cover:',set%mx_mcover,stat)
       call reg_alloc_mem(area_index, 3*set%mx_mcover,type_dbl)
       call stop_timer(tmr_std_allocation)
       do ind_cover=1,set%ng_cover
          nsx=nx_in_cover(ind_cover)
          nsy=ny_in_cover(ind_cover)
          nsz=nz_in_cover(ind_cover)
          nqx=1+mod(nx_o+nsx+nmodx-1,groups%ngcellx)
          nqy=1+mod(ny_o+nsy+nmody-1,groups%ngcelly)
          nqz=1+mod(nz_o+nsz+nmodz-1,groups%ngcellz)
          xadd=(nx_o+nsx-nqx)*dcellx
          yadd=(ny_o+nsy-nqy)*dcelly
          zadd=(nz_o+nsz-nqz)*dcellz
          ind_qart= set%lab_cell(ind_cover)
          if(groups%nm_group(ind_qart)>0) then
             set%n_ing_cover(ind_cover)=groups%nm_group(ind_qart)
             if(groups%icell_beg(ind_qart)+groups%nm_group(ind_qart)-1> mx_mcell) then
                call cq_abort('make_cs: x_atom_cover dim. exceeded', &
                     groups%icell_beg(ind_qart)+groups%nm_group(ind_qart)-1)
             endif
             do ni=1,groups%nm_group(ind_qart)
                set%xcover(set%icover_ibeg(ind_cover)+ni-1)= &
                     x_mem_cell(groups%icell_beg(ind_qart)+ni-1)+xadd
                set%ycover(set%icover_ibeg(ind_cover)+ni-1)= &
                     y_mem_cell(groups%icell_beg(ind_qart)+ni-1)+yadd
                set%zcover(set%icover_ibeg(ind_cover)+ni-1)= &
                     z_mem_cell(groups%icell_beg(ind_qart)+ni-1)+zadd
                set%ig_cover(set%icover_ibeg(ind_cover)+ni-1)= &
                     id_glob(groups%icell_beg(ind_qart)+ni-1)
             enddo
             if(ind_cover.lt.set%ng_cover) then
                set%icover_ibeg(ind_cover+1)=set%icover_ibeg(ind_cover)+ &
                     groups%nm_group(ind_qart)
             endif
          else
             if(ind_cover < set%ng_cover) then   !TM 26/Jun/2003
                set%icover_ibeg(ind_cover+1)=set%icover_ibeg(ind_cover)+ &
                     groups%nm_group(ind_qart)
             endif  ! (ind_cover < set%ng_cover)  !TM 26/Jun/2003
             set%n_ing_cover(ind_cover)=0
          endif ! nm_group>0
       enddo
    end if
    !if(members) then ! Added by TM 17/07/00
    !   do np=1,set%ng_cover  ! Loop over partitions in GCS
    !      if(set%n_ing_cover(np).gt.0) then  ! Are there atoms ?
    !         if(set%icover_ibeg(np)+set%n_ing_cover(np)-1.gt.&
    !              set%mx_mcover) then
    !            call cq_abort('make_cs: member index error', &
    !                 set%icover_ibeg(np)+set%n_ing_cover(np)-1)
    !         endif
    !      endif
    !   enddo
    !endif ! (members)
    call start_timer(tmr_std_allocation)
    deallocate(nx_in_cover,ny_in_cover,nz_in_cover,STAT=stat)
    if(stat/=0) call cq_abort("Error deallocating nx_in_cover: ",set%ng_cover,stat)
    call reg_dealloc_mem(area_index, 3*set%ng_cover,type_int)
    call stop_timer(tmr_std_allocation)
    call stop_timer(tmr_std_indexing)
    return
  end subroutine make_cs
!!***

!!****f* cover_module/make_iprim *
!!
!!  NAME
!!   make_iprim
!!  USAGE
!!
!!  PURPOSE
!!   sbrt make_iprim: creates iprim_group, which gives the covering set CC
!!   label for each primary set group - I'm pretty sure that this is
!!   right, but it's a really hard problem
!!
!!   Definitely right - it was broken a while back, but now fixed
!!  INPUTS
!!
!!
!!  USES
!!
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99 ?
!!  MODIFICATION HISTORY
!!   21/06/2001 dave
!!    Added ROBODoc header and cq_abort
!!   2008/05/16 ast
!!    Added timer
!!  SOURCE
!!
  subroutine make_iprim(set,prim)

    use basic_types
    use GenComms, ONLY: cq_abort
    use memory_module, ONLY: reg_alloc_mem, type_int
    use global_module, ONLY: area_index

    implicit none

    ! Passed variables
    type(cover_set) :: set
    type(primary_set) :: prim

    ! Local variables
    integer :: pr,nsx,nsy,nsz,ind_group,i,j,stat

    call start_timer(tmr_std_indexing)
    if(.NOT.ASSOCIATED(set%iprim_group)) then
       call start_timer(tmr_std_allocation)
       allocate(set%iprim_group(prim%mx_iprim),STAT=stat)
       if(stat/=0) call cq_abort('make_iprim: error allocating memory')
       call reg_alloc_mem(area_index,prim%mx_iprim,type_int)
       call stop_timer(tmr_std_allocation)
    endif
    pr = 0 ! Indexes primary members
    do i=1,prim%groups_on_node
       nsx = prim%idisp_primx(i)
       nsy = prim%idisp_primy(i)
       nsz = prim%idisp_primz(i)
       ind_group = (nsx+set%nspanlx)*set%ncovery*set%ncoverz+&
            (nsy+set%nspanly)*set%ncoverz+(nsz+set%nspanlz)+1
       if(prim%nm_nodgroup(i) > 0) then
          do j=1,prim%nm_nodgroup(i)
             pr = pr+1
             if(pr>prim%mx_iprim) call cq_abort('make_iprim: error in pr index',pr)
             set%iprim_group(pr) = ind_group
          enddo
       endif ! (prim%nm_nodgroup(i) > 0) then
    enddo
    call stop_timer(tmr_std_indexing)
  end subroutine make_iprim
!!***

!!****f* cover_module/send_ncover *
!!
!!  NAME
!!   send_ncover
!!  USAGE
!!
!!  PURPOSE
!!   Once make_cs has been run on every processor, this routine
!!   sends sizes of the CS to all other processors (and gathers theirs)
!!  INPUTS
!!
!!
!!  USES
!!
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99 ?
!!  MODIFICATION HISTORY
!!   21/06/2001 dave
!!    ROBODoc header and changed to gsum
!!   2008/08/12 10:49 dave
!!    Tidying: removed unnecessary variable ncovers and allocate/deallocate calls
!!  SOURCE
!!
! =====================================================================
!   sbrt send_ncover: Called once at the start of a run (after make_gcs),
!   to ensure that all nodes know all other nodes GCS sizes (ncoverx)
!   for transpose work
! ---------------------------------------------------------------------
  subroutine send_ncover(set,nnd)

    ! Module declaration
    use global_module, ONLY: numprocs
    use GenComms, ONLY: gsum, cq_abort

    implicit none

    ! Passed variables
    type(cover_set) :: set
    integer :: nnd ! Node no. starting from one

    set%ncover_rem = 0
    set%ncover_rem(3*(nnd-1)+1) = set%ncoverx
    set%ncover_rem(3*(nnd-1)+2) = set%ncovery
    set%ncover_rem(3*(nnd-1)+3) = set%ncoverz
    call gsum(set%ncover_rem, 3*numprocs)
    return
  end subroutine send_ncover
!!***

!!****f* cover_module/convert_primary *
!!
!!  NAME
!!   convert_primary
!!  USAGE
!!
!!  PURPOSE
!!   Given a primary set and a cutoff, creates the widths and left-spans
!!   of the CS in the specified groups and returns the origin of the CS
!!  INPUTS
!!
!!
!!  USES
!!
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99 ?
!!  MODIFICATION HISTORY
!!   21/06/2001 dave
!!    Added ROBODoc header and used cq_abort
!!  SOURCE
!!
  subroutine convert_primary(groups,prim,set,ncx_o,ncy_o,ncz_o)

    use datatypes
    use global_module
    use basic_types
    use numbers, ONLY: very_small
    use GenComms, ONLY: cq_abort

    implicit none

    ! Passed variables
    type(group_set)   :: groups
    type(primary_set) :: prim
    type(cover_set)   :: set
    integer :: ncx_o,ncy_o,ncz_o                ! Origin for cover set

    ! Local variables
    real(double) :: ro_x,ro_y,ro_z              ! Origin in reals
    real(double) :: ro_cx,ro_cy,ro_cz           ! CS origin in reals
    real(double) :: dx,dy,dz,dcx,dcy,dcz        ! Sizes of groups in reals
    real(double) :: lhx,lhy,lhz,lhcx,lhcy,lhcz  ! Left hand corners
    real(double) :: rhx,rhy,rhz,rhcx,rhcy,rhcz  ! Right hand corners
    real(double) :: cslx,csly,cslz              ! Covering set left

    if(.NOT.associated(prim%groups)) then
       call cq_abort('convert_primary: primary set groups pointer not assoc')
    endif
    ! First, find sizes of groups in terms of unit cell
    dx =  rcellx/real(prim%groups%ngcellx,double)   ! Check for a more portable way !
    dy =  rcelly/real(prim%groups%ngcelly,double)
    dz =  rcellz/real(prim%groups%ngcellz,double)
    dcx = rcellx/real(groups%ngcellx,double)
    dcy = rcelly/real(groups%ngcelly,double)
    dcz = rcellz/real(groups%ngcellz,double)
    ! Convert origin of prim to reals (add very_small to prevent ambiguity)
    ro_x = (real(prim%nx_origin-1,double))*dx + very_small
    ro_y = (real(prim%ny_origin-1,double))*dy + very_small
    ro_z = (real(prim%nz_origin-1,double))*dz + very_small
    ! Convert this to number of CS groups
    ro_cx = ro_x/dcx
    ro_cy = ro_y/dcy
    ro_cz = ro_z/dcz
    ! Round and offset (partition 1,1,1 has LH corner at 0.0,0.0,0.0)
    ! Accept implicit type conversion
    ncx_o = 1+anint(ro_cx)
    ncy_o = 1+anint(ro_cy)
    ncz_o = 1+anint(ro_cz)
    set%nx_origin = ncx_o
    set%ny_origin = ncy_o
    set%nz_origin = ncz_o
    ! Find CS origin in reals
    ro_cx = (real(ncx_o-1,double))*dcx
    ro_cy = (real(ncy_o-1,double))*dcy
    ro_cz = (real(ncz_o-1,double))*dcz
    ! Start by finding the left and right hand corners of the CS
    lhx = dx*(prim%nx_origin-1-prim%nleftx)-set%rcut
    lhy = dy*(prim%ny_origin-1-prim%nlefty)-set%rcut
    lhz = dz*(prim%nz_origin-1-prim%nleftz)-set%rcut
    rhx = dx*(prim%nx_origin-1-prim%nleftx+prim%nw_primx)+set%rcut
    rhy = dy*(prim%ny_origin-1-prim%nlefty+prim%nw_primy)+set%rcut
    rhz = dz*(prim%nz_origin-1-prim%nleftz+prim%nw_primz)+set%rcut
    ! The cover set left span must be offset from the CS origin
    cslx = (ro_cx - lhx)/dcx
    csly = (ro_cy - lhy)/dcy
    cslz = (ro_cz - lhz)/dcz
    ! Convert to integers
    set%nspanlx = ceiling(cslx)
    set%nspanly = ceiling(csly)
    set%nspanlz = ceiling(cslz)
    ! The width of the CS must be the distance from the LH corner to the
    ! RH corner - this way we pick the SMALLEST CS needed.
    set%ncoverx = ceiling(rhx/dcx)-(ncx_o-1-set%nspanlx)
    set%ncovery = ceiling(rhy/dcy)-(ncy_o-1-set%nspanly)
    set%ncoverz = ceiling(rhz/dcz)-(ncz_o-1-set%nspanlz)
    ! Calculate the corners of the CS for checks
    lhcx = dcx*(ncx_o-1-set%nspanlx)
    lhcy = dcy*(ncy_o-1-set%nspanly)
    lhcz = dcz*(ncz_o-1-set%nspanlz)
    rhcx = dcx*(ncx_o-1-set%nspanlx+set%ncoverx)
    rhcy = dcy*(ncy_o-1-set%nspanly+set%ncovery)
    rhcz = dcz*(ncz_o-1-set%nspanlz+set%ncoverz)
    ! First check that the LH and RH corners are large enough
    ! LH
    if(lhx+very_small<lhcx) then
       write(io_lun,*) 'CS too small; adjusting nspanlx',lhx,lhcx
       set%nspanlx = set%nspanlx+ceiling((lhcx-lhx)/dcx)
    endif
    if(lhy+very_small<lhcy) then
       write(io_lun,*) 'CS too small; adjusting nspanly',lhy,lhcy
       set%nspanly = set%nspanly+ceiling((lhcy-lhy)/dcy)
    endif
    if(lhz+very_small<lhcz) then
       write(io_lun,*) 'CS too small; adjusting nspanlz',lhz,lhcz
       set%nspanlz = set%nspanlz+ceiling((lhcz-lhz)/dcz)
    endif
    ! RH
    if(rhx>rhcx+very_small) then
       write(io_lun,*) 'CS too small; adjusting ncoverx',rhx,rhcx
       set%ncoverx = set%ncoverx+ceiling((rhx-rhcx)/dcx)
    endif
    if(rhy>rhcy+very_small) then
       write(io_lun,*) 'CS too small; adjusting ncovery',rhy,rhcy
       set%ncovery = set%ncovery+ceiling((rhy-rhcy)/dcy)
    endif
    if(rhz>rhcz+very_small) then
       write(io_lun,*) 'CS too small; adjusting ncoverz',rhz,rhcz
       set%ncoverz = set%ncoverz+ceiling((rhz-rhcz)/dcz)
    endif
    ! Now check that the CS isn't too large
    ! LH
    if(lhx-lhcx>dcx) then
       write(io_lun,*) 'CS too big; adjusting nspanlx',lhx,lhcx
       set%nspanlx = set%nspanlx-floor((lhx-lhcx)/dcx)
    endif
    if(lhy-lhcy>dcy) then
       write(io_lun,*) 'CS too big; adjusting nspanly',lhy,lhcy
       set%nspanly = set%nspanly-floor((lhy-lhcy)/dcy)
    endif
    if(lhz-lhcz>dcz) then
       write(io_lun,*) 'CS too big; adjusting nspanlz',lhz,lhcz
       set%nspanlz = set%nspanlz-floor((lhz-lhcz)/dcz)
    endif
    ! RH
    if(rhcx-rhx>dcx) then
       write(io_lun,*) 'CS too big; adjusting ncoverx',rhx,rhcx
       set%ncoverx = set%ncoverx-floor((rhx-rhcx)/dcx)
    endif
    if(rhcy-rhy>dcy) then
       write(io_lun,*) 'CS too big; adjusting ncovery',rhy,rhcy
       set%ncovery = set%ncovery-floor((rhy-rhcy)/dcy)
    endif
    if(rhcz-rhz>dcz) then
       write(io_lun,*) 'CS too big; adjusting ncoverz',rhz,rhcz
       set%ncoverz = set%ncoverz-floor((rhz-rhcz)/dcz)
    endif
    return
  end subroutine convert_primary
!!***

!!****f* cover_module/allocate_cs *
!!
!!  NAME
!!   allocate_cs
!!  USAGE
!!
!!  PURPOSE
!!   allocates memory to the CS derived type
!!  INPUTS
!!
!!
!!  USES
!!
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99
!!  MODIFICATION HISTORY
!!   21/06/2001 dave
!!    Added ROBODoc header and used cq_abort
!!   2008/05/16 ast
!!    Added timer
!!  SOURCE
!!
  subroutine allocate_cs(set,mx_node,members)

    ! Module usage
    use basic_types
    use GenComms, ONLY: cq_abort
    use global_module, ONLY: area_index
    use memory_module, ONLY: reg_alloc_mem, type_int

    implicit none

    ! Passed variables
    integer :: mx_node
    type(cover_set) :: set
    logical :: members

    ! Local variables
    integer :: stat

    call start_timer(tmr_std_allocation)
    allocate(set%lab_cell(set%mx_gcover),STAT=stat)
    if(stat/=0) then
       call cq_abort('allocate_cs: error(1)')
    endif
    allocate(set%lab_cover(set%mx_gcover),STAT=stat)
    if(stat/=0) then
       call cq_abort('allocate_cs: error(2)')
    endif
    !TM temporary ???
    allocate(set%inv_lab_cover(set%mx_gcover),STAT=stat)
    if(stat/=0) then
       call cq_abort('allocate_cs: error(3)')
    endif
    !TM temporary ???
    allocate(set%ncover_rem(3*mx_node),STAT=stat)
    if(stat/=0) then
       call cq_abort('allocate_cs: error(4)')
    endif
    call reg_alloc_mem(area_index,3*set%mx_gcover+3*mx_node,type_int)
    ! Be safe with pointer
    nullify(set%iprim_group)
    if(members) then
       allocate(set%n_ing_cover(set%mx_gcover),STAT=stat)
       if(stat/=0) then
          call cq_abort('allocate_cs: error(5)')
       endif
       allocate(set%icover_ibeg(set%mx_gcover),STAT=stat)
       if(stat/=0) then
          call cq_abort('allocate_cs: error(6)')
       endif
       call reg_alloc_mem(area_index,2*set%mx_gcover,type_int)
       ! These will be allocated later !
       !allocate(set%xcover(set%mx_mcover),STAT=stat)
       !if(stat/=0) then
       !   call cq_abort('allocate_cs: error(7)')
       !endif
       !allocate(set%ycover(set%mx_mcover),STAT=stat)
       !if(stat/=0) then
       !   call cq_abort('allocate_cs: error(8)')
       !endif
       !allocate(set%zcover(set%mx_mcover),STAT=stat)
       !if(stat/=0) then
       !   call cq_abort('allocate_cs: error(9)')
       !endif
    else
       nullify(set%n_ing_cover,set%icover_ibeg, &
            set%xcover,set%ycover,set%zcover)
    endif
    call stop_timer(tmr_std_allocation)
    return
  end subroutine allocate_cs
!!***

!!****f* cover_module/deallocate_cs *
!!
!!  NAME
!!   deallocate_cs
!!  USAGE
!!
!!  PURPOSE
!!   Deallocates memory from CS derived type
!!  INPUTS
!!
!!
!!  USES
!!
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99
!!  MODIFICATION HISTORY
!!   21/06/2001 dave
!!    Added ROBODoc header and cq_abort
!!   2008/05/16 ast
!!    Added timer
!!   2018/07/04 08:47 dave
!!    Added members logical to be consistent with allocate
!!   2018/09/06 10:35 dave
!!    Tidying and bug fix: members was bracketing ALL variables and
!!    stat wasn't initialised
!!  SOURCE
!!
  subroutine deallocate_cs(set,members)

    ! Module usage
    use basic_types
    use GenComms, ONLY: cq_abort
    use global_module, ONLY: area_index
    use memory_module, ONLY: reg_dealloc_mem, type_int

    implicit none

    ! Passed variables
    type(cover_set) :: set
    logical :: members

    ! Local variables
    integer :: stat

    stat = 0
    call start_timer(tmr_std_allocation)
    call reg_dealloc_mem(area_index,2*set%mx_gcover,type_int)
    call reg_dealloc_mem(area_index,3*set%mx_gcover,type_int)
    call reg_dealloc_mem(area_index,size(set%ncover_rem),type_int)
    if(ASSOCIATED(set%iprim_group)) then
       call reg_dealloc_mem(area_index,size(set%iprim_group),type_int)
       deallocate(set%iprim_group,STAT=stat)
       if(stat/=0) call cq_abort('deallocate_cs: error(1)')
    endif
    deallocate(set%ncover_rem, set%inv_lab_cover, set%lab_cover,set%lab_cell, STAT=stat)
    if(stat/=0) call cq_abort('deallocate_cs: error(2)')
    if(members) then
       deallocate(set%zcover,set%ycover,set%xcover, set%icover_ibeg,set%n_ing_cover, STAT=stat)
       if(stat/=0) call cq_abort('deallocate_cs: error(3)')
    endif
    call stop_timer(tmr_std_allocation)
    return
  end subroutine deallocate_cs
!!***

end module cover_module

max_stars_repo_name: negf/CONQUEST-release

max_issues_repo_name: negf/CONQUEST-release

max_forks_repo_name: negf/CONQUEST-release


text:
      SUBROUTINE PSLAED2( ICTXT, K, N, N1, NB, D, DROW, DCOL, Q, LDQ,
     $                    RHO, Z, W, DLAMDA, Q2, LDQ2, QBUF, CTOT, PSM,
     $                    NPCOL, INDX, INDXC, INDXP, INDCOL, COLTYP, NN,
     $                    NN1, NN2, IB1, IB2 )
*
*  -- ScaLAPACK auxiliary routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     December 31, 1998
*
*     .. Scalar Arguments ..
      INTEGER            DCOL, DROW, IB1, IB2, ICTXT, K, LDQ, LDQ2, N,
     $                   N1, NB, NN, NN1, NN2, NPCOL
      REAL               RHO
*     ..
*     .. Array Arguments ..
      INTEGER            COLTYP( * ), CTOT( 0: NPCOL-1, 4 ),
     $                   INDCOL( N ), INDX( * ), INDXC( * ), INDXP( * ),
     $                   PSM( 0: NPCOL-1, 4 )
      REAL               D( * ), DLAMDA( * ), Q( LDQ, * ),
     $                   Q2( LDQ2, * ), QBUF( * ), W( * ), Z( * )
*     ..
*
*  Purpose
*  =======
*
*  PSLAED2 sorts the two sets of eigenvalues together into a single
*  sorted set.  Then it tries to deflate the size of the problem.
*  There are two ways in which deflation can occur:  when two or more
*  eigenvalues are close together or if there is a tiny entry in the
*  Z vector.  For each such occurrence the order of the related secular
*  equation problem is reduced by one.
*
*  Arguments
*  =========
*
*  ICTXT  (global input) INTEGER
*         The BLACS context handle, indicating the global context of
*         the operation on the matrix. The context itself is global.
*
*  K      (output) INTEGER
*         The number of non-deflated eigenvalues, and the order of the
*         related secular equation. 0 <= K <=N.
*
*  N      (input) INTEGER
*         The dimension of the symmetric tridiagonal matrix.  N >= 0.
*
*  N1     (input) INTEGER
*         The location of the last eigenvalue in the leading sub-matrix.
*         min(1,N) < N1 < N.
*
*  NB      (global input) INTEGER
*          The blocking factor used to distribute the columns of the
*          matrix. NB >= 1.
*
*  D      (input/output) REAL array, dimension (N)
*         On entry, D contains the eigenvalues of the two submatrices to
*         be combined.
*         On exit, D contains the trailing (N-K) updated eigenvalues
*         (those which were deflated) sorted into increasing order.
*
*  DROW   (global input) INTEGER
*          The process row over which the first row of the matrix D is
*          distributed. 0 <= DROW < NPROW.
*
*  DCOL   (global input) INTEGER
*          The process column over which the first column of the
*          matrix D is distributed. 0 <= DCOL < NPCOL.
*
*  Q      (input/output) REAL array, dimension (LDQ, N)
*         On entry, Q contains the eigenvectors of two submatrices in
*         the two square blocks with corners at (1,1), (N1,N1)
*         and (N1+1, N1+1), (N,N).
*         On exit, Q contains the trailing (N-K) updated eigenvectors
*         (those which were deflated) in its last N-K columns.
*
*  LDQ    (input) INTEGER
*         The leading dimension of the array Q.  LDQ >= max(1,NQ).
*
*  RHO    (global input/output) REAL
*         On entry, the off-diagonal element associated with the rank-1
*         cut which originally split the two submatrices which are now
*         being recombined.
*         On exit, RHO has been modified to the value required by
*         PSLAED3.
*
*  Z      (global input) REAL array, dimension (N)
*         On entry, Z contains the updating vector (the last
*         row of the first sub-eigenvector matrix and the first row of
*         the second sub-eigenvector matrix).
*         On exit, the contents of Z have been destroyed by the updating
*         process.
*
*  DLAMDA (global output) REAL array, dimension (N)
*         A copy of the first K eigenvalues which will be used by
*         SLAED3 to form the secular equation.
*
*  W      (global output) REAL array, dimension (N)
*         The first k values of the final deflation-altered z-vector
*         which will be passed to SLAED3.
*
*  Q2     (output) REAL array, dimension (LDQ2, NQ)
*         A copy of the first K eigenvectors which will be used by
*
*  LDQ2    (input) INTEGER
*         The leading dimension of the array Q2.
*
*  QBUF   (workspace) REAL array, dimension 3*N
*
*  CTOT   (workspace) INTEGER array, dimension( NPCOL, 4)
*
*  PSM    (workspace) INTEGER array, dimension( NPCOL, 4)
*
*  NPCOL   (global input) INTEGER
*          The total number of columns over which the distributed
*           submatrix is distributed.
*
*  INDX   (workspace) INTEGER array, dimension (N)
*         The permutation used to sort the contents of DLAMDA into
*         ascending order.
*
*  INDXC  (output) INTEGER array, dimension (N)
*         The permutation used to arrange the columns of the deflated
*         Q matrix into three groups:  the first group contains non-zero
*         elements only at and above N1, the second contains
*         non-zero elements only below N1, and the third is dense.
*
*  INDXP  (workspace) INTEGER array, dimension (N)
*         The permutation used to place deflated values of D at the end
*         of the array.  INDXP(1:K) points to the nondeflated D-values
*         and INDXP(K+1:N) points to the deflated eigenvalues.
*
*  INDCOL (workspace) INTEGER array, dimension (N)
*
*  COLTYP (workspace/output) INTEGER array, dimension (N)
*         During execution, a label which will indicate which of the
*         following types a column in the Q2 matrix is:
*         1 : non-zero in the upper half only;
*         2 : dense;
*         3 : non-zero in the lower half only;
*         4 : deflated.
*
*  NN     (global output) INTEGER, the order of matrix U, (PSLAED1).
*  NN1    (global output) INTEGER, the order of matrix Q1, (PSLAED1).
*  NN2    (global output) INTEGER, the order of matrix Q2, (PSLAED1).
*  IB1    (global output) INTEGER, pointeur on Q1, (PSLAED1).
*  IB2    (global output) INTEGER, pointeur on Q2, (PSLAED1).
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               MONE, ZERO, ONE, TWO, EIGHT
      PARAMETER          ( MONE = -1.0E0, ZERO = 0.0E0, ONE = 1.0E0,
     $                   TWO = 2.0E0, EIGHT = 8.0E0 )
*     ..
*     .. Local Scalars ..
      INTEGER            COL, CT, I, IAM, IE1, IE2, IMAX, INFO, J, JJQ2,
     $                   JJS, JMAX, JS, K2, MYCOL, MYROW, N1P1, N2, NJ,
     $                   NJCOL, NJJ, NP, NPROCS, NPROW, PJ, PJCOL, PJJ
      REAL               C, EPS, S, T, TAU, TOL
*     ..
*     .. External Functions ..
      INTEGER            INDXG2L, INDXL2G, ISAMAX, NUMROC
      REAL               PSLAMCH, SLAPY2
      EXTERNAL           INDXG2L, INDXL2G, ISAMAX, NUMROC, PSLAMCH,
     $                   SLAPY2
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, BLACS_PINFO, INFOG1L, SCOPY,
     $                   SGERV2D, SGESD2D, SLAPST, SROT, SSCAL
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, MIN, MOD, SQRT
*     ..
*     .. External Functions ..
*     ..
*     .. Local Arrays ..
      INTEGER            PTT( 4 )
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( N.EQ.0 )
     $   RETURN
*
      CALL BLACS_PINFO( IAM, NPROCS )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
      NP = NUMROC( N, NB, MYROW, DROW, NPROW )
*
      N2 = N - N1
      N1P1 = N1 + 1
*
      IF( RHO.LT.ZERO ) THEN
         CALL SSCAL( N2, MONE, Z( N1P1 ), 1 )
      END IF
*
*     Normalize z so that norm(z) = 1.  Since z is the concatenation of
*     two normalized vectors, norm2(z) = sqrt(2).
*
      T = ONE / SQRT( TWO )
      CALL SSCAL( N, T, Z, 1 )
*
*     RHO = ABS( norm(z)**2 * RHO )
*
      RHO = ABS( TWO*RHO )
*
*     Calculate the allowable deflation tolerance
*
      IMAX = ISAMAX( N, Z, 1 )
      JMAX = ISAMAX( N, D, 1 )
      EPS = PSLAMCH( ICTXT, 'Epsilon' )
      TOL = EIGHT*EPS*MAX( ABS( D( JMAX ) ), ABS( Z( IMAX ) ) )
*
*     If the rank-1 modifier is small enough, no more needs to be done
*     except to reorganize Q so that its columns correspond with the
*     elements in D.
*
      IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
         K = 0
         GO TO 220
      END IF
*
*     If there are multiple eigenvalues then the problem deflates.  Here
*     the number of equal eigenvalues are found.  As each equal
*     eigenvalue is found, an elementary reflector is computed to rotate
*     the corresponding eigensubspace so that the corresponding
*     components of Z are zero in this new basis.
*
*
      CALL SLAPST( 'I', N, D, INDX, INFO )
*
      DO 10 I = 1, N1
         COLTYP( I ) = 1
   10 CONTINUE
      DO 20 I = N1P1, N
         COLTYP( I ) = 3
   20 CONTINUE
      COL = DCOL
      DO 40 I = 1, N, NB
         DO 30 J = 0, NB - 1
            IF( I+J.LE.N )
     $         INDCOL( I+J ) = COL
   30    CONTINUE
         COL = MOD( COL+1, NPCOL )
   40 CONTINUE
*
      K = 0
      K2 = N + 1
      DO 50 J = 1, N
         NJ = INDX( J )
         IF( RHO*ABS( Z( NJ ) ).LE.TOL ) THEN
*
*           Deflate due to small z component.
*
            K2 = K2 - 1
            COLTYP( NJ ) = 4
            INDXP( K2 ) = NJ
            IF( J.EQ.N )
     $         GO TO 80
         ELSE
            PJ = NJ
            GO TO 60
         END IF
   50 CONTINUE
   60 CONTINUE
      J = J + 1
      NJ = INDX( J )
      IF( J.GT.N )
     $   GO TO 80
      IF( RHO*ABS( Z( NJ ) ).LE.TOL ) THEN
*
*        Deflate due to small z component.
*
         K2 = K2 - 1
         COLTYP( NJ ) = 4
         INDXP( K2 ) = NJ
      ELSE
*
*        Check if eigenvalues are close enough to allow deflation.
*
         S = Z( PJ )
         C = Z( NJ )
*
*        Find sqrt(a**2+b**2) without overflow or
*        destructive underflow.
*
         TAU = SLAPY2( C, S )
         T = D( NJ ) - D( PJ )
         C = C / TAU
         S = -S / TAU
         IF( ABS( T*C*S ).LE.TOL ) THEN
*
*           Deflation is possible.
*
            Z( NJ ) = TAU
            Z( PJ ) = ZERO
            IF( COLTYP( NJ ).NE.COLTYP( PJ ) )
     $         COLTYP( NJ ) = 2
            COLTYP( PJ ) = 4
            CALL INFOG1L( NJ, NB, NPCOL, MYCOL, DCOL, NJJ, NJCOL )
            CALL INFOG1L( PJ, NB, NPCOL, MYCOL, DCOL, PJJ, PJCOL )
            IF( INDCOL( PJ ).EQ.INDCOL( NJ ) .AND. MYCOL.EQ.NJCOL ) THEN
               CALL SROT( NP, Q( 1, PJJ ), 1, Q( 1, NJJ ), 1, C, S )
            ELSE IF( MYCOL.EQ.PJCOL ) THEN
               CALL SGESD2D( ICTXT, NP, 1, Q( 1, PJJ ), NP, MYROW,
     $                       NJCOL )
               CALL SGERV2D( ICTXT, NP, 1, QBUF, NP, MYROW, NJCOL )
               CALL SROT( NP, Q( 1, PJJ ), 1, QBUF, 1, C, S )
            ELSE IF( MYCOL.EQ.NJCOL ) THEN
               CALL SGESD2D( ICTXT, NP, 1, Q( 1, NJJ ), NP, MYROW,
     $                       PJCOL )
               CALL SGERV2D( ICTXT, NP, 1, QBUF, NP, MYROW, PJCOL )
               CALL SROT( NP, QBUF, 1, Q( 1, NJJ ), 1, C, S )
            END IF
            T = D( PJ )*C**2 + D( NJ )*S**2
            D( NJ ) = D( PJ )*S**2 + D( NJ )*C**2
            D( PJ ) = T
            K2 = K2 - 1
            I = 1
   70       CONTINUE
            IF( K2+I.LE.N ) THEN
               IF( D( PJ ).LT.D( INDXP( K2+I ) ) ) THEN
                  INDXP( K2+I-1 ) = INDXP( K2+I )
                  INDXP( K2+I ) = PJ
                  I = I + 1
                  GO TO 70
               ELSE
                  INDXP( K2+I-1 ) = PJ
               END IF
            ELSE
               INDXP( K2+I-1 ) = PJ
            END IF
            PJ = NJ
         ELSE
            K = K + 1
            DLAMDA( K ) = D( PJ )
            W( K ) = Z( PJ )
            INDXP( K ) = PJ
            PJ = NJ
         END IF
      END IF
      GO TO 60
   80 CONTINUE
*
*     Record the last eigenvalue.
*
      K = K + 1
      DLAMDA( K ) = D( PJ )
      W( K ) = Z( PJ )
      INDXP( K ) = PJ
*
*     Count up the total number of the various types of columns, then
*     form a permutation which positions the four column types into
*     four uniform groups (although one or more of these groups may be
*     empty).
*
      DO 100 J = 1, 4
         DO 90 I = 0, NPCOL - 1
            CTOT( I, J ) = 0
   90    CONTINUE
         PTT( J ) = 0
  100 CONTINUE
      DO 110 J = 1, N
         CT = COLTYP( J )
         COL = INDCOL( J )
         CTOT( COL, CT ) = CTOT( COL, CT ) + 1
  110 CONTINUE
*
*     PSM(*) = Position in SubMatrix (of types 1 through 4)
*
      DO 120 COL = 0, NPCOL - 1
         PSM( COL, 1 ) = 1
         PSM( COL, 2 ) = 1 + CTOT( COL, 1 )
         PSM( COL, 3 ) = PSM( COL, 2 ) + CTOT( COL, 2 )
         PSM( COL, 4 ) = PSM( COL, 3 ) + CTOT( COL, 3 )
  120 CONTINUE
      PTT( 1 ) = 1
      DO 140 I = 2, 4
         CT = 0
         DO 130 J = 0, NPCOL - 1
            CT = CT + CTOT( J, I-1 )
  130    CONTINUE
         PTT( I ) = PTT( I-1 ) + CT
  140 CONTINUE
*
*     Fill out the INDXC array so that the permutation which it induces
*     will place all type-1 columns first, all type-2 columns next,
*     then all type-3's, and finally all type-4's.
*
      DO 150 J = 1, N
         JS = INDXP( J )
         COL = INDCOL( JS )
         CT = COLTYP( JS )
         I = INDXL2G( PSM( COL, CT ), NB, COL, DCOL, NPCOL )
         INDX( J ) = I
         INDXC( PTT( CT ) ) = I
         PSM( COL, CT ) = PSM( COL, CT ) + 1
         PTT( CT ) = PTT( CT ) + 1
  150 CONTINUE
      DO 160 J = 1, N
         JS = INDXP( J )
         JJS = INDXG2L( JS, NB, J, J, NPCOL )
         COL = INDCOL( JS )
         IF( COL.EQ.MYCOL ) THEN
            I = INDX( J )
            JJQ2 = INDXG2L( I, NB, J, J, NPCOL )
            CALL SCOPY( NP, Q( 1, JJS ), 1, Q2( 1, JJQ2 ), 1 )
         END IF
  160 CONTINUE
*
*
*     The deflated eigenvalues and their corresponding vectors go back
*     into the last N - K slots of D and Q respectively.
*
      CALL SCOPY( N, D, 1, Z, 1 )
      DO 170 J = K + 1, N
         JS = INDXP( J )
         I = INDX( J )
         D( I ) = Z( JS )
  170 CONTINUE
*
      PTT( 1 ) = 1
      DO 190 I = 2, 4
         CT = 0
         DO 180 J = 0, NPCOL - 1
            CT = CT + CTOT( J, I-1 )
  180    CONTINUE
         PTT( I ) = PTT( I-1 ) + CT
  190 CONTINUE
*
*
      IB1 = INDXC( 1 )
      IE1 = IB1
      IB2 = INDXC( PTT( 2 ) )
      IE2 = IB2
      DO 200 I = 2, PTT( 3 ) - 1
         IB1 = MIN( IB1, INDXC( I ) )
         IE1 = MAX( IE1, INDXC( I ) )
  200 CONTINUE
      DO 210 I = PTT( 2 ), PTT( 4 ) - 1
         IB2 = MIN( IB2, INDXC( I ) )
         IE2 = MAX( IE2, INDXC( I ) )
  210 CONTINUE
      NN1 = IE1 - IB1 + 1
      NN2 = IE2 - IB2 + 1
      NN = MAX( IE1, IE2 ) - MIN( IB1, IB2 ) + 1
  220 CONTINUE
      RETURN
*
*     End of PSLAED2
*
      END

ext: f

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack