text:
module response_matrix

  use netcdf
  use mpi

  implicit none

  public :: init_response_matrix, finish_response_matrix
  public :: read_response_matrix
  public :: response_matrix_initialized

  private

  logical :: response_matrix_initialized = .false.
  integer, parameter :: mat_unit = 70

#if defined MPI && defined ISO_C_BINDING
  integer :: window = MPI_WIN_NULL
#endif

contains

  subroutine init_response_matrix

    use linear_solve, only: lu_decomposition
    use fields_arrays, only: response_matrix
    use stella_layouts, only: vmu_lo
    use stella_layouts, only: iv_idx, is_idx
    use kt_grids, only: naky
    use extended_zgrid, only: iz_low, iz_up
    use extended_zgrid, only: neigen, ikxmod
    use extended_zgrid, only: nsegments
    use extended_zgrid, only: nzed_segment
    use extended_zgrid, only: periodic
    use job_manage, only: time_message
    use mp, only: proc0, job, mp_abort
    use run_parameters, only: mat_gen, lu_option_switch
    use run_parameters, only: lu_option_none, lu_option_local, lu_option_global
    use system_fortran, only: systemf
#if defined MPI && defined ISO_C_BINDING
    use, intrinsic :: iso_c_binding, only: c_ptr, c_f_pointer
    use mp, only: curr_focus, sgproc0, mp_comm, sharedsubprocs, scope, barrier
    use mpi
#endif

    implicit none

    integer :: iky, ie, iseg, iz
    integer :: ikx
    integer :: nz_ext, nresponse
    integer :: idx
    integer :: izl_offset, izup
#if defined MPI && defined ISO_C_BINDING
    integer :: prior_focus, ierr, nbytes_real, rec_len
    integer :: disp_unit = 1
    integer (kind=MPI_ADDRESS_KIND) :: win_size, real_size
    integer*8 :: cur_pos
    type(c_ptr) :: bptr, cptr
#endif
    real :: dum
    complex, dimension (:), allocatable :: phiext
    complex, dimension (:,:), allocatable :: gext
    logical :: debug = .false.
    character(100) :: message_dgdphi, message_QN, message_lu
    real, dimension (2) :: time_response_matrix_dgdphi
    real, dimension (2) :: time_response_matrix_QN
    real, dimension (2) :: time_response_matrix_lu

    ! Related to the saving of the the matrices in netcdf format
    character(len=15) :: fmt, job_str
    character(len=100) :: file_name
    integer :: istatus
    istatus = 0

    if (proc0.and.debug) then
      write (*,*) " "
      write (*,'(A)') "    ############################################################"
      write (*,'(A)') "                         RESPONSE MATRIX"
      write (*,'(A)') "    ############################################################"
    end if
    message_dgdphi = '     calculate dgdphi: '
    message_QN     = '     calculate QN:     '
    message_lu     = '     calculate LU:     '
    time_response_matrix_dgdphi = 0
    time_response_matrix_QN     = 0
    time_response_matrix_lu     = 0

!   All matrices handled by processor i_proc and job are stored
!   on a single file named: response_mat_job.iproc
    fmt = '(I5.5)'
    if (proc0.and.mat_gen) then

       call systemf('mkdir -p mat')

       write (job_str, '(I1.1)') job
       file_name = './mat/response_mat_'//trim(job_str)

       open(unit=mat_unit, status='replace', file=file_name, &
            position='rewind', action='write', form='unformatted')
       write(unit=mat_unit) naky
    end if

    if (response_matrix_initialized) return
    response_matrix_initialized = .true.

    if (.not.allocated(response_matrix)) allocate (response_matrix(naky))

#if defined ISO_C_BINDING && defined MPI
    inquire(iolength=rec_len) dum !this will return 2 or 8 for double
    if(rec_len.eq.1.or.rec_len.eq.4) then
      nbytes_real = 4
      real_size = 4_MPI_ADDRESS_KIND
    else if (rec_len.eq.2.or.rec_len.eq.8) then
      nbytes_real = 8
      real_size = 8_MPI_ADDRESS_KIND
    else
      call mp_abort('failure retrieving the size of a real')
    endif

!   Create a single shared memory window for all the response matrices and
!   permutation arrays.
!   Creating a window for each matrix/array would lead to performance
!   degradation on some clusters
    if(window.eq.MPI_WIN_NULL) then
      prior_focus = curr_focus
      call scope(sharedsubprocs)
      win_size = 0
      if(sgproc0) then
        do iky = 1, naky
          do ie = 1, neigen(iky)
            if (periodic(iky)) then
               nresponse = nsegments(ie,iky)*nzed_segment
            else
               nresponse = nsegments(ie,iky)*nzed_segment+1
            end if
            win_size =   win_size &
                       + int(nresponse,MPI_ADDRESS_KIND)*4_MPI_ADDRESS_KIND &
                       + int(nresponse**2,MPI_ADDRESS_KIND)*2*real_size
          enddo
        enddo
      endif
      call mpi_win_allocate_shared(win_size,disp_unit,MPI_INFO_NULL,mp_comm, &
                                   bptr,window,ierr)

      if(.not.sgproc0) then
        !make sure all the procs have the right memory address
        call mpi_win_shared_query(window, 0, win_size, disp_unit, bptr, ierr)
      endif
      call mpi_win_fence(0,window,ierr)

      !the following is a hack that allows us to perform pointer arithmetic in Fortran
      cur_pos = transfer(bptr,cur_pos)

      call scope(prior_focus)
    endif
#endif

    ! for a given ky and set of connected kx values
    ! give a unit impulse to phi at each zed location
    ! in the extended domain and solve for h(zed_extended,(vpa,mu,s))

    do iky = 1, naky

       if (proc0.and.mat_gen) THEN
          write(unit=mat_unit) iky, neigen(iky)
       end if

       ! the response matrix for each ky has neigen(ky)
       ! independent sets of connected kx values
       if (.not.associated(response_matrix(iky)%eigen)) &
            allocate (response_matrix(iky)%eigen(neigen(iky)))

       if(proc0.and.debug) then
         call time_message(.false., time_response_matrix_dgdphi, message_dgdphi)
       end if

       ! loop over the sets of connected kx values
       do ie = 1, neigen(iky)

          ! number of zeds x number of segments
          nz_ext = nsegments(ie,iky)*nzed_segment+1

          ! treat zonal mode specially to avoid double counting
          ! as it is periodic
          if (periodic(iky)) then
             nresponse = nz_ext-1
          else
             nresponse = nz_ext
          end if

          if (proc0.and.mat_gen) then
             write(unit=mat_unit) ie, nresponse
          end if

#if !defined ISO_C_BINDING || !defined MPI
          ! for each ky and set of connected kx values,
          ! must have a response matrix that is N x N
          ! with N = number of zeds per 2pi segment x number of 2pi segments
          if (.not.associated(response_matrix(iky)%eigen(ie)%zloc)) &
               allocate (response_matrix(iky)%eigen(ie)%zloc(nresponse,nresponse))

          ! response_matrix%idx is needed to keep track of permutations
          ! to the response matrix made during LU decomposition
          ! it will be input to LU back substitution during linear solve
          if (.not.associated(response_matrix(iky)%eigen(ie)%idx)) &
               allocate (response_matrix(iky)%eigen(ie)%idx(nresponse))
#else
          !exploit MPIs shared memory framework to reduce memory consumption of
          !response matrices

          if (.not.associated(response_matrix(iky)%eigen(ie)%zloc)) then
            cptr = transfer(cur_pos,cptr)
            call c_f_pointer(cptr,response_matrix(iky)%eigen(ie)%zloc,(/nresponse,nresponse/))
            cur_pos = cur_pos + nresponse**2*2*nbytes_real
          endif

          if (.not.associated(response_matrix(iky)%eigen(ie)%idx)) then
            cptr = transfer(cur_pos,cptr)
            call c_f_pointer(cptr,response_matrix(iky)%eigen(ie)%idx,(/nresponse/))
            cur_pos = cur_pos + nresponse*4
          endif
#endif

          allocate (gext(nz_ext,vmu_lo%llim_proc:vmu_lo%ulim_alloc))
          allocate (phiext(nz_ext))
          ! idx is the index in the extended zed domain
          ! that we are giving a unit impulse
          idx = 0

          ! loop over segments, starting with 1
          ! first segment is special because it has
          ! one more unique zed value than all others
          ! since domain is [z0-pi:z0+pi], including both endpoints
          ! i.e., one endpoint is shared with the previous segment
          iseg = 1
          ! ikxmod gives the kx corresponding to iseg,ie,iky
          ikx = ikxmod(iseg,ie,iky)
          izl_offset = 0
          ! avoid double-counting of periodic points for zonal mode (and other periodic modes)
          if (periodic(iky)) then
             izup = iz_up(iseg)-1
          else
             izup = iz_up(iseg)
          end if
          ! no need to obtain response to impulses at negative kx values
          do iz = iz_low(iseg), izup
             idx = idx + 1
             call get_dgdphi_matrix_column (iky, ikx, iz, ie, idx, nz_ext, nresponse, phiext, gext)
          end do
          ! once we have used one segments, remaining segments
          ! have one fewer unique zed point
          izl_offset = 1
          if (nsegments(ie,iky) > 1) then
             do iseg = 2, nsegments(ie,iky)
                ikx = ikxmod(iseg,ie,iky)
                do iz = iz_low(iseg)+izl_offset, iz_up(iseg)
                   idx = idx + 1
                   call get_dgdphi_matrix_column (iky, ikx, iz, ie, idx, nz_ext, nresponse, phiext, gext)
                end do
                if (izl_offset == 0) izl_offset = 1
             end do
          end if
          deallocate (gext,phiext)
       enddo
       !DSO - This ends parallelization over velocity space.
       !      At this point every processor has int dv dgdphi for a given ky
       !      and so the quasineutrality solve and LU decomposition can be
       !      parallelized locally if need be.
       !      This is preferable to parallelization over ky as the LU
       !      decomposition (and perhaps QN) will be dominated by the
       !      ky with the most connections

       if(proc0.and.debug) then
         call time_message(.true. , time_response_matrix_dgdphi, message_dgdphi)
         call time_message(.false., time_response_matrix_QN, message_QN)
       end if

#ifdef ISO_C_BINDING
       call mpi_win_fence(0,window,ierr)
#endif

       ! solve quasineutrality
       ! for local stella, this is a diagonal process, but global stella
       ! may require something more sophisticated

       ! loop over the sets of connected kx values
       do ie = 1, neigen(iky)
#if defined ISO_C_BINDING && defined MPI
         if(sgproc0) then
#endif
           ! number of zeds x number of segments
           nz_ext = nsegments(ie,iky)*nzed_segment+1

           ! treat zonal mode specially to avoid double counting
           ! as it is periodic
           if (periodic(iky)) then
             nresponse = nz_ext-1
           else
             nresponse = nz_ext
           end if

           allocate (phiext(nz_ext))

           do idx = 1, nresponse
             phiext(nz_ext) = 0.0
             phiext(:nresponse) = response_matrix(iky)%eigen(ie)%zloc(:,idx)
             call get_fields_for_response_matrix (phiext, iky, ie)

             ! next need to create column in response matrix from phiext
             ! negative sign because matrix to be inverted in streaming equation
             ! is identity matrix - response matrix
             ! add in contribution from identity matrix
             phiext(idx) = phiext(idx)-1.0
             response_matrix(iky)%eigen(ie)%zloc(:,idx) = -phiext(:nresponse)

           end do
           deallocate (phiext)
#if defined ISO_C_BINDING && defined MPI
         endif
#endif
       enddo

#ifdef ISO_C_BINDING
       call mpi_win_fence(0,window,ierr)
#endif

       if(proc0.and.debug) then
         call time_message(.true. , time_response_matrix_QN, message_QN)
         call time_message(.false., time_response_matrix_lu, message_lu)
       end if

       !now we have the full response matrix. Finally, perform its LU decomposition
#ifdef MPI
       select case (lu_option_switch)
       case (lu_option_global)
        call parallel_LU_decomposition_global(iky)
       case (lu_option_local)
#ifdef ISO_C_BINDING
         call parallel_LU_decomposition_local(iky)
#else
         call mp_abort('Stella must be built with HAS_ISO_BINDING in order to use local parallel LU decomposition.')
#endif
       case default
#endif
         do ie = 1, neigen(iky)
#ifdef ISO_C_BINDING
           if(sgproc0) then
#endif
             ! now that we have the reponse matrix for this ky and set of connected kx values
             !get the LU decomposition so we are ready to solve the linear system
             call lu_decomposition (response_matrix(iky)%eigen(ie)%zloc, &
                                    response_matrix(iky)%eigen(ie)%idx,dum)

#ifdef ISO_C_BINDING
           endif
#endif
        enddo
#ifdef MPI
       end select
#endif

       if(proc0.and.debug) then
         call time_message(.true., time_response_matrix_lu, message_lu)
       end if

       time_response_matrix_dgdphi = 0
       time_response_matrix_QN     = 0
       time_response_matrix_lu     = 0

       do ie = 1, neigen(iky)
          if (proc0.and.mat_gen) then
             write(unit=mat_unit) response_matrix(iky)%eigen(ie)%idx
             write(unit=mat_unit) response_matrix(iky)%eigen(ie)%zloc
          end if
       end do

       !if(proc0)  write (*,*) 'job', iky, iproc, response_matrix(iky)%eigen(1)%zloc(5,:)
    end do

#ifdef ISO_C_BINDING
    call mpi_win_fence(0,window,ierr)
#endif

    if (proc0.and.mat_gen) then
       close(unit=mat_unit)
    end if

    if (proc0.and.debug) then
      write (*,'(A)') "    ############################################################"
      write (*,'(A)') " "
    end if

  end subroutine init_response_matrix

  subroutine read_response_matrix

    use fields_arrays, only: response_matrix
    use common_types, only: response_matrix_type
    use kt_grids, only: naky
    use extended_zgrid, only: neigen
    use extended_zgrid, only: nsegments
    use extended_zgrid, only: nzed_segment
    use extended_zgrid, only: periodic
    use mp, only: proc0, job, broadcast, mp_abort

    implicit none

    integer :: iky, ie, nz_ext
    integer :: iky_dump, neigen_dump, naky_dump, nresponse_dump
    integer :: nresponse
    character(len=15) :: fmt, job_str
    character(len=100) :: file_name
    integer :: istatus, ie_dump, istat
    logical, parameter :: debug=.false.
    istatus = 0

!   All matrices handled for the job i_job are read
!   from a single file named: responst_mat.ijob by that
!   jobs root process

    if(proc0) then
      fmt = '(I5.5)'
      write (job_str, '(I1.1)') job
      file_name = './mat/response_mat.'//trim(job_str)

      open(unit=mat_unit, status='old', file=file_name, &
         action='read', form='unformatted', iostat=istat)
      if (istat /= 0) then
         print *, 'Error opening response_matrix by root processor for job ', job_str
      end if

      read(unit=mat_unit) naky_dump
      if(naky.ne.naky_dump) call mp_abort('mismatch in naky and naky_dump')
    endif

    if (.not.allocated(response_matrix)) allocate (response_matrix(naky))

    do iky = 1, naky
       if(proc0) then
         read(unit=mat_unit) iky_dump, neigen_dump
         if(iky_dump.ne.iky.or.neigen_dump.ne.neigen(iky)) &
           call mp_abort('mismatch in iky_dump/neigen_dump')
       endif

       if (.not.associated(response_matrix(iky)%eigen)) &
            allocate (response_matrix(iky)%eigen(neigen(iky)))

       ! loop over the sets of connected kx values
       do ie = 1, neigen(iky)
          ! number of zeds x number of segments
          nz_ext = nsegments(ie,iky)*nzed_segment+1

          ! treat zonal mode specially to avoid double counting
          ! as it is periodic
          if (periodic(iky)) then
             nresponse = nz_ext-1
          else
             nresponse = nz_ext
          end if

          if(proc0) then
            read(unit=mat_unit) ie_dump, nresponse_dump
            if(ie_dump.ne.ie.or.nresponse.ne.nresponse_dump) &
              call mp_abort('mismatch in ie/nresponse_dump')
          endif

          ! for each ky and set of connected kx values,
          ! must have a response matrix that is N x N
          ! with N = number of zeds per 2pi segment x number of 2pi segments
          if (.not.associated(response_matrix(iky)%eigen(ie)%zloc)) &
               allocate (response_matrix(iky)%eigen(ie)%zloc(nresponse,nresponse))

          ! response_matrix%idx is needed to keep track of permutations
          ! to the response matrix made during LU decomposition
          ! it will be input to LU back substitution during linear solve
          if (.not.associated(response_matrix(iky)%eigen(ie)%idx)) &
               allocate (response_matrix(iky)%eigen(ie)%idx(nresponse))
          if(proc0) then
            read(unit=mat_unit) response_matrix(iky)%eigen(ie)%idx
            read(unit=mat_unit) response_matrix(iky)%eigen(ie)%zloc
          endif

          call broadcast(response_matrix(iky)%eigen(ie)%idx)
          call broadcast(response_matrix(iky)%eigen(ie)%zloc)

       enddo
    enddo

    if (proc0) close (mat_unit)

    if (debug) then
       print *, 'File', file_name, ' successfully read by root proc for job: ', job_str
    end if
  end subroutine read_response_matrix

  subroutine get_dgdphi_matrix_column (iky, ikx, iz, ie, idx, nz_ext, nresponse, phiext, gext)

    use stella_layouts, only: vmu_lo
    use stella_layouts, only: iv_idx, imu_idx, is_idx
    use stella_time, only: code_dt
    use zgrid, only: delzed, nzgrid
    use extended_zgrid, only: periodic
    use species, only: spec
    use stella_geometry, only: gradpar, dbdzed
    use vpamu_grids, only: vpa, mu
    use vpamu_grids, only: maxwell_vpa, maxwell_mu, maxwell_fac
    use fields_arrays, only: response_matrix
    use gyro_averages, only: aj0x
    use run_parameters, only: driftkinetic_implicit
    use run_parameters, only: maxwellian_inside_zed_derivative
    use parallel_streaming, only: stream_tridiagonal_solve
    use parallel_streaming, only: stream_sign
    use run_parameters, only: zed_upwind, time_upwind
#if defined ISO_C_BINDING && defined MPI
    use mp, only: sgproc0
#endif

    implicit none

    integer, intent (in) :: iky, ikx, iz, ie, idx, nz_ext, nresponse
    complex, dimension (:), intent (in out) :: phiext
    complex, dimension (:,vmu_lo%llim_proc:), intent (in out) :: gext

    integer :: ivmu, iv, imu, is, ia
    integer :: izp, izm
    real :: mu_dbdzed_p, mu_dbdzed_m
    real :: fac, fac0, fac1, gyro_fac

    ia = 1

    if (.not.maxwellian_inside_zed_derivative) then
       ! get -vpa*b.gradz*Ze/T*F0*d<phi>/dz corresponding to unit impulse in phi
       do ivmu = vmu_lo%llim_proc, vmu_lo%ulim_proc
          ! initialize g to zero everywhere along extended zed domain
          gext(:,ivmu) = 0.0
          iv = iv_idx(vmu_lo,ivmu)
          imu = imu_idx(vmu_lo,ivmu)
          is = is_idx(vmu_lo,ivmu)

          ! give unit impulse to phi at this zed location
          ! and compute -vpa*b.gradz*Ze/T*d<phi>/dz*F0 (RHS of streaming part of GKE)

          ! NB:  assuming equal spacing in zed below
          ! here, fac = -dt*(1+alph_t)/2*vpa*Ze/T*F0*J0/dz
          ! b.gradz left out because needs to be centred in zed
          if (driftkinetic_implicit) then
             gyro_fac = 1.0
          else
             gyro_fac = aj0x(iky,ikx,iz,ivmu)
          end if

          ! 0.125 to account for two linear interpolations
          fac = -0.125*(1.+time_upwind)*code_dt*vpa(iv)*spec(is)%stm_psi0 &
               *gyro_fac*spec(is)%zt/delzed(0)*maxwell_vpa(iv,is)*maxwell_fac(is)

          ! In the following, gradpar and maxwell_mu are interpolated separately
          ! to ensure consistency to what is done in parallel_streaming.f90

          ! stream_sign < 0 corresponds to positive advection speed
          if (stream_sign(iv)<0) then
             if (iz > -nzgrid) then
                ! fac0 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at this zed index
                fac0 = fac*((1.+zed_upwind)*gradpar(iz) &
                            + (1.-zed_upwind)*gradpar(iz-1)) &
                          *((1.+zed_upwind)*maxwell_mu(ia,iz,imu,is) &
                            + (1.-zed_upwind)*maxwell_mu(ia,iz-1,imu,is))
                ! fac1 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at the zed index to the right of
                ! this one
                if (iz < nzgrid) then
                   fac1 = fac*((1.+zed_upwind)*gradpar(iz+1) &
                               + (1.-zed_upwind)*gradpar(iz)) &
                             *((1.+zed_upwind)*maxwell_mu(ia,iz+1,imu,is) &
                               + (1.-zed_upwind)*maxwell_mu(ia,iz,imu,is))
                else
                   fac1 = fac*((1.+zed_upwind)*gradpar(-nzgrid+1) &
                               + (1.-zed_upwind)*gradpar(nzgrid)) &
                             *((1.+zed_upwind)*maxwell_mu(ia,-nzgrid+1,imu,is) &
                               + (1.-zed_upwind)*maxwell_mu(ia,nzgrid,imu,is))
                end if
             else
                ! fac0 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at this zed index
                fac0 = fac*((1.+zed_upwind)*gradpar(iz) &
                            + (1.-zed_upwind)*gradpar(nzgrid-1)) &
                          *((1.+zed_upwind)*maxwell_mu(ia,iz,imu,is) &
                            + (1.-zed_upwind)*maxwell_mu(ia,nzgrid-1,imu,is))
                ! fac1 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at the zed index to the right of
                ! this one
                fac1 = fac*((1.+zed_upwind)*gradpar(iz+1)  &
                            + (1.-zed_upwind)*gradpar(iz)) &
                          *((1.+zed_upwind)*maxwell_mu(ia,iz+1,imu,is) &
                            + (1.-zed_upwind)*maxwell_mu(ia,iz,imu,is))
             end if
             gext(idx,ivmu) = fac0
             if (idx < nz_ext) gext(idx+1,ivmu) = -fac1
             ! zonal mode BC is periodic instead of zero, so must
             ! treat specially
             if (periodic(iky)) then
                if (idx == 1) then
                   gext(nz_ext,ivmu) = fac0
                else if (idx == nz_ext-1) then
                   gext(1,ivmu) = -fac1
                end if
             end if
          else
             if (iz < nzgrid) then
                ! fac0 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at this zed index
                fac0 = fac*((1.+zed_upwind)*gradpar(iz) &
                            + (1.-zed_upwind)*gradpar(iz+1)) &
                          *((1.+zed_upwind)*maxwell_mu(ia,iz,imu,is) &
                            + (1.-zed_upwind)*maxwell_mu(ia,iz+1,imu,is))
                ! fac1 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at the zed index to the left of
                ! this one
                if (iz > -nzgrid) then
                   fac1 = fac*((1.+zed_upwind)*gradpar(iz-1) &
                               + (1.-zed_upwind)*gradpar(iz)) &
                             *((1.+zed_upwind)*maxwell_mu(ia,iz-1,imu,is) &
                               + (1.-zed_upwind)*maxwell_mu(ia,iz,imu,is))
                else
                   fac1 = fac*((1.+zed_upwind)*gradpar(nzgrid-1) &
                               + (1.-zed_upwind)*gradpar(iz)) &
                             *((1.+zed_upwind)*maxwell_mu(ia,nzgrid-1,imu,is) &
                               + (1.-zed_upwind)*maxwell_mu(ia,iz,imu,is))
                end if
             else
                ! fac0 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at this zed index
                fac0 = fac*((1.+zed_upwind)*gradpar(iz) &
                            + (1.-zed_upwind)*gradpar(-nzgrid+1)) &
                          *((1.+zed_upwind)*maxwell_mu(ia,iz,imu,is) &
                            + (1.-zed_upwind)*maxwell_mu(ia,-nzgrid+1,imu,is))
                ! fac1 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at the zed index to the left of
                ! this one
                fac1 = fac*((1.+zed_upwind)*gradpar(iz-1) &
                            + (1.-zed_upwind)*gradpar(iz)) &
                          *((1.+zed_upwind)*maxwell_mu(ia,iz-1,imu,is) &
                            + (1.-zed_upwind)*maxwell_mu(ia,iz,imu,is))
             end if
             gext(idx,ivmu) = -fac0
             if (idx > 1) gext(idx-1,ivmu) = fac1
             ! zonal mode BC is periodic instead of zero, so must
             ! treat specially
             if (periodic(iky)) then
                if (idx == 1) then
                   gext(nz_ext,ivmu) = -fac0
                   gext(nz_ext-1,ivmu) = fac1
                else if (idx == 2) then
                   gext(nz_ext,ivmu) = fac1
                end if
             end if
          end if

          ! hack for now (duplicates much of the effort from sweep_zed_zonal)
          if (periodic(iky)) then
             call sweep_zed_zonal_response (iv, is, stream_sign(iv), gext(:,ivmu))
          else
             ! invert parallel streaming equation to get g^{n+1} on extended zed grid
             ! (I + (1+alph)/2*dt*vpa)*g_{inh}^{n+1} = RHS = gext
             call stream_tridiagonal_solve (iky, ie, iv, is, gext(:,ivmu))
          end if

       end do
    else
       ! get -vpa*b.gradz*Ze/T*F0*d<phi>/dz corresponding to unit impulse in phi
       do ivmu = vmu_lo%llim_proc, vmu_lo%ulim_proc
          ! initialize g to zero everywhere along extended zed domain
          gext(:,ivmu) = 0.0
          iv = iv_idx(vmu_lo,ivmu)
          imu = imu_idx(vmu_lo,ivmu)
          is = is_idx(vmu_lo,ivmu)

          ! give unit impulse to phi at this zed location
          ! and compute -vpa*b.gradz*Ze/T*d<phi>/dz*F0 (RHS of streaming part of GKE)

          ! NB:  assuming equal spacing in zed below
          ! here, fac = -dt*(1+alph_t)/2*vpa*Ze/T*F0*J0/dz
          ! b.gradz left out because needs to be centred in zed
          if (driftkinetic_implicit) then
             gyro_fac = 1.0
          else
             gyro_fac = aj0x(iky,ikx,iz,ivmu)
          end if

          fac = -0.25*(1.+time_upwind)*code_dt*vpa(iv)*spec(is)%stm_psi0 &
               *gyro_fac*spec(is)%zt*maxwell_vpa(iv,is)*maxwell_mu(ia,iz,imu,is)*maxwell_fac(is)

          mu_dbdzed_p = 1./delzed(0)+mu(imu)*dbdzed(ia,iz)*(1.+zed_upwind)
          mu_dbdzed_m = 1./delzed(0)+mu(imu)*dbdzed(ia,iz)*(1.-zed_upwind)

          ! stream_sign < 0 corresponds to positive advection speed
          if (stream_sign(iv)<0) then
             if (iz > -nzgrid) then
                ! fac0 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at this zed index
                fac0 = fac*((1.+zed_upwind)*gradpar(iz) &
                     + (1.-zed_upwind)*gradpar(iz-1))*mu_dbdzed_p
                ! fac1 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at the zed index to the right of
                ! this one
                if (iz < nzgrid) then
                   izp = iz+1
                else
                   izp = -nzgrid+1
                end if
                fac1 = fac*((1.+zed_upwind)*gradpar(izp) &
                     + (1.-zed_upwind)*gradpar(iz))*mu_dbdzed_m
             else
                ! fac0 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at this zed index
                fac0 = fac*((1.+zed_upwind)*gradpar(iz) &
                     + (1.-zed_upwind)*gradpar(nzgrid-1))*mu_dbdzed_p
                ! fac1 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at the zed index to the right of
                ! this one
                fac1 = fac*((1.+zed_upwind)*gradpar(iz+1) &
                     + (1.-zed_upwind)*gradpar(iz))*mu_dbdzed_m
             end if
             gext(idx,ivmu) = fac0
             if (idx < nz_ext) gext(idx+1,ivmu) = -fac1
             ! zonal mode BC is periodic instead of zero, so must
             ! treat specially
             if (periodic(iky)) then
                if (idx == 1) then
                   gext(nz_ext,ivmu) = fac0
                else if (idx == nz_ext-1) then
                   gext(1,ivmu) = -fac1
                end if
             end if
          else
             if (iz < nzgrid) then
                ! fac0 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at this zed index
                fac0 = fac*((1.+zed_upwind)*gradpar(iz) &
                     + (1.-zed_upwind)*gradpar(iz+1))*mu_dbdzed_p
                ! fac1 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at the zed index to the left of
                ! this one
                if (iz > -nzgrid) then
                   izm = iz-1
                else
                   izm = nzgrid-1
                end if
                fac1 = fac*((1.+zed_upwind)*gradpar(izm) &
                     + (1.-zed_upwind)*gradpar(iz))*mu_dbdzed_m
             else
                ! fac0 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at this zed index
                fac0 = fac*((1.+zed_upwind)*gradpar(iz) &
                     + (1.-zed_upwind)*gradpar(-nzgrid+1))*mu_dbdzed_p
                ! fac1 is the factor multiplying delphi on the RHS
                ! of the homogeneous GKE at the zed index to the left of
                ! this one
                fac1 = fac*((1.+zed_upwind)*gradpar(iz-1) &
                     + (1.-zed_upwind)*gradpar(iz))*mu_dbdzed_m
             end if
             gext(idx,ivmu) = -fac0
             if (idx > 1) gext(idx-1,ivmu) = fac1
             ! zonal mode BC is periodic instead of zero, so must
             ! treat specially
             if (periodic(iky)) then
                if (idx == 1) then
                   gext(nz_ext,ivmu) = -fac0
                   gext(nz_ext-1,ivmu) = fac1
                else if (idx == 2) then
                   gext(nz_ext,ivmu) = fac1
                end if
             end if
          end if

          ! hack for now (duplicates much of the effort from sweep_zed_zonal)
          if (periodic(iky)) then
             call sweep_zed_zonal_response (iv, is, stream_sign(iv), gext(:,ivmu))
          else
             ! invert parallel streaming equation to get g^{n+1} on extended zed grid
             ! (I + (1+alph)/2*dt*vpa)*g_{inh}^{n+1} = RHS = gext
             call stream_tridiagonal_solve (iky, ie, iv, is, gext(:,ivmu))
          end if

       end do
    end if

    ! we now have g on the extended zed domain at this ky and set of connected kx values
    ! corresponding to a unit impulse in phi at this location
    ! now integrate over velocities to get a square response matrix
    ! (this ends the parallelization over velocity space, so every core should have a
    !  copy of phiext)
    call integrate_over_velocity (gext, phiext, iky, ie)

#if !defined ISO_C_BINDING || !defined MPI
    response_matrix(iky)%eigen(ie)%zloc(:,idx) = phiext(:nresponse)
#else
    if(sgproc0) response_matrix(iky)%eigen(ie)%zloc(:,idx) = phiext(:nresponse)
#endif

  end subroutine get_dgdphi_matrix_column

! subroutine get_phi_matrix
! end subroutine get_phi_matrix

  subroutine sweep_zed_zonal_response (iv, is, sgn, g)

    use zgrid, only: nzgrid, delzed, nztot
    use run_parameters, only: zed_upwind, time_upwind
    use parallel_streaming, only: stream_c

    implicit none

    integer, intent (in) :: iv, is, sgn
    complex, dimension (:), intent (in out) :: g

    integer :: iz, iz1, iz2
    real :: fac1, fac2
    complex, dimension (:), allocatable :: gcf, gpi

    allocate (gpi(-nzgrid:nzgrid))
    allocate (gcf(-nzgrid:nzgrid))
    ! ky=0 is 2pi periodic (no extended zgrid)
    ! decompose into complementary function + particular integral
    ! zero BC for particular integral
    ! unit BC for complementary function (no source)
    if (sgn < 0) then
       iz1 = -nzgrid ; iz2 = nzgrid
    else
       iz1 = nzgrid ; iz2 = -nzgrid
    end if
    gpi(iz1) = 0. ; gcf(iz1) = 1.
    do iz = iz1-sgn, iz2, -sgn
       fac1 = 1.0+zed_upwind+sgn*(1.0+time_upwind)*stream_c(iz,iv,is)/delzed(0)
       fac2 = 1.0-zed_upwind-sgn*(1.0+time_upwind)*stream_c(iz,iv,is)/delzed(0)
       gpi(iz) = (-gpi(iz+sgn)*fac2 + 2.0*g(iz+nzgrid+1))/fac1
       gcf(iz) = -gcf(iz+sgn)*fac2/fac1
    end do
    ! g = g_PI + (g_PI(pi)/(1-g_CF(pi))) * g_CF
    g = gpi + (spread(gpi(iz2),1,nztot)/(1.-gcf(iz2)))*gcf
    deallocate (gpi, gcf)

  end subroutine sweep_zed_zonal_response

  subroutine integrate_over_velocity(g,phi,iky,ie)

    use stella_layouts, only: vmu_lo
    use species, only: nspec, spec
    use extended_zgrid, only: iz_low, iz_up
    use extended_zgrid, only: ikxmod
    use extended_zgrid, only: nsegments
    use kt_grids, only: zonal_mode, akx
    use vpamu_grids, only: integrate_species
    use gyro_averages, only: gyro_average
    use mp, only: sum_allreduce

    implicit none

    complex, dimension (:,vmu_lo%llim_proc:), intent (in) :: g
    complex, dimension (:), intent (out) :: phi
    integer, intent (in) :: iky, ie

    integer :: idx, iseg, ikx, iz, ia
    integer :: izl_offset
    real, dimension (nspec) :: wgt
    complex, dimension (:), allocatable :: g0

    ia = 1

    allocate (g0(vmu_lo%llim_proc:vmu_lo%ulim_alloc))

    wgt = spec%z*spec%dens_psi0
    phi = 0.

    idx = 0 ; izl_offset = 0
    iseg = 1
    ikx = ikxmod(iseg,ie,iky)
    if(zonal_mode(iky).and.abs(akx(ikx)) < epsilon(0.)) then
      phi(:) = 0.0
      return
    endif
    do iz = iz_low(iseg), iz_up(iseg)
       idx = idx + 1
        call gyro_average (g(idx,:), iky, ikx, iz, g0)
        call integrate_species (g0, iz, wgt, phi(idx),reduce_in=.false.)
    end do
    izl_offset = 1
    if (nsegments(ie,iky) > 1) then
       do iseg = 2, nsegments(ie,iky)
          ikx = ikxmod(iseg,ie,iky)
          do iz = iz_low(iseg)+izl_offset, iz_up(iseg)
             idx = idx + 1
             call gyro_average (g(idx,:), iky, ikx, iz, g0)
             call integrate_species (g0, iz, wgt, phi(idx),reduce_in=.false.)
          end do
          if (izl_offset == 0) izl_offset = 1
       end do
    end if

    call sum_allreduce(phi)

  end subroutine integrate_over_velocity

  subroutine get_fields_for_response_matrix (phi, iky, ie)

    use stella_layouts, only: vmu_lo
    use species, only: spec
    use species, only: has_electron_species
    use stella_geometry, only: dl_over_b
    use extended_zgrid, only: iz_low, iz_up
    use extended_zgrid, only: ikxmod
    use extended_zgrid, only: nsegments
    use kt_grids, only: zonal_mode, akx
    use dist_fn, only: adiabatic_option_switch
    use dist_fn, only: adiabatic_option_fieldlineavg
    use fields, only: gamtot3
    use fields_arrays, only: gamtot

    implicit none

    complex, dimension (:), intent (inout) :: phi
    integer, intent (in) :: iky, ie

    integer :: idx, iseg, ikx, iz, ia
    integer :: izl_offset
    complex, dimension (:), allocatable :: g0
    complex :: tmp

    ia = 1

    allocate (g0(vmu_lo%llim_proc:vmu_lo%ulim_alloc))

    idx = 0 ; izl_offset = 0
    iseg = 1
    ikx = ikxmod(iseg,ie,iky)
    if(zonal_mode(iky).and.abs(akx(ikx)) < epsilon(0.)) then
      phi(:) = 0.0
      return
    endif
    do iz = iz_low(iseg), iz_up(iseg)
       idx = idx + 1
       phi(idx) = phi(idx)/gamtot(iky,ikx,iz)
    end do
    izl_offset = 1
    if (nsegments(ie,iky) > 1) then
       do iseg = 2, nsegments(ie,iky)
          ikx = ikxmod(iseg,ie,iky)
          do iz = iz_low(iseg)+izl_offset, iz_up(iseg)
             idx = idx + 1
             phi(idx) = phi(idx)/gamtot(iky,ikx,iz)
          end do
          if (izl_offset == 0) izl_offset = 1
       end do
    end if

    if (.not.has_electron_species(spec) .and. &
         adiabatic_option_switch == adiabatic_option_fieldlineavg) then
       if (zonal_mode(iky)) then
          ! no connections for ky = 0
          iseg = 1
          tmp = sum(dl_over_b(ia,:)*phi)
          phi = phi + tmp*gamtot3(ikxmod(1,ie,iky),:)
       end if
    end if

    deallocate (g0)

  end subroutine get_fields_for_response_matrix

  subroutine finish_response_matrix

    use fields_arrays, only: response_matrix
#if !defined ISO_C_BINDING

    implicit none

#else
    use mpi

    implicit none

    integer :: ierr

    if (window.ne.MPI_WIN_NULL) call mpi_win_free(window,ierr)
#endif

    if (allocated(response_matrix)) deallocate (response_matrix)
    response_matrix_initialized = .false.

  end subroutine finish_response_matrix

#ifdef MPI
!----------------------------------------------------------!
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !
! !!!!!!!!!!!!!!PARALLEL LU DECOMPOSITIONS!!!!!!!!!!!!!!!! !
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !
!----------------------------------------------------------!

#ifdef ISO_C_BINDING

  !this subroutine parallelizes the LU decomposition on a single
  !node using MPIs shared memory interface
  !It also splits up jtwist the independent matrices across nodes
  !Ideal speed up: cores_per_node*min(jtwist,ncores)
  subroutine parallel_LU_decomposition_local (iky)

    use, intrinsic :: iso_c_binding, only: c_ptr, c_f_pointer
    use fields_arrays, only: response_matrix
    use mp, only: barrier, broadcast, sum_allreduce
    use mp, only: mp_comm, scope, allprocs, sharedprocs, curr_focus
    use mp, only: scrossdomprocs, sgproc0, mp_abort
    use mp, only: job, iproc, proc0, nproc, numnodes, inode
    use job_manage, only: njobs
    use extended_zgrid, only: neigen
    use mpi
    use linear_solve, only: imaxloc

    implicit none

    integer, intent (in) :: iky

    integer, dimension (:,:), allocatable :: eig_limits
    integer, dimension (:), allocatable :: job_list
    integer, dimension (:), allocatable :: row_limits
    logical, dimension (:,:), allocatable :: node_jobs

    real, parameter :: zero = 1.0e-20
    real, dimension (:), allocatable :: vv
    complex, dimension (:), allocatable :: dum
    complex, dimension (:,:), pointer :: lu

    type(c_ptr) :: bptr

    logical :: needs_send = .false.

    integer :: prior_focus, nodes_on_job
    integer :: ijob, i,j,k,ie,n, rec_len
    integer :: imax, jroot, neig, ierr, win, nroot
    integer (kind=MPI_ADDRESS_KIND) :: win_size, real_size
    integer :: rdiv, rmod
    integer :: ediv, emod
    integer :: disp_unit = 1
    real :: dmax, tmp

    prior_focus = curr_focus

    call scope (sharedprocs)

    allocate (node_jobs(numnodes,njobs)); node_jobs = .false.
    allocate (job_list(nproc)); job_list = 0
    allocate (row_limits(0:nproc))
    allocate (eig_limits(0:numnodes,njobs)); eig_limits = 0

    inquire(iolength=rec_len) zero !this will return 2 or 8 for double
    if(rec_len.eq.1.or.rec_len.eq.4) then
      real_size = 4_MPI_ADDRESS_KIND
    else if (rec_len.eq.2.or.rec_len.eq.8) then
      real_size = 8_MPI_ADDRESS_KIND
    else
      call mp_abort('failure retrieving the size of a real')
    endif

    job_list(iproc+1) = job
    call sum_allreduce(job_list)

    if(proc0) then
      do j = 1, nproc
        node_jobs(inode+1,job_list(j)+1) = .true. !create a map of which nodes have which jobs
      enddo
    endif

    !make sure all processors have this map
    call scope (allprocs)
    call mpi_allreduce &
         (MPI_IN_PLACE, node_jobs,size(node_jobs), MPI_LOGICAL, MPI_LOR, mp_comm, ierr)
    call scope (sharedprocs)

    do ijob = 0, njobs-1
      jroot=-1
      do j = 1, nproc
        if(job_list(j).eq.ijob) then
          jroot = j-1 !the first processor on this job will be the root process
          exit
        endif
      enddo

      if(jroot.eq.-1) cycle !no processors on this node are on this job

      if(iproc.eq.jroot) neig = neigen(iky)

      ! broadcast number of matrices
      call broadcast(neig,jroot)

      ! split up neig across nodes that have the current job
      nodes_on_job = count(node_jobs(:,ijob+1))
      ediv = neig/nodes_on_job
      emod = mod(neig,nodes_on_job)

      eig_limits(0,ijob+1) = 1
      do j = 1, numnodes
        if(node_jobs(j,ijob+1)) then
          eig_limits(j,ijob+1) = eig_limits(j-1,ijob+1) + ediv
          if(emod.gt.0) then
            eig_limits(j,ijob+1) = eig_limits(j,ijob+1) + 1
            emod = emod - 1
          endif
        else
          eig_limits(j,ijob+1) = eig_limits(j-1,ijob+1)
        endif
      enddo

      do ie = eig_limits(inode,ijob+1), eig_limits(inode+1,ijob+1)-1
        win_size = 0
        if(iproc.eq.jroot) then
          needs_send = .true.
          n = size(response_matrix(iky)%eigen(ie)%idx)
          win_size = int(n*n,MPI_ADDRESS_KIND)*2*real_size !complex size
        endif

        !broadcast size of matrix
        call broadcast(n, jroot)

        allocate (dum(n))
        allocate (vv(n))

        !allocate the window
        call mpi_win_allocate_shared(win_size,disp_unit,MPI_INFO_NULL,mp_comm,bptr,win,ierr)

        if(iproc.ne.jroot) then
          !make sure all the procs have the right memory address
          call mpi_win_shared_query(win, jroot, win_size, disp_unit, bptr, ierr)
        endif

        ! bind this c_ptr to our fortran matrix
        call c_f_pointer(bptr,lu,(/n,n/))

        !load the matrix
        if(iproc.eq.jroot) lu = response_matrix(iky)%eigen(ie)%zloc

        !syncronize the processors
        call mpi_win_fence(0,win,ierr)

        ! All the processors have the matrix.
        ! Now perform LU decomposition
        vv = maxval(cabs(lu),dim=2)
        if (any(vv==0.0)) &
           write (*,*) 'singular matrix in lu_decomposition on job ', job, ', process ', iproc
        vv = 1.0/vv
        do j = 1, n
           !divide up the work using row_limits
           rdiv = (n-j)/nproc
           rmod = mod(n-j,nproc)
           row_limits(0) = j+1
           if(rdiv.eq.0) then
             row_limits(rmod+1:) = -1
             do k=1,rmod
               row_limits(k)  = row_limits(k-1) + 1
             enddo
           else
             do k=1,nproc
               row_limits(k) = row_limits(k-1) + rdiv
               if(k.le.rmod) row_limits(k) = row_limits(k) + 1
             enddo
           endif

           !pivot if needed
           dmax = -1.0
           do k = j, n
             tmp = vv(k)*abs(lu(k,j))
             if(tmp.gt.dmax) then
               dmax = tmp
               imax = k
             endif
           enddo
!          imax = (j-1) + imaxloc(vv(j:n)*cabs(lu(j:n,j)))

           if(iproc.eq.jroot) then
             response_matrix(iky)%eigen(ie)%idx(j) = imax
             if (j /= imax) then
               dum = lu(imax,:)
               lu(imax,:) = lu(j,:)
               lu(j,:) = dum
               vv(imax) = vv(j)
             end if
             if (lu(j,j)==0.0) lu(j,j) = zero
           else
             if (j /= imax) vv(imax) = vv(j)
           endif

           call mpi_win_fence(0,win,ierr)

           !get the lead multiplier
           do i = row_limits(iproc), row_limits(iproc+1)-1
             lu(i,j) = lu(i,j)/lu(j,j)
           enddo

           call mpi_win_fence(0,win,ierr)

           do k=row_limits(iproc), row_limits(iproc+1)-1
             do i = j+1,n
                lu(i,k) = lu(i,k) - lu(i,j)*lu(j,k)
             enddo
           enddo

           call mpi_win_fence(0,win,ierr)
        enddo
        !LU decomposition ends here

        !copy the decomposed matrix over
        if(iproc.eq.jroot) response_matrix(iky)%eigen(ie)%zloc = lu

        call mpi_win_free(win,ierr)
        deallocate (vv,dum)
      enddo
    enddo

    call scope(scrossdomprocs)

    !copy all the matrices across all nodes
    if(sgproc0) then
      do ie = 1, neigen(iky)
        nroot = 0
        if(needs_send.and. &
          (ie.ge.eig_limits(inode,job+1).and.ie.lt.eig_limits(inode+1,job+1))) nroot = iproc
        !first let processors know who is sending the data
        call sum_allreduce(nroot)
        !now send the data
        call broadcast(response_matrix(iky)%eigen(ie)%zloc,nroot)
        call broadcast(response_matrix(iky)%eigen(ie)%idx, nroot)
      enddo
    endif

    call scope(prior_focus)

    deallocate (node_jobs,job_list,row_limits,eig_limits)
  end subroutine parallel_LU_decomposition_local

#endif /* ISO_C_BINDING */

  !this subroutine parallelizes the LU decomposition across
  !all cores. Ideal speed up: ncores
  subroutine parallel_LU_decomposition_global (iky)

    use fields_arrays, only: response_matrix
    use mp, only: barrier, broadcast, sum_allreduce
    use mp, only: mp_comm, scope, allprocs, sharedprocs, curr_focus
    use mp, only: job, iproc, proc0, nproc, mpicmplx
#ifdef ISO_C_BINDING
    use mp, only: sgproc0, scrossdomprocs
#endif
    use job_manage, only: njobs
    use extended_zgrid, only: neigen
    use mpi
    use linear_solve, only: imaxloc

    implicit none

    integer, intent (in) :: iky

    integer, dimension (:), allocatable :: job_roots,  eig_roots
    integer, dimension (:), allocatable :: row_limits, eig_limits
    integer, dimension (MPI_STATUS_SIZE) :: status

    real, parameter :: zero = 1.0e-20
    real, dimension (:), allocatable :: vv
    complex, dimension (:), allocatable :: dum
    integer, dimension (:), allocatable :: idx
    complex, dimension (:,:), allocatable :: lu

    integer :: sproc
    logical :: sproc0

    integer :: eig_comm, ceig_comm !c for 'cross'
    integer :: ieig_core, ceig_core, eig_cores
    integer :: ncomm
    integer :: prior_focus
    integer :: ie, ie_hi, r_lo, r_hi
    integer :: ijob, i,j,k,n, n_send, rsize
    integer :: imax, neig, ierr
    integer :: istage, nstage
    integer :: rdiv, rmod
    integer :: ediv, emod

    real :: dmax, tmp

    prior_focus = curr_focus

    sproc  = iproc
    sproc0 = proc0

    call scope (allprocs)

    allocate (job_roots(0:njobs-1)); job_roots = 0

    if(sproc0) job_roots(job) = iproc

    call sum_allreduce(job_roots)

    do ijob = 0, njobs-1
      if(job.eq.ijob.and.sproc0) then
        neig = neigen(iky)
      endif

      ! broadcast number of matrices for this job
      call broadcast(neig,job_roots(ijob))

      !set up communicator for cores working on a single matrix
      call mpi_comm_split(mp_comm,mod(iproc,neig),iproc,eig_comm,ierr)
      call mpi_comm_size(eig_comm,eig_cores,ierr)
      call mpi_comm_rank(eig_comm,ieig_core,ierr)

      !set up a communicator that crosses the previous one
      call mpi_comm_split(mp_comm,ieig_core,iproc,ceig_comm,ierr)
      call mpi_comm_rank(ceig_comm,ceig_core,ierr)

      call mpi_bcast(ceig_core,1,MPI_INT,0,eig_comm,ierr)

      ncomm = min(neig,nproc) !number of communicators

      allocate (eig_roots(0:ncomm-1)); eig_roots = 0
      allocate (eig_limits(0:ncomm))
      allocate (row_limits(0:eig_cores))

      if(ieig_core.eq.0) eig_roots(ceig_core) = iproc

      call sum_allreduce(eig_roots)

      ! split up neigen across cores
      ediv = neig/ncomm
      emod = mod(neig,ncomm)

      !how many stages will the LU decomposition take?
      nstage = ediv
      if(emod.gt.0) nstage = nstage + 1

      !determine which parts of neigen this communicator processes
      eig_limits(0) = 1
      do j = 1, ncomm
        eig_limits(j) = eig_limits(j-1) + ediv
        if(j.le.emod) then
          eig_limits(j) = eig_limits(j) + 1
        endif
      enddo

      do istage = 0, nstage-1
        !transfer the data from job root to root of subcommunicator
        do j = 0,ncomm-1
          ie    = eig_limits(j) + istage
          ie_hi = eig_limits(j+1) - 1
          if(ie.gt.ie_hi) cycle

          if(iproc.eq.job_roots(ijob).and.iproc.eq.eig_roots(j)) then !no need for data transfer
            n=size(response_matrix(iky)%eigen(ie)%idx)
            allocate(lu(n,n))
            lu=response_matrix(iky)%eigen(ie)%zloc
          else if(iproc.eq.job_roots(ijob)) then !send data to subroots
            !send size of matrix
            n_send=size(response_matrix(iky)%eigen(ie)%idx)
            call mpi_send(n_send,1,MPI_INT,eig_roots(j),j,mp_comm,ierr)
            !send matrix
            call mpi_send(response_matrix(iky)%eigen(ie)%zloc, &
                          n_send*n_send,mpicmplx,eig_roots(j),nproc+j,mp_comm,ierr)
          else if(iproc.eq.eig_roots(j)) then !subroot gets the data
            !receive size of matrix
            call mpi_recv(n,1,MPI_INT,job_roots(ijob),j,mp_comm,status, ierr)
            allocate(lu(n,n))
            !receive matrix
            call mpi_recv(lu,n*n,mpicmplx,job_roots(ijob),nproc+j,mp_comm,status,ierr)
          endif
        enddo

        if(istage.ge.(eig_limits(ceig_core+1)-eig_limits(ceig_core))) cycle !nothing for this communicator to do

        !broadcast matrix and its size across the communicator
        call mpi_bcast(n,1,MPI_INT,0,eig_comm,ierr)
        if(.not.allocated(lu)) allocate (lu(n,n))
        if(.not.allocated(vv)) allocate (vv(n))

        call mpi_bcast(lu,n*n,mpicmplx,0,eig_comm,ierr)

        allocate (dum(n))
        allocate (idx(n))

        ! All the processors have the matrix.
        ! Now perform LU decomposition
        vv = maxval(cabs(lu),dim=2)
        if (any(vv==0.0)) &
          write (*,*) 'singular matrix in lu_decomposition on job ', job, ', process ', iproc
        vv = 1.0/vv
        do j = 1, n
          !divide up the work using row_limits
          rdiv = (n-j)/eig_cores
          rmod = mod(n-j,eig_cores)
          row_limits(0) = j+1
          if(rdiv.eq.0) then
            row_limits(rmod+1:) = -1
            do k=1,rmod
              row_limits(k)  = row_limits(k-1) + 1
            enddo
          else
            do k=1,eig_cores
              row_limits(k) = row_limits(k-1) + rdiv
              if(k.le.rmod) row_limits(k) = row_limits(k) + 1
            enddo
          endif

          !pivot if needed
          dmax = -1.0
          do k = j, n
            tmp = vv(k)*abs(lu(k,j))
            if(tmp.gt.dmax) then
              dmax = tmp
              imax = k
            endif
          enddo
!         imax = (j-1) + imaxloc(vv(j:n)*cabs(lu(j:n,j)))
          if (j /= imax) then
            dum = lu(imax,:)
            lu(imax,:) = lu(j,:)
            lu(j,:) = dum
            vv(imax) = vv(j)
          end if
          if(ieig_core.eq.0) idx(j) = imax

           !get the lead multiplier
           if (lu(j,j)==0.0) lu(j,j) = zero
           do i = j+1, n
             lu(i,j) = lu(i,j)/lu(j,j)
           enddo

           r_lo = row_limits(ieig_core)
           r_hi = row_limits(ieig_core+1)-1

           do k=r_lo, r_hi
             do i = j+1,n
                lu(i,k) = lu(i,k) - lu(i,j)*lu(j,k)
             enddo
           enddo

           do i=0,eig_cores-1
             r_lo = row_limits(i)
             r_hi = row_limits(i+1)-1
             rsize = (r_hi-r_lo+1)*(n-j)
             if(r_lo.gt.r_hi) cycle
             !call mpi_bcast(lu(j+1:n,r_lo:r_hi),rsize,mpicmplx,i,eig_comm,ierr)
             do k=r_lo,r_hi
               call mpi_bcast(lu(j+1:n,k),n-j,mpicmplx,i,eig_comm,ierr)
             enddo
           enddo
        enddo
        !LU decomposition ends here

        !copy the decomposed matrix over
        do j = 0,ncomm-1

          ie = eig_limits(j) + istage
          ie_hi = eig_limits(j+1) - 1
          if(ie.gt.ie_hi) cycle

          if(iproc.eq.job_roots(ijob).and.iproc.eq.eig_roots(j)) then !no need for data transfer
            response_matrix(iky)%eigen(ie)%zloc = lu
            response_matrix(iky)%eigen(ie)%idx  = idx
          else if(iproc.eq.eig_roots(j)) then !subroot sends the data
            !send indices
            call mpi_send(idx,n,MPI_INT,job_roots(ijob),j,mp_comm,ierr)
            !send matrix
            call mpi_send(lu,n*n,mpicmplx,job_roots(ijob),nproc+j,mp_comm,ierr)
          else if(iproc.eq.job_roots(ijob)) then !receive data from subroot
            !receive indices
            call mpi_recv(response_matrix(iky)%eigen(ie)%idx, &
                          n,MPI_INT,eig_roots(j),j,mp_comm,status,ierr)
            !receive matrix
            call mpi_recv(response_matrix(iky)%eigen(ie)%zloc, &
                          n*n,mpicmplx,eig_roots(j),nproc+j,mp_comm,status,ierr)
          endif
        enddo
        deallocate (vv, lu, idx, dum)
      enddo
      deallocate (eig_roots, eig_limits, row_limits)
    enddo

#ifdef ISO_C_BINDING
    if(sgproc0) then
      call scope(scrossdomprocs)
      !copy all the matrices across all nodes
      do ie = 1, neigen(iky)
        call broadcast(response_matrix(iky)%eigen(ie)%zloc)
        call broadcast(response_matrix(iky)%eigen(ie)%idx)
      enddo
    endif

    call scope(prior_focus)
#else
    call scope(prior_focus)

    !copy all the matrices across all nodes
    do ie = 1, neigen(iky)
      call broadcast(response_matrix(iky)%eigen(ie)%zloc)
      call broadcast(response_matrix(iky)%eigen(ie)%idx)
    enddo
#endif
    deallocate (job_roots)
  end subroutine parallel_LU_decomposition_global

#endif /* MPI */

end module response_matrix

ext: fpp


text:
  !=========================================================================================
  !Copyright (c) 2009-2015, The Regents of the University of Massachusetts, Amherst.
  !E. Polizzi research lab
  !All rights reserved.
  !
  !Redistribution and use in source and binary forms, with or without modification,
  !are permitted provided that the following conditions are met:
  !
  !1. Redistributions of source code must retain the above copyright notice, this list of conditions
  !   and the following disclaimer.
  !2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
  !   and the following disclaimer in the documentation and/or other materials provided with the distribution.
  !3. Neither the name of the University nor the names of its contributors may be used to endorse or promote
  !    products derived from this software without specific prior written permission.
  !
  !THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
  !BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  !ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  !EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  !SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  !LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
  !IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  !==========================================================================================

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!! SPARSE PRIMITIVES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!! DOUBLE PRECISION VERSION

!!!!! List of routines

!daddcsr, zaddcsr           :: add two square csr matrices and return a new one (C=alpha*A+beta*B)
!zdaddcsr                   :: same as zaddcsr but A, B are real;  alpha, beta, C are complex

!dcsr2csr_up,  zcsr2csr_up  :: extract upper triangular part of a full csr matrix
!dcsr2csr_low, zcsr2csr_low :: extract lower triangular part of a full csr matrix

!dcsrmm, zcsrmm             :: mat-vec multiplication =>  b=alpha*A*x+beta*b
!zhcsrmm                    :: mat-vec multiplication =>  b=alpha*A*x+beta*b (A is complex Hermitian if UPLO=L,U)
!dzcsrmm                    :: same as zcsrmm but A is real;  alpha, beta, b are complex

!dcsr_transpose, zcsr_transpose :: transpose a csr matrix
!zcsr_htranspose                :: transpose conjugate a csr matrix

!dcsr_uplo_to_csr, zcsr_uplo_to_csr  :: convert a upper or lower csr matrix into a full csr one
!zhcsr_uplo_to_csr                   :: convert a upper or lower (Hermitian) csr matrix into a full csr one

subroutine daddcsr(N,opt,alpha,sa,isa,jsa,beta,sb,isb,jsb,sc,isc,jsc)
  !  Purpose
  !  =======
  !
  !  Addition of two (square) sparse CSR matrices A and B and return a new (square) CSR matrix C
  !  C=alpha*A+beta*B
  !
  !  A,B,C are real double precision
  !
  !  Arguments
  !  =========
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,B,C.  N >= 0.
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs addition, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine computes only values of the array isc of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value ic(n+1) - 1 is the actual
  !         number of the elements in the arrays sc and jsc.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array jsc
  !
  !         if opt=3, the routine has been called previously with the parameter opt=2. On exit it will
  !         return the array sc
  !
  !
  !  alpha  (input) DOUBLE PRECISION
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa DOUBLE PRECISION)
  !
  !  beta  (input) DOUBLE PRECISION
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb DOUBLE PRECISION)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc DOUBLE PRECISION).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  double precision,dimension(*) :: sa,sb,sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  double precision :: alpha,beta
  integer :: opt,N

  integer :: i,ka,kb,kc,ja,jb
  logical :: next,row,col,val

  row=.false.
  col=.false.
  val=.false.

  if (opt==0) then
     row=.true.
     col=.true.
     val=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     col=.true.
  elseif (opt==3) then
     val=.true.
  endif

  kc=1
  if (row) isc(1)=kc
  do i=1,N
     ka=isa(i)
     kb=isb(i)

     next=.false.
     do while (.not.next)

        if (ka<=isa(i+1)-1) then
           ja=jsa(ka)
        else
           ja=n+1
        endif

        if (kb<=isb(i+1)-1) then
           jb=jsb(kb)
        else
           jb=n+1
        endif

        if (ja==jb) then
           if (val) sc(kc)=alpha*sa(ka)+beta*sb(kb)
           if (col) jsc(kc)=ja
           ka=ka+1
           kb=kb+1
           kc=kc+1
        elseif (ja<jb) then
           if (val) sc(kc)=alpha*sa(ka)
           if (col) jsc(kc)=ja
           ka=ka+1
           kc=kc+1
        elseif (jb<ja) then
           if (val) sc(kc)=beta*sb(kb)
           if (col) jsc(kc)=jb
           kb=kb+1
           kc=kc+1
        end if

        if ((ka>isa(i+1)-1).and.(kb>isb(i+1)-1)) next=.true.

     end do

     if (row) isc(i+1)=kc

  end do

end subroutine daddcsr

subroutine zaddcsr(N,opt,alpha,sa,isa,jsa,beta,sb,isb,jsb,sc,isc,jsc)
  !  Purpose
  !  =======
  !
  !  Addition of two (square) sparse CSR matrices A and B and return a new (square) CSR matrix C
  !  C=alpha*A+beta*B
  !
  !  A,B,C are complex double precision
  !
  !  Arguments
  !  =========
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,B,C.  N >= 0.
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs addition, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine computes only values of the array isc of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value ic(n+1) - 1 is the actual
  !         number of the elements in the arrays sc and jsc.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array jsc
  !
  !         if opt=3, the routine has been called previously with the parameter opt=2. On exit it will
  !         return the array sc
  !
  !
  !  alpha  (input) COMPLEX DOUBLE PRECISION
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa COMPLEX DOUBLE PRECISION)
  !
  !  beta  (input) COMPLEX DOUBLE PRECISION
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb COMPLEX DOUBLE PRECISION)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc COMPLEX DOUBLE PRECISION).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  complex(kind=kind(1.0d0)),dimension(*) :: sa,sb,sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  complex(kind=kind(1.0d0)) :: alpha,beta
  integer :: opt,N

  integer :: i,ka,kb,kc,ja,jb
  logical :: next,row,col,val

  row=.false.
  col=.false.
  val=.false.

  if (opt==0) then
     row=.true.
     col=.true.
     val=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     col=.true.
  elseif (opt==3) then
     val=.true.
  endif

  kc=1
  if (row) isc(1)=kc

  do i=1,N
     ka=isa(i)
     kb=isb(i)

     next=.false.
     do while (.not.next)

        if (ka<=isa(i+1)-1) then
           ja=jsa(ka)
        else
           ja=n+1
        endif

        if (kb<=isb(i+1)-1) then
           jb=jsb(kb)
        else
           jb=n+1
        endif

        if (ja==jb) then
           if (val) sc(kc)=alpha*sa(ka)+beta*sb(kb)
           if (col) jsc(kc)=ja
           ka=ka+1
           kb=kb+1
           kc=kc+1
        elseif (ja<jb) then
           if (val) sc(kc)=alpha*sa(ka)
           if (col) jsc(kc)=ja
           ka=ka+1
           kc=kc+1
        elseif (jb<ja) then
           if (val) sc(kc)=beta*sb(kb)
           if (col) jsc(kc)=jb
           kb=kb+1
           kc=kc+1
        end if

        if ((ka>isa(i+1)-1).and.(kb>isb(i+1)-1)) next=.true.

     end do

     if (row) isc(i+1)=kc

  end do

end subroutine zaddcsr

subroutine zdaddcsr(N,opt,alpha,sa,isa,jsa,beta,sb,isb,jsb,sc,isc,jsc)
  !  Purpose
  !  =======
  !
  !  Addition of two (square) sparse CSR matrices A and B and return a new (square) CSR matrix C
  !  C=alpha*A+beta*B
  !
  !  C is complex double precision, A and B are real double precision, (alpha,beta are COMPLEX)
  !
  !  Arguments
  !  =========
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,B,C.  N >= 0.
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs addition, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine computes only values of the array isc of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value ic(n+1) - 1 is the actual
  !         number of the elements in the arrays sc and jsc.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array jsc
  !
  !         if opt=3, the routine has been called previously with the parameter opt=2. On exit it will
  !         return the array sc
  !
  !
  !  alpha  (input) COMPLEX DOUBLE PRECISION
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa DOUBLE PRECISION)
  !
  !  beta  (input) COMPLEX DOUBLE PRECISION
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb DOUBLE PRECISION)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc COMPLEX DOUBLE PRECISION).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  double precision,dimension(*) :: sa,sb
  complex(kind=kind(1.0d0)),dimension(*) :: sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  complex(kind=kind(1.0d0)) :: alpha,beta
  integer :: opt,N
  integer :: i,ka,kb,kc,ja,jb
  logical :: next,row,col,val

  row=.false.
  col=.false.
  val=.false.

  if (opt==0) then
     row=.true.
     col=.true.
     val=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     col=.true.
  elseif (opt==3) then
     val=.true.
  endif

  kc=1
  if (row) isc(1)=kc

  do i=1,N
     ka=isa(i)
     kb=isb(i)

     next=.false.
     do while (.not.next)

        if (ka<=isa(i+1)-1) then
           ja=jsa(ka)
        else
           ja=n+1
        endif

        if (kb<=isb(i+1)-1) then
           jb=jsb(kb)
        else
           jb=n+1
        endif

        if (ja==jb) then
           if (val) sc(kc)=alpha*sa(ka)+beta*sb(kb)
           if (col) jsc(kc)=ja
           ka=ka+1
           kb=kb+1
           kc=kc+1
        elseif (ja<jb) then
           if (val) sc(kc)=alpha*sa(ka)
           if (col) jsc(kc)=ja
           ka=ka+1
           kc=kc+1
        elseif (jb<ja) then
           if (val) sc(kc)=beta*sb(kb)
           if (col) jsc(kc)=jb
           kb=kb+1
           kc=kc+1
        end if

        if ((ka>isa(i+1)-1).and.(kb>isb(i+1)-1)) next=.true.

     end do

     if (row) isc(i+1)=kc

  end do

end subroutine zdaddcsr

subroutine dcsr2csr_up(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the upper triangular part (including diagonal) of a real CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr upper format may be overestimated
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  double precision,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)>=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine dcsr2csr_up

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine zcsr2csr_up(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the upper triangular part (including diagonal) of a real CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr upper format may be overestimated
  !
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)>=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine zcsr2csr_up

subroutine dcsr2csr_low(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the lower triangular part (including diagonal) of a real CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr lower format may be overestimated
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  double precision,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)<=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine dcsr2csr_low
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine zcsr2csr_low(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the lower triangular part (including diagonal) of a COMPLEX CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr lower format may be overestimated
  !
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)<=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine zcsr2csr_low

SUBROUTINE dcsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      REAL DOUBLE PRECISION version  !
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) DOUBLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE PRECISION)
  !
  !  X     (input) DOUBLE PRECISION
  !        matrix of size (Mxrhs)
  !
  !  beta  (input) DOUBLE PRECISION
  !
  !  B     (input/output) DOUBLE PRECISION
  !        matrix of size (Nxrhs). On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  double precision :: alpha,beta
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  double precision,dimension(M,*):: x
  double precision,dimension(N,*) ::b
!!!!!!!!

  integer ::i,j,k
  double precision,parameter :: DZERO=0.0d0

!!!!!!!! Initialization

if (beta/=DZERO) then
call DSCAL(N*rhs,b,beta,1)
else
b(1:N,1:rhs)=DZERO
end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
           end do
        end do

     elseif (TRANS=='T') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

     end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  elseif ((UPLO=='L').or.(UPLO=='U')) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
              if (j/=i) b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

  end if

end SUBROUTINE dcsrmm

SUBROUTINE zcsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      COMPLEX DOUBLE PRECISION version
  !
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A   -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !       if TRANS='C' o(A)=A^H -- transpose conjugate mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) DOUBLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE PRECISION)
  !
  !  X     (input) DOUBLE PRECISION
  !        matrix of size (Mxrhs)
  !
  !  beta  (input) DOUBLE PRECISION
  !
  !  B     (input/output) DOUBLE PRECISION
  !        matrix of size (Nxrhs). On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  complex (kind=kind(1.0d0)) :: alpha,beta
  complex (kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex (kind=kind(1.0d0)),dimension(M,*):: x
  complex (kind=kind(1.0d0)),dimension(N,*) ::b
!!!!!!!!

  integer ::i,j,k
  complex (kind=kind(1.0d0)),parameter :: ZEROC=(0.0d0,0.0d0)

!!!!!!!! Initialization

if (beta/=ZEROC) then
call ZSCAL(N*rhs,b,beta,1)
else
b(1:N,1:rhs)=ZEROC
end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
           end do
        end do

     elseif (TRANS=='T') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

     elseif (TRANS=='C') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*conjg(a(k))*x(i,1:rhs)
           end do
        end do

     end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  elseif ((UPLO=='L').or.(UPLO=='U')) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

if ((TRANS=='N').or.(TRANS=='T')) then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
              if (j/=i) b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

 elseif (TRANS=='C') then

         do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*conjg(a(k))*x(i,1:rhs)
              if (j/=i) b(i,1:rhs)=b(i,1:rhs)+alpha*conjg(a(k))*x(j,1:rhs)
           end do
        end do

     end if

  end if

END SUBROUTINE zcsrmm

SUBROUTINE zhcsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      Note: A is complex Hermitian if UPLO='L or UPLO='U'
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric/hermitian CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric/hermitian CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A   -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !       if TRANS='C' o(A)=A^H -- transpose conjugate mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) COMPLEX DOUBLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX DOUBLE PRECISION)
  !
  !  X     (input) COMPLEX DOUBLE PRECISION
  !        matrix of size (Mxrhs)
  !
  !  beta  (input) COMPLEX DOUBLE PRECISION
  !
  !  B     (input/output) COMPLEX DOUBLE PRECISION
  !        matrix of size (Nxrhs). On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  complex(kind=kind(1.0d0)) :: alpha,beta
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),dimension(M,*) :: x
  complex(kind=kind(1.0d0)),dimension(N,*) :: b
!!!!!!!!

  integer ::i,j,k
  double precision,parameter :: DZERO=0.0d0
  complex(kind=kind(1.0d0)),parameter :: ZEROC=(DZERO,DZERO)

!!!!!!!! Initialization

if (beta/=ZEROC) then
call ZSCAL(N*rhs,b,beta,1)
else
b(1:N,1:rhs)=ZEROC
end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
           end do
        end do

     elseif (TRANS=='T') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

     elseif (TRANS=='C') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*conjg(a(k))*x(i,1:rhs)
           end do
        end do

     end if

  elseif ((UPLO=='L').or.(UPLO=='U')) then

if ((TRANS=='N').or.(TRANS=='C')) then

     do i=1,N
        do k=ia(i),ia(i+1)-1
           j=ja(k)
           b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
   if (j/=i) b(j,1:rhs)=b(j,1:rhs)+alpha*conjg(a(k))*x(i,1:rhs)
        end do
     end do

elseif (TRANS=='T') then

         do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*conjg(a(k))*x(i,1:rhs)
              if (j/=i) b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
           end do
        end do

     end if

  end if

end SUBROUTINE zhcsrmm

SUBROUTINE dzcsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      DOUBLE PRECISION version
  !
  !      Note: A is real
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A   -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !       if TRANS='C' o(A)=A^H -- transpose conjugate mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) DOUBLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE PRECISION)
  !
  !  X     (input) DOUBLE PRECISION
  !        matrix of size (Mxrhs)
  !
  !  beta  (input) DOUBLE PRECISION
  !
  !  B     (input/output) DOUBLE PRECISION
  !        matrix of size (Nxrhs). On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  complex (kind=kind(1.0d0)) :: alpha,beta
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex (kind=kind(1.0d0)),dimension(M,*):: x
  complex (kind=kind(1.0d0)),dimension(N,*) ::b
!!!!!!!!

  integer ::i,j,k
  complex (kind=kind(1.0d0)),parameter :: ZEROC=(0.0d0,0.0d0)

!!!!!!!! Initialization

if (beta/=ZEROC) then
call ZSCAL(N*rhs,b,beta,1)
else
b(1:N,1:rhs)=ZEROC
end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
           end do
        end do

     elseif ((TRANS=='T').or.(TRANS=='C')) then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

     end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  elseif ((UPLO=='L').or.(UPLO=='U')) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
              if (j/=i) b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

  end if

END SUBROUTINE dzcsrmm

subroutine dcsr_transpose(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine transpose a double real CSR matrix A
  !      into a new CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  double precision,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
  integer :: i,j,k,istart
  integer,dimension(1:N) :: tab
!!!!

  !! #of elements by row
  tab=0
  istart=N
  do i=1,N
     do j=ia(i),ia(i+1)-1
        tab(ja(j))=tab(ja(j))+1
        if (ja(j)<istart) istart=ja(j)
     end do
  end do

  !! find sia
  do i=1,istart-1
     sia(i)=0
  end do
  sia(istart)=1
  do i=istart+1,N+1
     sia(i)=sia(i-1)+tab(i-1)
  end do

  tab=sia(1:N)
  !! find jsa,sa
  do i=1,N
     do k=ia(i),ia(i+1)-1
        sja(tab(ja(k)))=i
        sa(tab(ja(k)))=a(k)
        tab(ja(k))=tab(ja(k))+1
     end do
  end do

end subroutine dcsr_transpose

subroutine zcsr_transpose(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine performs the transpose conjugate of
  !      a double complex CSR matrix A into a new CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k,istart
  integer,dimension(1:N) :: tab
!!!!

  !! #of elements by row
  tab=0
  istart=N
  do i=1,N
     do j=ia(i),ia(i+1)-1
        tab(ja(j))=tab(ja(j))+1
        if (ja(j)<istart) istart=ja(j)
     end do
  end do

  !! find sia
  do i=1,istart-1
     sia(i)=0
  end do
  sia(istart)=1
  do i=istart+1,N+1
     sia(i)=sia(i-1)+tab(i-1)
  end do

  tab=sia(1:N)
  !! find jsa,sa
  do i=1,N
     do k=ia(i),ia(i+1)-1
        sja(tab(ja(k)))=i
  !      if (ja(k)/=i) then
           sa(tab(ja(k)))=a(k)
  !      else
  !         sa(tab(ja(k)))=a(k)
  !      end if
        tab(ja(k))=tab(ja(k))+1
     end do
  end do

end subroutine zcsr_transpose

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine zcsr_htranspose(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine performs the transpose conjugate of
  !      a double complex CSR matrix A into a new CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k,istart
  integer,dimension(1:N) :: tab
!!!!

  !! #of elements by row
  tab=0
  istart=N
  do i=1,N
     do j=ia(i),ia(i+1)-1
        tab(ja(j))=tab(ja(j))+1
        if (ja(j)<istart) istart=ja(j)
     end do
  end do

  !! find sia
  do i=1,istart-1
     sia(i)=0
  end do
  sia(istart)=1
  do i=istart+1,N+1
     sia(i)=sia(i-1)+tab(i-1)
  end do

  tab=sia(1:N)
  !! find jsa,sa
  do i=1,N
     do k=ia(i),ia(i+1)-1
        sja(tab(ja(k)))=i
  !      if (ja(k)/=i) then
           sa(tab(ja(k)))=conjg(a(k))
  !      else
  !         sa(tab(ja(k)))=a(k)
  !      end if
        tab(ja(k))=tab(ja(k))+1
     end do
  end do

end subroutine zcsr_htranspose

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine dcsr_uplo_to_csr(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a double real upper or lower CSR matrix A
  !      into a new full CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  include 'f90_noruntime_interface.fi'
  integer :: N
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  double precision,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  double precision,dimension(:),pointer ::sb
  integer,dimension(:),pointer :: sib
  integer,dimension(:),pointer :: sjb
  double precision,parameter :: DONE=1.0d0,DZERO=0.0d0
!!!!
  nnza=ia(n+1)-1

  call wallocate_1d(sb, nnza, info)
  call wallocate_1i(sib, N+1, info)
  call wallocate_1i(sjb, nnza, info)

!!! transpose
  call dcsr_transpose(N,a,ia,ja,sb,sib,sjb)
!!! nullify diagonal
  do i=1,N
     do k=sib(i),sib(i+1)-1
        if (sjb(k)==i) sb(k)=DZERO
     end do
  enddo
  !! add the two matrices to form the full matrix
  opt=0
  call daddcsr(N,opt,DONE,a,ia,ja,DONE,sb,sib,sjb,sa,sia,sja)

  call wdeallocate_1d(sb)
  call wdeallocate_1i(sib)
  call wdeallocate_1i(sjb)

end subroutine dcsr_uplo_to_csr

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine zcsr_uplo_to_csr(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a complex symetric upper or lower CSR matrix A
  !      into a new full CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  include 'f90_noruntime_interface.fi'
  integer :: N
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  complex(kind=kind(1.0d0)),dimension(:),pointer ::sb
  integer,dimension(:),pointer :: sib
  integer,dimension(:),pointer :: sjb
  double precision,parameter :: DONE=1.0d0,DZERO=0.0d0
  complex(kind=kind(1.0d0)),parameter :: ONEC=(DONE,DZERO),ZEROC=(DZERO,DZERO)
!!!!
  nnza=ia(n+1)-1

  call wallocate_1z(sb, nnza, info)
  call wallocate_1i(sib, N+1, info)
  call wallocate_1i(sjb, nnza, info)

!!! transpose
  call zcsr_transpose(N,a,ia,ja,sb,sib,sjb)
!!! nullify diagonal
  do i=1,N
     do k=sib(i),sib(i+1)-1
        if (sjb(k)==i) sb(k)=ZEROC
     end do
  enddo
  !! add the two matrices to form the full matrix
  opt=0
  call zaddcsr(N,opt,ONEC,a,ia,ja,ONEC,sb,sib,sjb,sa,sia,sja)

  call wdeallocate_1z(sb)
  call wdeallocate_1i(sib)
  call wdeallocate_1i(sjb)

end subroutine zcsr_uplo_to_csr

subroutine zhcsr_uplo_to_csr(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a complex Hermitian upper or lower CSR matrix A
  !      into a new full CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  include 'f90_noruntime_interface.fi'
  integer :: N
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  complex(kind=kind(1.0d0)),dimension(:),pointer ::sb
  integer,dimension(:),pointer :: sib
  integer,dimension(:),pointer :: sjb
  double precision,parameter :: DONE=1.0d0,DZERO=0.0d0
  complex(kind=kind(1.0d0)),parameter :: ONEC=(DONE,DZERO),ZEROC=(DZERO,DZERO)
!!!!
  nnza=ia(n+1)-1

  call wallocate_1z(sb, nnza, info)
  call wallocate_1i(sib, N+1, info)
  call wallocate_1i(sjb, nnza, info)

!!! transpose
  call zcsr_htranspose(N,a,ia,ja,sb,sib,sjb)
!!! nullify diagonal
  do i=1,N
     do k=sib(i),sib(i+1)-1
        if (sjb(k)==i) sb(k)=ZEROC
     end do
  enddo
  !! add the two matrices to form the full matrix
  opt=0
  call zaddcsr(N,opt,ONEC,a,ia,ja,ONEC,sb,sib,sjb,sa,sia,sja)

  call wdeallocate_1z(sb)
  call wdeallocate_1i(sib)
  call wdeallocate_1i(sjb)

end subroutine zhcsr_uplo_to_csr

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

max_stars_repo_name: lcirrott/Kratos

max_issues_repo_name: PFEM/feast

max_forks_repo_name: PFEM/feast


text:
C----------------------------------------------------------------------
        subroutine ftghdn(iunit,hdunum)

C       return the number of the current header data unit.  The
C       first HDU (the primary array) is number 1.

C       iunit   i  fortran unit number
C       hdunum  i  returned number of the current HDU
C
C       written by Wm Pence, HEASARC/GSFC, March, 1993

        integer iunit,hdunum

C       COMMON BLOCK DEFINITIONS:--------------------------------------------
        integer nb,ne
        parameter (nb = 20)
        parameter (ne = 200)
        integer bufnum,chdu,hdutyp,maxhdu,hdstrt,hdend,nxthdr,dtstrt
        integer nxtfld
        logical wrmode
        common/ft0001/bufnum(199),chdu(nb),hdutyp(nb),maxhdu(nb),
     &  wrmode(nb),hdstrt(nb,ne),hdend(nb),nxthdr(nb),dtstrt(nb),nxtfld
C       END OF COMMON BLOCK DEFINITIONS-----------------------------------

        hdunum=chdu(bufnum(iunit))
        end

ext: f

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf