text:
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X   libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X   Portions of this code were written by
! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X   Copyright 2006-2010.
! H0 X
! H0 X   These portions of the source code are released under the GNU General
! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X   If you would like to license the source code under different terms,
! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X   Portions of this code were written by Noam Bernstein as part of
! H0 X   his employment for the U.S. Government, and are not subject
! H0 X   to copyright in the USA.
! H0 X
! H0 X
! H0 X   When using this software, please cite the following reference:
! H0 X
! H0 X   http://www.libatoms.org
! H0 X
! H0 X  Additional contributions by
! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

!X
!X RS_SparseMatrix module
!X
!% Module for sparce matrix math.
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

module RS_SparseMatrix_module

use System_module
use MPI_Context_module
use linearalgebra_module
use Atoms_module
use ScaLAPACK_module
use Matrix_module
implicit none
private

public :: RS_SparseMatrixL
type RS_SparseMatrixL
  integer :: N = 0, N_dense_rows = 0
  integer :: n_blocks = 0, data_size = 0
  integer, allocatable :: block_size(:), dense_row_of_row(:)
  integer, allocatable :: row_indices(:), data_ptrs(:)
  integer, allocatable :: col(:)

  integer :: cur_row = 1, cur_col_offset = 0
end type RS_SparseMatrixL

public :: RS_SparseMatrixD
type RS_SparseMatrixD
  type(RS_SparseMatrixL) :: l
  real(dp), allocatable :: data(:)
end type RS_SParseMatrixD

public :: RS_SparseMatrixZ
type RS_SparseMatrixZ
  type(RS_SparseMatrixL) :: l
  complex(dp), allocatable :: data(:)
end type RS_SParseMatrixZ

public :: Initialise
interface Initialise
  module procedure RS_SparseMatrixL_Initialise_at, RS_SparseMatrixL_Initialise_prod
  module procedure RS_SparseMatrixD_Initialise_at, RS_SparseMatrixD_Initialise_prod
  module procedure RS_SparseMatrixZ_Initialise_at, RS_SparseMatrixZ_Initialise_prod
end interface Initialise

public :: Finalise
interface Finalise
  module procedure RS_SparseMatrixL_Finalise
  module procedure RS_SparseMatrixD_Finalise
  module procedure RS_SparseMatrixZ_Finalise
end interface Finalise

public :: Wipe
interface Wipe
  module procedure RS_SparseMatrixL_Wipe
  module procedure RS_SparseMatrixD_Wipe
  module procedure RS_SparseMatrixZ_Wipe
end interface Wipe

public :: Zero
interface Zero
  module procedure RS_SparseMatrixD_Zero
  module procedure RS_SparseMatrixZ_Zero
end interface Zero

public :: Print
interface Print
  module procedure RS_SparseMatrixD_Print
  module procedure RS_SparseMatrixZ_Print
end interface Print

public :: Print_simple
interface Print_simple
  module procedure RS_SparseMatrixD_Print_simple
!  module procedure RS_SparseMatrixZ_Print
end interface Print_simple

public :: copy
interface copy
  module procedure copy_dd_dsp
  module procedure copy_zd_zsp
  module procedure copy_dsp_dd
  module procedure copy_zsp_zd
end interface

public :: matrix_product_sub
interface matrix_product_sub
  module procedure matrix_product_sub_dsp_dsp_dsp, matrix_product_sub_dden_dden_dsp
  module procedure matrix_product_sub_zsp_zsp_zsp
  module procedure matrix_product_sub_dden_dden_zsp, matrix_product_sub_zden_dden_zsp, matrix_product_sub_zden_zden_zsp
  module procedure matrix_product_sub_zden_zden_dsp
end interface matrix_product_sub

public :: add_block
interface add_block
  module procedure RS_SparseMatrixD_add_block
  module procedure RS_SparseMatrixZ_add_block
end interface

public :: partial_TraceMult
interface partial_TraceMult
  module procedure RS_SparseMatrix_partial_TraceMult_dden_dsp
  module procedure RS_SparseMatrix_partial_TraceMult_dden_zsp
  module procedure RS_SparseMatrix_partial_TraceMult_zden_dsp
  module procedure RS_SparseMatrix_partial_TraceMult_zden_zsp
  module procedure RS_SparseMatrix_partial_TraceMult_dsp_zden
  module procedure RS_SparseMatrix_partial_TraceMult_dsp_dden
  module procedure RS_SparseMatrix_partial_TraceMult_zsp_zden
end interface partial_TraceMult

public :: TraceMult
interface TraceMult
  module procedure RS_SparseMatrix_TraceMult_dden_dsp
  module procedure RS_SparseMatrix_TraceMult_dden_zsp
  module procedure RS_SparseMatrix_TraceMult_zden_dsp
  module procedure RS_SparseMatrix_TraceMult_zden_zsp
  module procedure RS_SparseMatrix_TraceMult_dsp_zden
  module procedure RS_SparseMatrix_TraceMult_dsp_dden
  module procedure RS_SparseMatrix_TraceMult_zsp_zden
end interface TraceMult

public :: multDiagRL
interface multDiagRL
  module procedure RS_SparseMatrixD_multDiagRL_d, RS_SparseMatrixZ_multDiagRL_d
end interface multDiagRL

public :: check_sparse
interface check_sparse
  module procedure check_sparse_layout
end interface check_sparse

interface get_dense_block
  module procedure get_dense_blockD, get_dense_blockZ
end interface get_dense_block

contains

subroutine RS_SparseMatrixL_Initialise_at(this, at, first_orb_of_atom, cutoff, mpi_obj)
  type(RS_SparseMatrixL), intent(inout) :: this
  type(Atoms), intent(in) :: at
  integer, intent(in) :: first_orb_of_atom(:)
  real(dp), intent(in), optional :: cutoff
  type(MPI_context), intent(in), optional :: mpi_obj

  integer :: i, ji, j
  integer :: cur_data_pos, cur_row_offset
  integer, allocatable :: neighbour_sorted(:)
  real(dp), allocatable :: neighbour_distance_sorted(:)
  integer :: n_uniq

  real(dp) :: my_cutoff

  if (present(mpi_obj)) then
    call System_abort("No support for parallel sparse matrix")
  end if

  my_cutoff = HUGE(1.0_dp)
  if (present(cutoff)) my_cutoff = cutoff

  call Wipe(this)

  this%N = at%N
  if (size(first_orb_of_atom) /= at%N+1) call system_abort("RS_SparseMatrixL_Initialise_at size(first_orb_of_atom) = " // size(first_orb_of_atom) // "  != at%N+1 = " // (at%N+1))
  this%N_dense_rows = first_orb_of_atom(at%N+1)-1

  allocate(this%block_size(at%N))
  call ALLOC_TRACE("RS_L_init block_size", size(this%block_size)*INTEGER_SIZE)
  allocate(this%dense_row_of_row(at%N))
  call ALLOC_TRACE("RS_L_init dense_row_of_row", size(this%dense_row_of_row)*INTEGER_SIZE)
  this%block_size(:) = first_orb_of_atom(2:at%N+1)-first_orb_of_atom(1:at%N)
  this%dense_row_of_row(:) = first_orb_of_atom(1:at%N)

  allocate(this%row_indices(at%N+1))
  call ALLOC_TRACE("RS_L_init row_indices", size(this%row_indices)*INTEGER_SIZE)

  this%n_blocks = 0
  this%data_size = 0

  this%row_indices = 0
  do i=1, at%N

    allocate(neighbour_sorted(n_neighbours(at, i)))
    allocate(neighbour_distance_sorted(n_neighbours(at, i)))
    do ji=1, n_neighbours(at, i)
      neighbour_sorted(ji) = neighbour(at, i, ji, neighbour_distance_sorted(ji))
    end do
    call sort_array(neighbour_sorted, neighbour_distance_sorted)

    n_uniq = uniq_minval(neighbour_sorted, neighbour_distance_sorted)

    if(n_uniq == 0) call system_abort("atoms with no neighbours, not even itself!?!?!?")

    do ji=1, n_uniq
      j = neighbour_sorted(ji)
      if (ji > 1) then
	if (j == neighbour_sorted(ji-1)) cycle
      end if

      ! insert data for regular element
      if (neighbour_distance_sorted(ji) <= my_cutoff) then
	this%row_indices(i+1) = this%row_indices(i+1) + 1
	this%n_blocks = this%n_blocks + 1
	this%data_size = this%data_size + this%block_size(i)*this%block_size(j)
      end if

    end do

    deallocate(neighbour_sorted)
    deallocate(neighbour_distance_sorted)

  end do

  this%row_indices(1) = 1
  do i=2, at%N+1
    this%row_indices(i) = this%row_indices(i) + this%row_indices(i-1)
  end do

  if (this%n_blocks == 0) then
    allocate(this%col(1))
    allocate(this%data_ptrs(1))
  else
    allocate(this%col(this%n_blocks))
    allocate(this%data_ptrs(this%n_blocks))
  end if
  call ALLOC_TRACE("RS_L_init col", size(this%col)*INTEGER_SIZE)
  call ALLOC_TRACE("RS_L_init data_ptrs", size(this%data_ptrs)*INTEGER_SIZE)

  cur_data_pos = 1
  do i=1, at%N

    allocate(neighbour_sorted(n_neighbours(at, i)))
    allocate(neighbour_distance_sorted(n_neighbours(at, i)))
    do ji=1, n_neighbours(at, i)
      neighbour_sorted(ji) = neighbour(at, i, ji, neighbour_distance_sorted(ji))
    end do
    call sort_array(neighbour_sorted, neighbour_distance_sorted)

    n_uniq = uniq_minval(neighbour_sorted, neighbour_distance_sorted)

    if(n_uniq == 0) call system_abort("atoms with no neighbours, not even itself!?!?!?")

    cur_row_offset = 0

    do ji=1, n_uniq
      j = neighbour_sorted(ji)
      if (ji > 1) then
	if (j == neighbour_sorted(ji-1)) cycle
      end if

      ! insert data for regular element
      if (neighbour_distance_sorted(ji) <= my_cutoff) then
	this%col(this%row_indices(i)+cur_row_offset) = j
	this%data_ptrs(this%row_indices(i)+cur_row_offset) = cur_data_pos
	cur_row_offset = cur_row_offset + 1
	cur_data_pos = cur_data_pos + this%block_size(i)*this%block_size(j)
      end if

    end do

    deallocate(neighbour_sorted)
    deallocate(neighbour_distance_sorted)

  end do
end subroutine RS_SparseMatrixL_Initialise_at

function uniq_minval(i_data, r_data)
  integer, intent(inout) :: i_data(:)
  real(dp), intent(inout) :: r_data(:)
  integer :: uniq_minval

  real(dp) :: min_r_data
  integer :: data_len, i, ii, imin, imax

  data_len = size(i_data)
  if (data_len == 0) then
    uniq_minval=0
    return
  end if

  ii = 1
  do i=1, data_len
    if (i == 1) then
      imin = i
      cycle
    end if
    if (i_data(i) /= i_data(i-1)) then
      imax = i-1
      min_r_data = minval(r_data(imin:imax))
      i_data(ii) = i_data(imin)
      r_data(ii) = min_r_data
      ii = ii + 1
      imin = i
    end if
  end do
  min_r_data = minval(r_data(imin:data_len))
  i_data(ii) = i_data(imin)
  r_data(ii) = min_r_data

  uniq_minval = ii
end function uniq_minval

subroutine RS_SparseMatrixD_Initialise_at(this, at, first_orb_of_atom, cutoff, mpi_obj)
  type(RS_SparseMatrixD), intent(inout) :: this
  type(Atoms), intent(in) :: at
  integer, intent(in) :: first_orb_of_atom(:)
  real(dp), intent(in), optional :: cutoff
  type(MPI_context), intent(in), optional :: mpi_obj

  call Wipe(this)

  call Initialise(this%l, at, first_orb_of_atom, cutoff, mpi_obj)

  allocate(this%data(this%l%data_size))
  call ALLOC_TRACE("RS_D_init data", size(this%data)*REAL_SIZE)
end subroutine RS_SparseMatrixD_Initialise_at

subroutine RS_SparseMatrixZ_Initialise_at(this, at, first_orb_of_atom, cutoff, mpi_obj)
  type(RS_SparseMatrixZ), intent(inout) :: this
  type(Atoms), intent(in) :: at
  integer, intent(in) :: first_orb_of_atom(:)
  real(dp), intent(in), optional :: cutoff
  type(MPI_context), intent(in), optional :: mpi_obj

  call Wipe(this)

  call Initialise(this%l, at, first_orb_of_atom, cutoff, mpi_obj)

  allocate(this%data(this%l%data_size))
  call ALLOC_TRACE("RS_Z_init data", size(this%data)*COMPLEX_SIZE)
end subroutine RS_SparseMatrixZ_Initialise_at

subroutine RS_SparseMatrixL_Initialise_prod(this, al, bl, mpi_obj)
  type(RS_SparseMatrixL), intent(inout) :: this
  type(RS_SparseMatrixL), intent(in) :: al, bl
  type(MPI_context), intent(in), optional :: mpi_obj

  integer is, kk, ks, jj, js
  integer cur_entry, p, pd
  logical, allocatable :: col_used(:)
  integer, allocatable :: col_list(:)

  if (present(mpi_obj)) then
    call System_abort("No support for mpi_obj in sparse matrix")
  end if

  this%N = al%N
  this%N_dense_rows = al%N_dense_rows
  allocate(this%block_size(this%N))
  call ALLOC_TRACE("RS_L_init_product block_size", size(this%block_size)*INTEGER_SIZE)
  allocate(this%dense_row_of_row(this%N))
  call ALLOC_TRACE("RS_L_init_product dense_row_of_row", size(this%dense_row_of_row)*INTEGER_SIZE)
  this%block_size = al%block_size
  this%dense_row_of_row = al%dense_row_of_row

  allocate(this%row_indices(this%N+1))
  call ALLOC_TRACE("RS_L_init_product row_indices", size(this%row_indices)*INTEGER_SIZE)
  this%row_indices = 0
  this%n_blocks = 0
  this%data_size = 0

  allocate(col_used(al%N))
  allocate(col_list(al%N))
  col_used = .false.

  do is=1, al%N
    do kk=al%row_indices(is), al%row_indices(is+1)-1
      ks = al%col(kk)
      cur_entry = 0
      do jj=bl%row_indices(ks), bl%row_indices(ks+1)-1
	js = bl%col(jj)
	if (.not.col_used(js)) then
	  cur_entry = cur_entry + 1
	  col_list(cur_entry) = js
	  col_used(js) = .true.
	  this%n_blocks = this%n_blocks + 1
	  this%data_size = this%data_size + al%block_size(is)*bl%block_size(js)
	end if
      end do
    end do
    do jj=1, cur_entry
      col_used(col_list(jj)) = .false.
    end do
  end do

  allocate(this%col(this%n_blocks))
  allocate(this%data_ptrs(this%n_blocks))
  call ALLOC_TRACE("RS_L_init_product col", size(this%col)*INTEGER_SIZE)
  call ALLOC_TRACE("RS_L_init_product data_ptrs", size(this%data_ptrs)*INTEGER_SIZE)

  p = 1
  pd = 1
  do is=1, al%N
    this%row_indices(is) = p
    do kk=al%row_indices(is), al%row_indices(is+1)-1
      ks = al%col(kk)
      cur_entry = 0
      do jj=bl%row_indices(ks), bl%row_indices(ks+1)-1
	js = bl%col(jj)
	if (.not.col_used(js)) then
	  cur_entry = cur_entry + 1
	  col_list(cur_entry) = js
	  col_used(js) = .true.

	  this%col(p) = js
	  this%data_ptrs(p) = pd
	  p = p + 1
	  pd = pd + al%block_size(is)*bl%block_size(js)
	end if
      end do
    end do
    do jj=1, cur_entry
      col_used(col_list(jj)) = .false.
    end do
  end do
  this%row_indices(al%N+1) = p

  deallocate(col_used)
  deallocate(col_list)
end subroutine RS_SparseMatrixL_Initialise_prod

subroutine RS_SparseMatrixD_Initialise_prod(this, al, bl, mpi_obj)
  type(RS_SparseMatrixD), intent(inout) :: this
  type(RS_SparseMatrixL), intent(in) :: al, bl
  type(MPI_context), intent(in), optional :: mpi_obj

  call Wipe(this)

  call Initialise(this%l, al, bl)

  allocate(this%data(this%l%data_size))
  call ALLOC_TRACE("RS_D_init_product data", size(this%data)*REAL_SIZE)

end subroutine RS_SparseMatrixD_Initialise_prod

subroutine RS_SparseMatrixZ_Initialise_prod(this, al, bl, mpi_obj)
  type(RS_SparseMatrixZ), intent(inout) :: this
  type(RS_SparseMatrixL), intent(in) :: al, bl
  type(MPI_context), intent(in), optional :: mpi_obj

  call Wipe(this)

  call Initialise(this%l, al, bl)

  allocate(this%data(this%l%data_size))
  call ALLOC_TRACE("RS_Z_init_product data", size(this%data)*COMPLEX_SIZE)

end subroutine RS_SparseMatrixZ_Initialise_prod

subroutine RS_SparseMatrixL_Finalise(this)
  type(RS_SparseMatrixL), intent(inout) :: this

  call Wipe(this)
end subroutine RS_SparseMatrixL_Finalise

subroutine RS_SparseMatrixD_Finalise(this)
  type(RS_SparseMatrixD), intent(inout) :: this

  call Wipe(this)
end subroutine RS_SparseMatrixD_Finalise

subroutine RS_SparseMatrixZ_Finalise(this)
  type(RS_SparseMatrixZ), intent(inout) :: this

  call Wipe(this)
end subroutine RS_SparseMatrixZ_Finalise

subroutine RS_SparseMatrixL_Wipe(this)
  type(RS_SparseMatrixL), intent(inout) :: this

  call DEALLOC_TRACE("RS_L_wipe block_size", size(this%block_size)*INTEGER_SIZE)
  if (allocated(this%block_size)) deallocate(this%block_size)
  call DEALLOC_TRACE("RS_L_wipe dense_row_of_row", size(this%dense_row_of_row)*INTEGER_SIZE)
  if (allocated(this%dense_row_of_row)) deallocate(this%dense_row_of_row)
  call DEALLOC_TRACE("RS_L_wipe row_indices", size(this%row_indices)*INTEGER_SIZE)
  if (allocated(this%row_indices)) deallocate(this%row_indices)
  call DEALLOC_TRACE("RS_L_wipe data_ptrs", size(this%data_ptrs)*INTEGER_SIZE)
  if (allocated(this%data_ptrs)) deallocate(this%data_ptrs)
  call DEALLOC_TRACE("RS_L_wipe col", size(this%col)*INTEGER_SIZE)
  if (allocated(this%col)) deallocate(this%col)

  this%N = 0
  this%N_dense_rows = 0
  this%n_blocks = 0
  this%data_size = 0
end subroutine RS_SparseMatrixL_Wipe

subroutine RS_SparseMatrixD_Wipe(this)
  type(RS_SparseMatrixD), intent(inout) :: this

  call Wipe(this%l)
  call DEALLOC_TRACE("RS_D_wipe data", size(this%data)*REAL_SIZE)
  if (allocated(this%data)) deallocate(this%data)

end subroutine RS_SparseMatrixD_Wipe

subroutine RS_SparseMatrixZ_Wipe(this)
  type(RS_SparseMatrixZ), intent(inout) :: this

  call Wipe(this%l)
  call DEALLOC_TRACE("RS_Z_wipe data", size(this%data)*COMPLEX_SIZE)
  if (allocated(this%data)) deallocate(this%data)

end subroutine RS_SparseMatrixZ_Wipe

subroutine RS_SparseMatrixD_Zero(this, d_mask, od_mask)
  type(RS_SparseMatrixD),    intent(inout)           :: this
  logical, intent(in), optional :: d_mask(:), od_mask(:)

  integer i, ji, j, block_nr, block_nc

  if (present(d_mask) .or. present(od_mask)) then
    do i=1, this%l%N
      block_nr = this%l%block_size(i)
      do ji=this%l%row_indices(i), this%l%row_indices(i+1)-1
	j = this%l%col(ji)
	block_nc = this%l%block_size(j)
	if (present(d_mask) .and. i == j) then
	  if (d_mask(i)) this%data(this%l%data_ptrs(ji):this%l%data_ptrs(ji)+block_nr*block_nc-1) = 0.0_dp
	end if
	if (present(od_mask) .and. i /= j) then
	  if (od_mask(i) .or. od_mask(j)) this%data(this%l%data_ptrs(ji):this%l%data_ptrs(ji)+block_nr*block_nc-1) = 0.0_dp
	end if
      end do
    end do
  else
    this%data = 0.0_dp
  end if

end subroutine RS_SparseMatrixD_Zero

subroutine RS_SparseMatrixZ_Zero(this, d_mask, od_mask)
  type(RS_SparseMatrixZ),    intent(inout)           :: this
  logical, intent(in), optional :: d_mask(:), od_mask(:)

  integer i, ji, j, block_nr, block_nc

  if (present(d_mask) .or. present(od_mask)) then
    do i=1, this%l%N
      block_nr = this%l%block_size(i)
      do ji=this%l%row_indices(i), this%l%row_indices(i+1)-1
	j = this%l%col(ji)
	block_nc = this%l%block_size(j)
	if (present(d_mask) .and. i == j) then
	  if (d_mask(i)) this%data(this%l%data_ptrs(ji):this%l%data_ptrs(ji)+block_nr*block_nc-1) = 0.0_dp
	end if
	if (present(od_mask) .and. i /= j) then
	  if (od_mask(i) .or. od_mask(j)) this%data(this%l%data_ptrs(ji):this%l%data_ptrs(ji)+block_nr*block_nc-1) = 0.0_dp
	end if
      end do
    end do
  else
    this%data = 0.0_dp
  end if

end subroutine RS_SparseMatrixZ_Zero

subroutine RS_SparseMatrixD_Print(this,file)
  type(RS_SparseMatrixD),    intent(in)           :: this
  type(Inoutput), intent(inout),optional:: file

  integer :: i, j
  integer :: block_nr, block_nc

  call Print ('RS_SparseMatrixD : N ' // this%l%N, file=file)

  call Print ('RS_SparseMatrixD : n_blocks data_size ' // this%l%n_blocks // " " // this%l%data_size, file=file)

  do i=1, this%l%N
    call Print ('RS_SparseMatrixD : row, block_size dense_row ' // i // " " // this%l%block_size(i) // " " // &
					 this%l%dense_row_of_row(i), file=file)
    do j=this%l%row_indices(i), this%l%row_indices(i+1)-1
      call Print ("RS_SparseMatrixD entry j "// j // " col " // this%l%col(j) //  " data_ptr " // this%l%data_ptrs(j), &
	file=file)
      block_nr = this%l%block_size(i)
      block_nc = this%l%block_size(this%l%col(j))
      call Print(reshape( &
	this%data(this%l%data_ptrs(j):this%l%data_ptrs(j)+block_nr*block_nc-1),(/block_nr,block_nc/) ), &
	file=file)
    end do
  end do

end subroutine RS_SparseMatrixD_Print

subroutine RS_SparseMatrixZ_Print(this,file)
  type(RS_SparseMatrixZ),    intent(in)           :: this
  type(Inoutput), intent(inout),optional:: file

  integer :: i, j
  integer :: block_nr, block_nc

  call Print('RS_SparseMatrixZ : ', file=file)

  call Print ('RS_SparseMatrixZ : N ' //  this%l%N, file=file)

  call Print ('RS_SparseMatrixZ : n_blocks data_size ' // this%l%n_blocks // this%l%data_size, file=file)

  do i=1, this%l%N
    call Print ('RS_SparseMatrixZ : row, block_size dense_row ' // i // " " //  this%l%block_size(i) //  " " // &
					 this%l%dense_row_of_row(i), file=file)
    do j=this%l%row_indices(i), this%l%row_indices(i+1)-1
      call Print ("RS_SparseMatrixZ entry j " // j // " col " // this%l%col(j) //  " data_ptr " // this%l%data_ptrs(j), file=file)
      block_nr = this%l%block_size(i)
      block_nc = this%l%block_size(this%l%col(j))
      call Print(reshape( &
	this%data(this%l%data_ptrs(j):this%l%data_ptrs(j)+block_nr*block_nc-1),(/block_nr,block_nc/) ), file=file)
    end do
  end do

end subroutine RS_SparseMatrixZ_Print

subroutine RS_SparseMatrixD_Print_simple(this,file)
  type(RS_SparseMatrixD),    intent(in)           :: this
  type(Inoutput), intent(inout),optional:: file

  integer :: i, tj, j, ii, jj
  integer :: block_nr, block_nc
  integer :: n_entries

  n_entries = 0
  do i=1, this%l%N
    block_nr = this%l%block_size(i)
    do j=this%l%row_indices(i), this%l%row_indices(i+1)-1
      block_nc = this%l%block_size(this%l%col(j))
      n_entries = n_entries + block_nr*block_nc
    end do
  end do

  call print(this%l%N // " " // n_entries, file=file)

  do i=1, this%l%N
    block_nr = this%l%block_size(i)
    do tj=this%l%row_indices(i), this%l%row_indices(i+1)-1
      j = this%l%col(tj)
      do ii=1, block_nr
      do jj=1, block_nc
	call print((this%l%dense_row_of_row(i)+ii-1) // " " // (this%l%dense_row_of_row(j)+jj-1) // " " // &
	  this%data(this%l%data_ptrs(tj)+ii-1+(jj-1)*block_nr), file=file)
      end do
      end do
    end do
  end do
end subroutine RS_SparseMatrixD_Print_simple

subroutine copy_dsp_dd(this, dense)
  type(RS_SparseMatrixD), intent(in) :: this
  real(dp), intent(inout), dimension(:,:) :: dense

  integer :: i, tj, j, ii, jj
  integer :: block_nr, block_nc
  integer :: n_entries

  !if (any(shape(dense) /= (/this%l%N, this%l%N/))) call system_abort("size of dense matrix wrong: expecting "//this%l%N)

  n_entries = 0
  do i=1, this%l%N
    block_nr = this%l%block_size(i)
    do j=this%l%row_indices(i), this%l%row_indices(i+1)-1
      block_nc = this%l%block_size(this%l%col(j))
      n_entries = n_entries + block_nr*block_nc
    end do
  end do

  do i=1, this%l%N
    block_nr = this%l%block_size(i)
    do tj=this%l%row_indices(i), this%l%row_indices(i+1)-1
      j = this%l%col(tj)
      do ii=1, block_nr
        do jj=1, block_nc
  	       dense(this%l%dense_row_of_row(i)+ii-1, &
                 this%l%dense_row_of_row(j)+jj-1) = &
  	             this%data(this%l%data_ptrs(tj)+ii-1+(jj-1)*block_nr)
        end do
      end do
    end do
  end do
end subroutine copy_dsp_dd

subroutine copy_zsp_zd(this, dense)
  type(RS_SparseMatrixZ), intent(in) :: this
  complex(dp), intent(inout), dimension(:,:) :: dense

  integer :: i, tj, j, ii, jj
  integer :: block_nr, block_nc
  integer :: n_entries

  if (any(shape(dense) /= (/this%l%N, this%l%N/))) call system_abort("size of dense matrix wrong: expecting "//this%l%N)

  n_entries = 0
  do i=1, this%l%N
    block_nr = this%l%block_size(i)
    do j=this%l%row_indices(i), this%l%row_indices(i+1)-1
      block_nc = this%l%block_size(this%l%col(j))
      n_entries = n_entries + block_nr*block_nc
    end do
  end do

  do i=1, this%l%N
    block_nr = this%l%block_size(i)
    do tj=this%l%row_indices(i), this%l%row_indices(i+1)-1
      j = this%l%col(tj)
      do ii=1, block_nr
        do jj=1, block_nc
  	       dense(this%l%dense_row_of_row(i)+ii-1, &
                 this%l%dense_row_of_row(j)+jj-1) = &
  	             this%data(this%l%data_ptrs(tj)+ii-1+(jj-1)*block_nr)
        end do
      end do
    end do
  end do
end subroutine copy_zsp_zd

subroutine copy_dd_dsp(a, b, d_mask, od_mask)
  type(MatrixD), intent(inout) :: a
  type(RS_SparseMatrixD), intent(in) :: b
  logical, optional :: d_mask(:), od_mask(:)

  integer is, jj, js, id, jd

  a%data = 0.0_dp
  do is=1, b%l%N
    id = b%l%dense_row_of_row(is)
    do jj=b%l%row_indices(is), b%l%row_indices(is+1)-1
      js = b%l%col(jj)
      if (present(d_mask)) then
	if ((is == js) .and. .not. d_mask(is)) cycle
      endif
      if (present(od_mask)) then
	if ((is /= js) .and. .not. od_mask(is) .and. .not. od_mask(js)) cycle
      endif
      jd = b%l%dense_row_of_row(js)
      a%data(id:id+b%l%block_size(is)-1,jd:jd+b%l%block_size(js)-1) = &
	reshape(b%data(b%l%data_ptrs(jj):b%l%data_ptrs(jj)+b%l%block_size(is)*b%l%block_size(js)-1), &
	        (/b%l%block_size(is), b%l%block_size(js)/) )
    end do
  end do
end subroutine copy_dd_dsp

subroutine copy_zd_zsp(a, b, d_mask, od_mask)
  type(MatrixZ), intent(inout) :: a
  type(RS_SparseMatrixZ), intent(in) :: b
  logical, optional :: d_mask(:), od_mask(:)

  integer is, jj, js, id, jd

  a%data = 0.0_dp
  do is=1, b%l%N
    id = b%l%dense_row_of_row(is)
    do jj=b%l%row_indices(is), b%l%row_indices(is+1)-1
      js = b%l%col(jj)
      if (present(d_mask)) then
	if ((is == js) .and. .not. d_mask(is)) cycle
      endif
      if (present(od_mask)) then
	if ((is /= js) .and. .not. od_mask(is) .and. .not. od_mask(js)) cycle
      endif
      jd = b%l%dense_row_of_row(js)
      a%data(id:id+b%l%block_size(is)-1,jd:jd+b%l%block_size(js)-1) = &
	reshape(b%data(b%l%data_ptrs(jj):b%l%data_ptrs(jj)+b%l%block_size(is)*b%l%block_size(js)-1), &
	        (/b%l%block_size(is), b%l%block_size(js)/) )
    end do
  end do
end subroutine copy_zd_zsp

subroutine matrix_product_sub_dsp_dsp_dsp(c,a,b, a_transpose, b_transpose)
  type(RS_SparseMatrixD), intent(inout) :: c
  type(RS_SparseMatrixD), intent(in) :: a
  type(RS_SparseMatrixD), intent(in) :: b
  logical, intent(in), optional :: a_transpose, b_transpose

  integer is, kk, ks, jj_b, jj_c
  integer block_ni, block_nj, block_nk
  integer max_block_size
  real(dp), allocatable :: a_block(:,:)

  c%data = 0.0_dp

  if (present(a_transpose)) then
    if (a_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if
  if (present(b_transpose)) then
    if (b_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if

  max_block_size = maxval(a%l%block_size)
  allocate(a_block(max_block_size,max_block_size))

  do is=1, a%l%N
    block_ni = a%l%block_size(is)
    do kk=a%l%row_indices(is), a%l%row_indices(is+1)-1
      ks = a%l%col(kk)
      block_nk = a%l%block_size(ks)
      a_block(1:block_ni,1:block_nk) = reshape(a%data(a%l%data_ptrs(kk):a%l%data_ptrs(kk)+block_ni*block_nk-1), &
					       (/block_ni,block_nk/))

      jj_c = c%l%row_indices(is)
      jj_b = b%l%row_indices(ks)
      do while (jj_c < c%l%row_indices(is+1) .and. jj_b < b%l%row_indices(ks+1))
	if (b%l%col(jj_b) == c%l%col(jj_c)) then
	  block_nj = b%l%block_size(b%l%col(jj_b))
	  c%data(c%l%data_ptrs(jj_c):c%l%data_ptrs(jj_c)+block_ni*block_nj-1) = &
	    c%data(c%l%data_ptrs(jj_c):c%l%data_ptrs(jj_c)+block_ni*block_nj-1) + reshape( &
	      matmul( a_block(1:block_ni,1:block_nk), &
		      reshape(b%data(b%l%data_ptrs(jj_b):b%l%data_ptrs(jj_b)+block_nk*block_nj-1), &
			      (/block_nk,block_nj/)) ), &
	    (/block_ni*block_nj/) )
	  jj_b = jj_b + 1
	  jj_c = jj_c + 1
	else if (b%l%col(jj_b) < c%l%col(jj_c)) then
	  jj_b = jj_b + 1
	else
	  jj_c = jj_c + 1
	end if
      end do

    end do
  end do

  deallocate(a_block)

end subroutine matrix_product_sub_dsp_dsp_dsp

subroutine matrix_product_sub_zsp_zsp_zsp(c,a,b, a_transpose, a_conjugate, b_transpose, b_conjugate)
  type(RS_SparseMatrixZ), intent(inout) :: c
  type(RS_SparseMatrixZ), intent(in) :: a
  type(RS_SparseMatrixZ), intent(in) :: b
  logical, intent(in), optional :: a_transpose, a_conjugate, b_transpose, b_conjugate

  integer is, kk, ks, jj_b, jj_c
  integer block_ni, block_nj, block_nk
  integer max_block_size
  complex(dp), allocatable :: a_block(:,:)

  c%data = 0.0_dp

  if (present(a_transpose)) then
    if (a_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if
  if (present(b_transpose)) then
    if (b_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if
  if (present(a_conjugate)) then
    if (a_conjugate) call system_abort("Can't do sparse = sparse * sparse with conjugate matrices")
  end if
  if (present(b_conjugate)) then
    if (b_conjugate) call system_abort("Can't do sparse = sparse * sparse with conjugate matrices")
  end if

  max_block_size = maxval(a%l%block_size)
  allocate(a_block(max_block_size,max_block_size))

  do is=1, a%l%N
    block_ni = a%l%block_size(is)
    do kk=a%l%row_indices(is), a%l%row_indices(is+1)-1
      ks = a%l%col(kk)
      block_nk = a%l%block_size(ks)
      a_block(1:block_ni,1:block_nk) = reshape(a%data(a%l%data_ptrs(kk):a%l%data_ptrs(kk)+block_ni*block_nk-1), &
					       (/block_ni,block_nk/))

      jj_c = c%l%row_indices(is)
      jj_b = b%l%row_indices(ks)
      do while (jj_c < c%l%row_indices(is+1) .and. jj_b < b%l%row_indices(ks+1))
	if (b%l%col(jj_b) == c%l%col(jj_c)) then
	  block_nj = b%l%block_size(b%l%col(jj_b))
	  c%data(c%l%data_ptrs(jj_c):c%l%data_ptrs(jj_c)+block_ni*block_nj-1) = &
	    c%data(c%l%data_ptrs(jj_c):c%l%data_ptrs(jj_c)+block_ni*block_nj-1) + reshape( &
	      matmul( a_block(1:block_ni,1:block_nk), &
		      reshape(b%data(b%l%data_ptrs(jj_b):b%l%data_ptrs(jj_b)+block_nk*block_nj-1), &
			      (/block_nk,block_nj/)) ), &
	    (/block_ni*block_nj/) )
	  jj_b = jj_b + 1
	  jj_c = jj_c + 1
	else if (b%l%col(jj_b) < c%l%col(jj_c)) then
	  jj_b = jj_b + 1
	else
	  jj_c = jj_c + 1
	end if
      end do

    end do
  end do

  deallocate(a_block)

end subroutine matrix_product_sub_zsp_zsp_zsp

subroutine matrix_product_sub_dden_dden_dsp(c,a,b,a_transpose, b_transpose, diag_mask, offdiag_mask)
  type(MatrixD), intent(inout) :: c
  type(MatrixD), intent(in) :: a
  type(RS_SparseMatrixD), intent(in) :: b
  logical, intent(in), optional :: a_transpose, b_transpose
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)

  integer id, is, kd, ks, jj, jd, js
  integer block_ni, block_nj, block_nk

  c%data = 0.0_dp

  if (present(a_transpose)) then
    if (a_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if
  if (present(b_transpose)) then
    if (b_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if

  id = 1
  is = 1
  do while (id <= a%N)
    block_ni = b%l%block_size(is)
    kd = 1
    ks = 1
    do while (kd <= a%M)
      block_nk = b%l%block_size(ks)
      do jj=b%l%row_indices(ks), b%l%row_indices(ks+1)-1
	js = b%l%col(jj)
	jd = b%l%dense_row_of_row(js)
	block_nj = b%l%block_size(js)

	if (present(diag_mask)) then
	  if (js == ks .and. .not. diag_mask(ks)) cycle
	end if
	if (present(offdiag_mask)) then
	  if (js /= ks .and. .not. offdiag_mask(ks) .and. .not. offdiag_mask(js)) cycle
	end if

	c%data(id:id+block_ni-1,jd:jd+block_nj-1) = &
	  c%data(id:id+block_ni-1,jd:jd+block_nj-1) + &
	  matmul( a%data(id:id+block_ni-1,kd:kd+block_nk-1), &
		  reshape(b%data(b%l%data_ptrs(jj):b%l%data_ptrs(jj)+block_nk*block_nj-1),  &
			  (/block_nk,block_nj/)) )
      end do
      kd = kd + b%l%block_size(ks)
      ks = ks + 1
    end do
    id = id + b%l%block_size(is)
    is = is + 1
  end do
end subroutine

subroutine matrix_product_sub_dden_dden_zsp(c,a,b,a_transpose, b_transpose, diag_mask, offdiag_mask)
  type(MatrixD), intent(inout) :: c
  type(MatrixD), intent(in) :: a
  type(RS_SparseMatrixZ), intent(in) :: b
  logical, intent(in), optional :: a_transpose, b_transpose
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)

  integer id, is, kd, ks, jj, jd, js
  integer block_ni, block_nj, block_nk

  c%data = 0.0_dp

  if (present(a_transpose)) then
    if (a_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if
  if (present(b_transpose)) then
    if (b_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if

  id = 1
  is = 1
  do while (id <= a%N)
    block_ni = b%l%block_size(is)
    kd = 1
    ks = 1
    do while (kd <= a%M)
      block_nk = b%l%block_size(ks)
      do jj=b%l%row_indices(ks), b%l%row_indices(ks+1)-1
	js = b%l%col(jj)
	jd = b%l%dense_row_of_row(js)
	block_nj = b%l%block_size(js)

        if (present(diag_mask)) then
          if (js == ks .and. .not. diag_mask(ks)) cycle
        end if
        if (present(offdiag_mask)) then
          if (js /= ks .and. .not. offdiag_mask(ks) .and. .not. offdiag_mask(js)) cycle
        end if

	c%data(id:id+block_ni-1,jd:jd+block_nj-1) = &
	  c%data(id:id+block_ni-1,jd:jd+block_nj-1) + &
	  matmul( a%data(id:id+block_ni-1,kd:kd+block_nk-1), &
		  reshape(real(b%data(b%l%data_ptrs(jj):b%l%data_ptrs(jj)+block_nk*block_nj-1)),  &
			  (/block_nk,block_nj/)) )
      end do
      kd = kd + b%l%block_size(ks)
      ks = ks + 1
    end do
    id = id + b%l%block_size(is)
    is = is + 1
  end do
end subroutine

subroutine matrix_product_sub_zden_dden_zsp(c,a,b,a_transpose, b_transpose, b_conjugate, diag_mask, offdiag_mask)
  type(MatrixZ), intent(inout) :: c
  type(MatrixD), intent(in) :: a
  type(RS_SparseMatrixZ), intent(in) :: b
  logical, intent(in), optional :: a_transpose, b_transpose, b_conjugate
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)

  integer id, is, kd, ks, jj, jd, js
  integer block_ni, block_nj, block_nk

  c%data = 0.0_dp

  if (present(a_transpose)) then
    if (a_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if
  if (present(b_transpose)) then
    if (b_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if
  if (present(b_conjugate)) then
    if (b_conjugate) call system_abort("Can't do sparse = sparse * sparse with conjugate matrices")
  end if

  id = 1
  is = 1
  do while (id <= a%N)
    block_ni = b%l%block_size(is)
    kd = 1
    ks = 1
    do while (kd <= a%M)
      block_nk = b%l%block_size(ks)
      do jj=b%l%row_indices(ks), b%l%row_indices(ks+1)-1
	js = b%l%col(jj)
	jd = b%l%dense_row_of_row(js)
	block_nj = b%l%block_size(js)

        if (present(diag_mask)) then
          if (js == ks .and. .not. diag_mask(ks)) cycle
        end if
        if (present(offdiag_mask)) then
          if (js /= ks .and. .not. offdiag_mask(ks) .and. .not. offdiag_mask(js)) cycle
        end if

	c%data(id:id+block_ni-1,jd:jd+block_nj-1) = &
	  c%data(id:id+block_ni-1,jd:jd+block_nj-1) + &
	  matmul( a%data(id:id+block_ni-1,kd:kd+block_nk-1), &
		  reshape(b%data(b%l%data_ptrs(jj):b%l%data_ptrs(jj)+block_nk*block_nj-1),  &
			  (/block_nk,block_nj/)) )
      end do
      kd = kd + b%l%block_size(ks)
      ks = ks + 1
    end do
    id = id + b%l%block_size(is)
    is = is + 1
  end do
end subroutine

subroutine matrix_product_sub_zden_zden_zsp(c,a,b,a_transpose, a_conjugate, b_transpose, b_conjugate, diag_mask, offdiag_mask)
  type(MatrixZ), intent(inout) :: c
  type(MatrixZ), intent(in) :: a
  type(RS_SparseMatrixZ), intent(in) :: b
  logical, intent(in), optional :: a_transpose, a_conjugate, b_transpose, b_conjugate
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)

  integer id, is, kd, ks, jj, jd, js
  integer block_ni, block_nj, block_nk

  c%data = 0.0_dp

  if (present(a_transpose)) then
    if (a_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if
  if (present(b_transpose)) then
    if (b_transpose) call system_abort("Can't do sparse = sparse * sparse with transposed matrices")
  end if
  if (present(a_conjugate)) then
    if (a_conjugate) call system_abort("Can't do sparse = sparse * sparse with conjugate matrices")
  end if
  if (present(b_conjugate)) then
    if (b_conjugate) call system_abort("Can't do sparse = sparse * sparse with conjugate matrices")
  end if

  id = 1
  is = 1
  do while (id <= a%N)
    block_ni = b%l%block_size(is)
    kd = 1
    ks = 1
    do while (kd <= a%M)
      block_nk = b%l%block_size(ks)
      do jj=b%l%row_indices(ks), b%l%row_indices(ks+1)-1
	js = b%l%col(jj)
	jd = b%l%dense_row_of_row(js)
	block_nj = b%l%block_size(js)

        if (present(diag_mask)) then
          if (js == ks .and. .not. diag_mask(ks)) cycle
        end if
        if (present(offdiag_mask)) then
          if (js /= ks .and. .not. offdiag_mask(ks) .and. .not. offdiag_mask(js)) cycle
        end if

	c%data(id:id+block_ni-1,jd:jd+block_nj-1) = &
	  c%data(id:id+block_ni-1,jd:jd+block_nj-1) + &
	  matmul( a%data(id:id+block_ni-1,kd:kd+block_nk-1), &
		  reshape(b%data(b%l%data_ptrs(jj):b%l%data_ptrs(jj)+block_nk*block_nj-1),  &
			  (/block_nk,block_nj/)) )
      end do
      kd = kd + b%l%block_size(ks)
      ks = ks + 1
    end do
    id = id + b%l%block_size(is)
    is = is + 1
  end do
end subroutine

subroutine matrix_product_sub_zden_zden_dsp(c,a,b,a_transpose, a_conjugate, b_transpose, diag_mask, offdiag_mask)
  type(MatrixZ), intent(inout) :: c
  type(MatrixZ), intent(in) :: a
  type(RS_SparseMatrixD), intent(in) :: b
  logical, intent(in), optional :: a_transpose, a_conjugate, b_transpose
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)

  integer id, is, kd, ks, jj, jd, js
  integer block_ni, block_nj, block_nk

  if (present(a_transpose)) then
    if (a_transpose) call system_abort("Can't do dense = dense * sparse with transposed matrices")
  end if
  if (present(b_transpose)) then
    if (b_transpose) call system_abort("Can't do dense = dense * sparse with transposed matrices")
  end if
  if (present(a_conjugate)) then
    if (a_conjugate) call system_abort("Can't do dense = dense * sparse with conjugate matrices")
  end if

  c%data = 0.0_dp

  id = 1
  is = 1
  do while (id <= a%N)
    block_ni = b%l%block_size(is)
    kd = 1
    ks = 1
    do while (kd <= a%M)
      block_nk = b%l%block_size(ks)
      do jj=b%l%row_indices(ks), b%l%row_indices(ks+1)-1
	js = b%l%col(jj)
	jd = b%l%dense_row_of_row(js)
	block_nj = b%l%block_size(js)

        if (present(diag_mask)) then
          if (js == ks .and. .not. diag_mask(ks)) cycle
        end if
        if (present(offdiag_mask)) then
          if (js /= ks .and. .not. offdiag_mask(ks) .and. .not. offdiag_mask(js)) cycle
        end if

	c%data(id:id+block_ni-1,jd:jd+block_nj-1) = &
	  c%data(id:id+block_ni-1,jd:jd+block_nj-1) + &
	  matmul( a%data(id:id+block_ni-1,kd:kd+block_nk-1), &
		  reshape(b%data(b%l%data_ptrs(jj):b%l%data_ptrs(jj)+block_nk*block_nj-1),  &
			  (/block_nk,block_nj/)) )
      end do
      kd = kd + b%l%block_size(ks)
      ks = ks + 1
    end do
    id = id + b%l%block_size(is)
    is = is + 1
  end do
end subroutine

subroutine RS_SparseMatrixD_add_block(this, block, block_nr, block_nc, at_row, at_col)
  type(RS_SparseMatrixD), intent(inout) :: this
  real(dp), intent(in) :: block(:,:)
  integer, intent(in) :: block_nr, block_nc
  integer, intent(in), optional :: at_row, at_col

  integer n_cols
  integer cur_col, col_offset

  if (.not. present(at_row) .or. .not. present(at_col)) then
    call System_abort("need at_row and at_col for add_block to RS_SparseMatrixD")
  end if

  if (at_row > this%l%N .or. at_col > this%l%N) then
      call system_abort("RS_SparseMatrixD_add_block tried to add block outside of matrix bounds "//at_row//","//at_col//" "//this%l%N)
  end if

  if (this%l%cur_row /= at_row) then
    this%l%cur_row = at_row
    this%l%cur_col_offset = 0
  end if

  cur_col = this%l%col(this%l%row_indices(this%l%cur_row)+this%l%cur_col_offset)
  if (cur_col == at_col) then
    call add_block_d(this, this%l%cur_row, this%l%cur_col_offset, block_nr, block_nc, block)
    return
  else if (cur_col > at_col) then
    do col_offset = this%l%cur_col_offset-1, 0, -1
      cur_col = this%l%col(this%l%row_indices(this%l%cur_row)+col_offset)
      if (cur_col == at_col) then
	this%l%cur_col_offset = col_offset
	call add_block_d(this, this%l%cur_row, this%l%cur_col_offset, block_nr, block_nc, block)
	return
      end if
    end do
  else
    n_cols = this%l%row_indices(this%l%cur_row+1)-this%l%row_indices(this%l%cur_row)
    do col_offset = this%l%cur_col_offset+1, n_cols-1
      cur_col = this%l%col(this%l%row_indices(this%l%cur_row)+col_offset)
      if (cur_col == at_col) then
	this%l%cur_col_offset = col_offset
	call add_block_d(this, this%l%cur_row, this%l%cur_col_offset, block_nr, block_nc, block)
	return
      end if
    end do
  end if

end subroutine RS_SparseMatrixD_add_block

subroutine RS_SparseMatrixZ_add_block(this, block, block_nr, block_nc, at_row, at_col)
  type(RS_SparseMatrixZ), intent(inout) :: this
  complex(dp), intent(in) :: block(:,:)
  integer, intent(in) :: block_nr, block_nc
  integer, intent(in) :: at_row, at_col

  integer n_cols
  integer cur_col, col_offset

  if (at_row > this%l%N .or. at_col > this%l%N) then
      call system_abort("RS_SparseMatrixZ_add_block_tried to add block outside of matrix bounds "//at_row//","//at_col//" "//this%l%N)
  end if

  if (this%l%cur_row /= at_row) then
    this%l%cur_row = at_row
    this%l%cur_col_offset = 0
  end if

  cur_col = this%l%col(this%l%row_indices(this%l%cur_row)+this%l%cur_col_offset)
  if (cur_col == at_col) then
    call add_block_z(this, this%l%cur_row, this%l%cur_col_offset, block_nr, block_nc, block)
    return
  else if (cur_col > at_col) then
    do col_offset = this%l%cur_col_offset-1, 0, -1
      cur_col = this%l%col(this%l%row_indices(this%l%cur_row)+col_offset)
      if (cur_col == at_col) then
	this%l%cur_col_offset = col_offset
	call add_block_z(this, this%l%cur_row, this%l%cur_col_offset, block_nr, block_nc, block)
	return
      end if
    end do
  else
    n_cols = this%l%row_indices(this%l%cur_row+1)-this%l%row_indices(this%l%cur_row)
    do col_offset = this%l%cur_col_offset+1, n_cols-1
      cur_col = this%l%col(this%l%row_indices(this%l%cur_row)+col_offset)
      if (cur_col == at_col) then
	this%l%cur_col_offset = col_offset
	call add_block_z(this, this%l%cur_row, this%l%cur_col_offset, block_nr, block_nc, block)
	return
      end if
    end do
  end if

end subroutine RS_SparseMatrixZ_add_block

subroutine add_block_d(this, cur_row, cur_col_offset, in_nr, in_nc, block)
  type(RS_SparseMatrixD), intent(inout) :: this
  integer, intent(in) :: cur_row, cur_col_offset
  integer, intent(in) :: in_nr, in_nc
  real(dp), intent(in) :: block(:,:)

  integer ptr_i, this_nr, this_nc

  ptr_i = this%l%row_indices(cur_row)+cur_col_offset
  this_nr = this%l%block_size(cur_row)
  this_nc = this%l%block_size(this%l%col(ptr_i))
  if (in_nr /= this_nr .or. in_nc /= this_nc) then
    call System_abort("add_block_d tried to add block of wrong shape in_nr,nc "//in_nr// " "//in_nc//" sp_nr,nc "//this_nr//" "//this_nc)
  end if
  this%data(this%l%data_ptrs(ptr_i):this%l%data_ptrs(ptr_i)+this_nr*this_nc-1) = &
    this%data(this%l%data_ptrs(ptr_i):this%l%data_ptrs(ptr_i)+this_nr*this_nc-1) + &
    reshape( block(1:in_nr,1:in_nc), (/this_nr*this_nc/) )
end subroutine add_block_d

subroutine add_block_z(this, cur_row, cur_col_offset, in_nr, in_nc, block)
  type(RS_SparseMatrixZ), intent(inout) :: this
  integer, intent(in) :: cur_row, cur_col_offset
  integer, intent(in) :: in_nr, in_nc
  complex(dp), intent(in) :: block(:,:)

  integer ptr_i, this_nr, this_nc

  ptr_i = this%l%row_indices(cur_row)+cur_col_offset
  this_nr = this%l%block_size(cur_row)
  this_nc = this%l%block_size(this%l%col(ptr_i))
  if (in_nr /= this_nr .or. in_nc /= this_nc) then
    call System_abort("add_block_z tried to add block of wrong shape in_nr,nc "//in_nr// " "//in_nc//" sp_nr,nc "//this_nr//" "//this_nc)
  end if
  this%data(this%l%data_ptrs(ptr_i):this%l%data_ptrs(ptr_i)+this_nr*this_nc-1) = &
    this%data(this%l%data_ptrs(ptr_i):this%l%data_ptrs(ptr_i)+this_nr*this_nc-1) + &
    reshape( block(1:in_nr,1:in_nc), (/this_nr*this_nc/) )
end subroutine add_block_z

function RS_SparseMatrix_partial_TraceMult_dden_dsp(a, b, a_T, b_T, diag_mask, offdiag_mask) result(v)
  type(MatrixD), intent(in) :: a
  type(RS_SparseMatrixD), intent(in) ::  b
  logical, intent(in), optional :: a_T, b_T
  logical, intent(in), optional, target :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: v(a%N)

  logical u_a_T, u_b_T
  integer i, ji, j, tt
  integer block_ni, block_nj, max_block_size
  real(dp), allocatable :: a_block(:,:), b_block(:,:)
  logical, pointer :: d_mask(:), od_mask(:)
  logical use_a_ij_transpose
  logical have_any

  max_block_size = maxval(b%l%block_size)
  allocate(a_block(max_block_size,max_block_size))
  allocate(b_block(max_block_size,max_block_size))
  if (present(diag_mask)) then
    d_mask => diag_mask
  else
    allocate(d_mask(b%l%N))
    d_mask = .true.
  end if
  if (present(offdiag_mask)) then
    od_mask => offdiag_mask
  else
    allocate(od_mask(b%l%N))
    od_mask = .true.
  end if

  u_a_T = optional_default(.false., a_T)
  u_b_T = optional_default(.false., b_T)

  use_a_ij_transpose = .false.
  if ((u_a_T .and. u_b_T) .or. (.not. u_a_T .and. .not. u_b_T)) use_a_ij_transpose = .true.

  v = 0.0_dp
  do i=1, b%l%N
    block_ni = b%l%block_size(i)
    do ji=b%l%row_indices(i), b%l%row_indices(i+1)-1
      j = b%l%col(ji)
      if ((i == j .and. .not. d_mask(i)) .or. ((i /= j) .and. .not. od_mask(i) .and. .not. od_mask(j))) cycle
      block_nj = b%l%block_size(j)
      have_any = get_dense_block(a%data, a%scalapack_info_obj, use_a_ij_transpose, i, j, b%l, block_ni, block_nj, a_block)

      if (have_any) then
	b_block(1:block_ni,1:block_nj) = reshape(b%data(b%l%data_ptrs(ji):b%l%data_ptrs(ji)+block_ni*block_nj-1),(/block_ni,block_nj/))

	if (u_b_T) then
	  do tt=1, block_nj
	    v(b%l%dense_row_of_row(i):b%l%dense_row_of_row(i)+block_ni-1) = &
	      v(b%l%dense_row_of_row(i):b%l%dense_row_of_row(i)+block_ni-1) + &
	      a_block(1:block_ni,tt)*b_block(1:block_ni,tt)
	  end do
	else
	  do tt=1, block_ni
	    v(b%l%dense_row_of_row(j):b%l%dense_row_of_row(j)+block_nj-1) = &
	      v(b%l%dense_row_of_row(j):b%l%dense_row_of_row(j)+block_nj-1) + &
	      a_block(tt,1:block_nj)*b_block(tt,1:block_nj)
	  end do
	end if
      end if

    end do
  end do

  deallocate(a_block)
  deallocate(b_block)
  if (.not.present(diag_mask)) deallocate(d_mask)
  if (.not.present(offdiag_mask)) deallocate(od_mask)

end function RS_SparseMatrix_partial_TraceMult_dden_dsp

function RS_SparseMatrix_partial_TraceMult_dden_zsp(a, b, a_T, b_H, diag_mask, offdiag_mask) result(v)
  type(MatrixD), intent(in) :: a
  type(RS_SparseMatrixZ), intent(in) ::  b
  logical, intent(in), optional :: a_T, b_H
  logical, intent(in), optional, target :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: v(a%N)

  logical u_a_T, u_b_H
  integer i, jj, j, tt
  integer block_ni, block_nj, max_block_size
  real(dp), allocatable :: a_block(:,:)
  complex(dp), allocatable :: b_block(:,:)
  logical, pointer :: d_mask(:), od_mask(:)
  logical :: use_a_ij_transpose
  logical have_any

  max_block_size = maxval(b%l%block_size)
  allocate(a_block(max_block_size,max_block_size))
  allocate(b_block(max_block_size,max_block_size))
  if (present(diag_mask)) then
    d_mask => diag_mask
  else
    allocate(d_mask(b%l%N))
    d_mask = .true.
  end if
  if (present(offdiag_mask)) then
    od_mask => offdiag_mask
  else
    allocate(od_mask(b%l%N))
    od_mask = .true.
  end if

  u_a_T = optional_default(.false., a_T)
  u_b_H = optional_default(.false., b_H)

  use_a_ij_transpose = .false.
  if ((u_a_T .and. u_b_H) .or. (.not. u_a_T .and. .not. u_b_H)) use_a_ij_transpose = .true.

  v = 0.0_dp
  do i=1, b%l%N
    block_ni = b%l%block_size(i)
    do jj=b%l%row_indices(i), b%l%row_indices(i+1)-1
      j = b%l%col(jj)
      if ((i == j .and. .not. d_mask(i)) .or. ((i /= j) .and. .not. od_mask(i) .and. .not. od_mask(j))) cycle
      block_nj = b%l%block_size(j)
      have_any = get_dense_block(a%data, a%scalapack_info_obj, use_a_ij_transpose, i, j, b%l, block_ni, block_nj, a_block)

      if (have_any) then
	b_block(1:block_ni,1:block_nj) = reshape(b%data(b%l%data_ptrs(jj):b%l%data_ptrs(jj)+block_ni*block_nj-1),(/block_ni,block_nj/))
	if (u_b_H) b_block = conjg(b_block)

	if (u_b_H) then
	  do tt=1, block_nj
	    v(b%l%dense_row_of_row(i):b%l%dense_row_of_row(i)+block_ni-1) = &
	      v(b%l%dense_row_of_row(i):b%l%dense_row_of_row(i)+block_ni-1) + &
	      a_block(1:block_ni,tt)*b_block(1:block_ni,tt)
	  end do
	else
	  do tt=1, block_ni
	    v(b%l%dense_row_of_row(j):b%l%dense_row_of_row(j)+block_nj-1) = &
	      v(b%l%dense_row_of_row(j):b%l%dense_row_of_row(j)+block_nj-1) + &
	      a_block(tt,1:block_nj)*b_block(tt,1:block_nj)
	  end do
	end if
      end if

    end do
  end do

  deallocate(a_block)
  deallocate(b_block)
  if (.not.present(diag_mask)) deallocate(d_mask)
  if (.not.present(offdiag_mask)) deallocate(od_mask)

end function RS_SparseMatrix_partial_TraceMult_dden_zsp

function RS_SparseMatrix_partial_TraceMult_zden_dsp(a, b, a_H, b_T, diag_mask, offdiag_mask) result(v)
  type(MatrixZ), intent(in) :: a
  type(RS_SparseMatrixD), intent(in) ::  b
  logical, intent(in), optional :: a_H, b_T
  logical, intent(in), optional, target :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: v(a%N)

  logical u_a_H, u_b_T
  integer i, jj, j, tt
  integer block_ni, block_nj, max_block_size
  complex(dp), allocatable :: a_block(:,:)
  real(dp), allocatable :: b_block(:,:)
  logical, pointer :: d_mask(:), od_mask(:)
  logical :: use_a_ij_transpose
  logical have_any

  max_block_size = maxval(b%l%block_size)
  allocate(a_block(max_block_size,max_block_size))
  allocate(b_block(max_block_size,max_block_size))
  if (present(diag_mask)) then
    d_mask => diag_mask
  else
    allocate(d_mask(b%l%N))
    d_mask = .true.
  end if
  if (present(offdiag_mask)) then
    od_mask => offdiag_mask
  else
    allocate(od_mask(b%l%N))
    od_mask = .true.
  end if

  u_a_H = optional_default(.false., a_H)
  u_b_T = optional_default(.false., b_T)

  use_a_ij_transpose = .false.
  if ((u_a_H .and. u_b_T) .or. (.not. u_a_H .and. .not. u_b_T)) use_a_ij_transpose = .true.

  v = 0.0_dp
  do i=1, b%l%N
    block_ni = b%l%block_size(i)
    do jj=b%l%row_indices(i), b%l%row_indices(i+1)-1
      j = b%l%col(jj)
      if ((i == j .and. .not. d_mask(i)) .or. ((i /= j) .and. .not. od_mask(i) .and. .not. od_mask(j))) cycle
      block_nj = b%l%block_size(j)
      have_any = get_dense_block(a%data, a%scalapack_info_obj, use_a_ij_transpose, i, j, b%l, block_ni, block_nj, a_block)
      if (u_a_H) a_block = conjg(a_block)

      if (have_any) then
	b_block(1:block_ni,1:block_nj) = reshape(b%data(b%l%data_ptrs(jj):b%l%data_ptrs(jj)+block_ni*block_nj-1),(/block_ni,block_nj/))

	if (u_b_T) then
	  do tt=1, block_nj
	    v(b%l%dense_row_of_row(i):b%l%dense_row_of_row(i)+block_ni-1) = &
	      v(b%l%dense_row_of_row(i):b%l%dense_row_of_row(i)+block_ni-1) + &
	      a_block(1:block_ni,tt)*b_block(1:block_ni,tt)
	  end do
	else
	  do tt=1, block_ni
	    v(b%l%dense_row_of_row(j):b%l%dense_row_of_row(j)+block_nj-1) = &
	      v(b%l%dense_row_of_row(j):b%l%dense_row_of_row(j)+block_nj-1) + &
	      a_block(tt,1:block_nj)*b_block(tt,1:block_nj)
	  end do
	end if
      end if

    end do
  end do

  deallocate(a_block)
  deallocate(b_block)
  if (.not.present(diag_mask)) deallocate(d_mask)
  if (.not.present(offdiag_mask)) deallocate(od_mask)

end function RS_SparseMatrix_partial_TraceMult_zden_dsp

function RS_SparseMatrix_partial_TraceMult_zden_zsp(a, b, a_H, b_H, diag_mask, offdiag_mask) result(v)
  type(MatrixZ), intent(in) :: a
  type(RS_SparseMatrixZ), intent(in) ::  b
  logical, intent(in), optional :: a_H, b_H
  logical, intent(in), optional, target :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: v(a%N)

  logical u_a_H, u_b_H
  integer i, jj, j, tt
  integer block_ni, block_nj, max_block_size
  complex(dp), allocatable :: a_block(:,:), b_block(:,:)
  logical, pointer :: d_mask(:), od_mask(:)
  logical :: use_a_ij_transpose
  logical have_any

  max_block_size = maxval(b%l%block_size)
  allocate(a_block(max_block_size,max_block_size))
  allocate(b_block(max_block_size,max_block_size))
  if (present(diag_mask)) then
    d_mask => diag_mask
  else
    allocate(d_mask(b%l%N))
    d_mask = .true.
  end if
  if (present(offdiag_mask)) then
    od_mask => offdiag_mask
  else
    allocate(od_mask(b%l%N))
    od_mask = .true.
  end if

  u_a_H = optional_default(.false., a_H)
  u_b_H = optional_default(.false., b_H)

  use_a_ij_transpose = .false.
  if ((u_a_H .and. u_b_H) .or. (.not. u_a_H .and. .not. u_b_H)) use_a_ij_transpose = .true.

  v = 0.0_dp
  do i=1, b%l%N
    block_ni = b%l%block_size(i)
    do jj=b%l%row_indices(i), b%l%row_indices(i+1)-1
      j = b%l%col(jj)
      if ((i == j .and. .not. d_mask(i)) .or. ((i /= j) .and. .not. od_mask(i) .and. .not. od_mask(j))) cycle
      block_nj = b%l%block_size(j)
      have_any = get_dense_block(a%data, a%scalapack_info_obj, use_a_ij_transpose, i, j, b%l, block_ni, block_nj, a_block)
      if (u_a_H) a_block = conjg(a_block)

      if (have_any) then
	b_block(1:block_ni,1:block_nj) = reshape(b%data(b%l%data_ptrs(jj):b%l%data_ptrs(jj)+block_ni*block_nj-1),(/block_ni,block_nj/))
	if (u_b_H) b_block = conjg(b_block)

	if (u_b_H) then
	  do tt=1, block_nj
	    v(b%l%dense_row_of_row(i):b%l%dense_row_of_row(i)+block_ni-1) = &
	      v(b%l%dense_row_of_row(i):b%l%dense_row_of_row(i)+block_ni-1) + &
	      a_block(1:block_ni,tt)*b_block(1:block_ni,tt)
	  end do
	else
	  do tt=1, block_ni
	    v(b%l%dense_row_of_row(j):b%l%dense_row_of_row(j)+block_nj-1) = &
	      v(b%l%dense_row_of_row(j):b%l%dense_row_of_row(j)+block_nj-1) + &
	      a_block(tt,1:block_nj)*b_block(tt,1:block_nj)
	  end do
	end if
      end if

    end do
  end do

  deallocate(a_block)
  deallocate(b_block)
  if (.not.present(diag_mask)) deallocate(d_mask)
  if (.not.present(offdiag_mask)) deallocate(od_mask)

end function RS_SparseMatrix_partial_TraceMult_zden_zsp

function RS_SparseMatrix_partial_TraceMult_dsp_zden(a, b, a_T, b_H, diag_mask, offdiag_mask) result(v)
  type(RS_SparseMatrixD), intent(in) ::  a
  type(MatrixZ), intent(in) :: b
  logical, intent(in), optional :: a_T, b_H
  logical, intent(in), optional, target :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: v(b%N)

  logical u_a_T, u_b_H
  integer i, jj, j, tt
  integer block_ni, block_nj, max_block_size
  real(dp), allocatable :: a_block(:,:)
  complex(dp), allocatable :: b_block(:,:)
  logical, pointer :: d_mask(:), od_mask(:)
  logical use_b_ij_transpose
  logical have_any

  max_block_size = maxval(a%l%block_size)
  allocate(b_block(max_block_size,max_block_size))
  allocate(a_block(max_block_size,max_block_size))
  if (present(diag_mask)) then
    d_mask => diag_mask
  else
    allocate(d_mask(a%l%N))
    d_mask = .true.
  end if
  if (present(offdiag_mask)) then
    od_mask => offdiag_mask
  else
    allocate(od_mask(a%l%N))
    od_mask = .true.
  end if

  u_a_T = optional_default(.false., a_T)
  u_b_H = optional_default(.false., b_H)

  use_b_ij_transpose = .false.
  if ((u_a_T .and. u_b_H) .or. (.not. u_a_T .and. .not. u_b_H)) use_b_ij_transpose = .true.

  v = 0.0_dp
  do i=1, a%l%N
    block_ni = a%l%block_size(i)
    do jj=a%l%row_indices(i), a%l%row_indices(i+1)-1
      j = a%l%col(jj)
      if ((i == j .and. .not. d_mask(i)) .or. ((i /= j) .and. .not. od_mask(i) .and. .not. od_mask(j))) cycle
      block_nj = a%l%block_size(j)
      have_any = get_dense_block(b%data, b%scalapack_info_obj, use_b_ij_transpose, i, j, a%l, block_ni, block_nj, b_block)
      if (u_b_H) b_block = conjg(b_block)

      if (have_any) then
	a_block(1:block_ni,1:block_nj) = reshape(a%data(a%l%data_ptrs(jj):a%l%data_ptrs(jj)+block_ni*block_nj-1),(/block_ni,block_nj/))

	if (u_a_T) then
	  do tt=1, block_ni
	    v(a%l%dense_row_of_row(j):a%l%dense_row_of_row(j)+block_nj-1) = &
	      v(a%l%dense_row_of_row(j):a%l%dense_row_of_row(j)+block_nj-1) + &
	      a_block(tt,1:block_nj)*b_block(tt,1:block_nj)
	  end do
	else
	  do tt=1, block_nj
	    v(a%l%dense_row_of_row(i):a%l%dense_row_of_row(i)+block_ni-1) = &
	      v(a%l%dense_row_of_row(i):a%l%dense_row_of_row(i)+block_ni-1) + &
	      a_block(1:block_ni,tt)*b_block(1:block_ni,tt)
	  end do
	end if
      end if

    end do
  end do

  deallocate(b_block)
  deallocate(a_block)
  if (.not.present(diag_mask)) deallocate(d_mask)
  if (.not.present(offdiag_mask)) deallocate(od_mask)

end function RS_SparseMatrix_partial_TraceMult_dsp_zden

function RS_SparseMatrix_partial_TraceMult_zsp_zden(a, b, a_H, b_H, diag_mask, offdiag_mask) result(v)
  type(RS_SparseMatrixZ), intent(in) ::  a
  type(MatrixZ), intent(in) :: b
  logical, intent(in), optional :: a_H, b_H
  logical, intent(in), optional, target :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: v(b%N)

  logical u_a_H, u_b_H
  integer i, jj, j, tt
  integer block_ni, block_nj, max_block_size
  complex(dp), allocatable :: a_block(:,:)
  complex(dp), allocatable :: b_block(:,:)
  logical, pointer :: d_mask(:), od_mask(:)
  logical use_b_ij_transpose
  logical have_any

  max_block_size = maxval(a%l%block_size)
  allocate(b_block(max_block_size,max_block_size))
  allocate(a_block(max_block_size,max_block_size))
  if (present(diag_mask)) then
    d_mask => diag_mask
  else
    allocate(d_mask(a%l%N))
    d_mask = .true.
  end if
  if (present(offdiag_mask)) then
    od_mask => offdiag_mask
  else
    allocate(od_mask(a%l%N))
    od_mask = .true.
  end if

  u_a_H = optional_default(.false., a_H)
  u_b_H = optional_default(.false., b_H)

  use_b_ij_transpose = .false.
  if ((u_a_H .and. u_b_H) .or. (.not. u_a_H .and. .not. u_b_H)) use_b_ij_transpose = .true.

  v = 0.0_dp
  do i=1, a%l%N
    block_ni = a%l%block_size(i)
    do jj=a%l%row_indices(i), a%l%row_indices(i+1)-1
      j = a%l%col(jj)
      if ((i == j .and. .not. d_mask(i)) .or. ((i /= j) .and. .not. od_mask(i) .and. .not. od_mask(j))) cycle
      block_nj = a%l%block_size(j)
      have_any = get_dense_block(b%data, b%scalapack_info_obj, use_b_ij_transpose, i, j, a%l, block_ni, block_nj, b_block)
      if (u_b_H) b_block = conjg(b_block)

      if (have_any) then
	a_block(1:block_ni,1:block_nj) = reshape(a%data(a%l%data_ptrs(jj):a%l%data_ptrs(jj)+block_ni*block_nj-1),(/block_ni,block_nj/))
	if (u_a_H) a_block = conjg(a_block)

	if (u_a_H) then
	  do tt=1, block_ni
	    v(a%l%dense_row_of_row(j):a%l%dense_row_of_row(j)+block_nj-1) = &
	      v(a%l%dense_row_of_row(j):a%l%dense_row_of_row(j)+block_nj-1) + &
	      a_block(tt,1:block_nj)*b_block(tt,1:block_nj)
	  end do
	else
	  do tt=1, block_nj
	    v(a%l%dense_row_of_row(i):a%l%dense_row_of_row(i)+block_ni-1) = &
	      v(a%l%dense_row_of_row(i):a%l%dense_row_of_row(i)+block_ni-1) + &
	      a_block(1:block_ni,tt)*b_block(1:block_ni,tt)
	  end do
	end if
      end if

    end do
  end do

  deallocate(b_block)
  deallocate(a_block)
  if (.not.present(diag_mask)) deallocate(d_mask)
  if (.not.present(offdiag_mask)) deallocate(od_mask)

end function RS_SparseMatrix_partial_TraceMult_zsp_zden

function RS_SparseMatrix_partial_TraceMult_dsp_dden(a, b, a_T, b_T, diag_mask, offdiag_mask) result(v)
  type(RS_SparseMatrixD), intent(in) ::  a
  type(MatrixD), intent(in) :: b
  logical, intent(in), optional :: a_T, b_T
  logical, intent(in), optional, target :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: v(b%N)

  logical u_a_T, u_b_T
  integer i, jj, j, tt
  integer block_ni, block_nj, max_block_size
  real(dp), allocatable :: a_block(:,:)
  real(dp), allocatable :: b_block(:,:)
  logical, pointer :: d_mask(:), od_mask(:)
  logical use_b_ij_transpose
  logical have_any

  max_block_size = maxval(a%l%block_size)
  allocate(b_block(max_block_size,max_block_size))
  allocate(a_block(max_block_size,max_block_size))
  if (present(diag_mask)) then
    d_mask => diag_mask
  else
    allocate(d_mask(a%l%N))
    d_mask = .true.
  end if
  if (present(offdiag_mask)) then
    od_mask => offdiag_mask
  else
    allocate(od_mask(a%l%N))
    od_mask = .true.
  end if

  u_a_T = optional_default(.false., a_T)
  u_b_T = optional_default(.false., b_T)

  use_b_ij_transpose = .false.
  if ((u_a_T .and. u_b_T) .or. (.not. u_a_T .and. .not. u_b_T)) use_b_ij_transpose = .true.

  v = 0.0_dp
  do i=1, a%l%N
    block_ni = a%l%block_size(i)
    do jj=a%l%row_indices(i), a%l%row_indices(i+1)-1
      j = a%l%col(jj)
      if ((i == j .and. .not. d_mask(i)) .or. ((i /= j) .and. .not. od_mask(i) .and. .not. od_mask(j))) cycle
      block_nj = a%l%block_size(j)
      have_any = get_dense_block(b%data, b%scalapack_info_obj, use_b_ij_transpose, i, j, a%l, block_ni, block_nj, b_block)

      if (have_any) then
	a_block(1:block_ni,1:block_nj) = reshape(a%data(a%l%data_ptrs(jj):a%l%data_ptrs(jj)+block_ni*block_nj-1),(/block_ni,block_nj/))

	if (u_a_T) then
	  do tt=1, block_ni
	    v(a%l%dense_row_of_row(j):a%l%dense_row_of_row(j)+block_nj-1) = &
	      v(a%l%dense_row_of_row(j):a%l%dense_row_of_row(j)+block_nj-1) + &
	      a_block(tt,1:block_nj)*b_block(tt,1:block_nj)
	  end do
	else
	  do tt=1, block_nj
	    v(a%l%dense_row_of_row(i):a%l%dense_row_of_row(i)+block_ni-1) = &
	      v(a%l%dense_row_of_row(i):a%l%dense_row_of_row(i)+block_ni-1) + &
	      a_block(1:block_ni,tt)*b_block(1:block_ni,tt)
	  end do
	end if
      end if

    end do
  end do

  deallocate(b_block)
  deallocate(a_block)
  if (.not.present(diag_mask)) deallocate(d_mask)
  if (.not.present(offdiag_mask)) deallocate(od_mask)

end function RS_SparseMatrix_partial_TraceMult_dsp_dden

function get_dense_blockD(a_data, a_scalapack_info_obj, use_a_ij_transpose, i, j, b_l, block_ni, block_nj, a_block)
  real(dp), intent(in) :: a_data(:,:)
  type(Matrix_ScaLAPACK_Info), intent(in) :: a_scalapack_info_obj
  logical, intent(in) :: use_a_ij_transpose
  integer, intent(in) :: i, j
  type(RS_SparseMatrixL), intent(in) :: b_l
  integer, intent(in) :: block_ni, block_nj
  real(dp), intent(out) :: a_block(:,:)
  logical :: get_dense_blockD

  integer io, jo, g_i, g_j, l_i, l_j

  if (a_scalapack_info_obj%active) then
    get_dense_blockD = .false.
    do io = 0, block_ni-1
    do jo = 0, block_nj-1
      g_i = b_l%dense_row_of_row(i) + io
      g_j = b_l%dense_row_of_row(j) + jo
      if (use_a_ij_transpose) then
	call coords_global_to_local(a_scalapack_info_obj, g_j, g_i, l_i, l_j)
	if (l_i > 0 .and. l_j > 0) then
	  a_block(io+1,jo+1) = a_data(l_i,l_j)
	  get_dense_blockD = .true.
	else
	  a_block(io+1,jo+1) = 0.0_dp
	end if
      else
	call coords_global_to_local(a_scalapack_info_obj, g_i, g_j, l_i, l_j)
	if (l_i > 0 .and. l_j > 0) then
	  a_block(io+1,jo+1) = a_data(l_i,l_j)
	  get_dense_blockD = .true.
	else
	  a_block(io+1,jo+1) = 0.0_dp
	end if
      end if
    end do
    end do
  else
    get_dense_blockD = .true.
    if (use_a_ij_transpose) then
      a_block(1:block_ni,1:block_nj) = transpose(a_data(b_l%dense_row_of_row(j):b_l%dense_row_of_row(j)+block_nj-1, &
				                        b_l%dense_row_of_row(i):b_l%dense_row_of_row(i)+block_ni-1))
    else
      a_block(1:block_ni,1:block_nj) = a_data(b_l%dense_row_of_row(i):b_l%dense_row_of_row(i)+block_ni-1, &
		                              b_l%dense_row_of_row(j):b_l%dense_row_of_row(j)+block_nj-1)
    end if
  end if

end function get_dense_blockD

function get_dense_blockZ(a_data, a_scalapack_info_obj, use_a_ij_transpose, i, j, b_l, block_ni, block_nj, a_block)
  complex(dp), intent(in) :: a_data(:,:)
  type(Matrix_ScaLAPACK_Info), intent(in) :: a_scalapack_info_obj
  logical, intent(in) :: use_a_ij_transpose
  integer, intent(in) :: i, j
  type(RS_SparseMatrixL), intent(in) :: b_l
  integer, intent(in) :: block_ni, block_nj
  complex(dp), intent(out) :: a_block(:,:)
  logical :: get_dense_blockZ

  integer io, jo, g_i, g_j, l_i, l_j

  if (a_scalapack_info_obj%active) then
    get_dense_blockZ = .false.
    do io = 0, block_ni-1
    do jo = 0, block_nj-1
      g_i = b_l%dense_row_of_row(i) + io
      g_j = b_l%dense_row_of_row(j) + jo
      if (use_a_ij_transpose) then
	call coords_global_to_local(a_scalapack_info_obj, g_j, g_i, l_i, l_j)
	if (l_i > 0 .and. l_j > 0) then
	  a_block(io+1,jo+1) = a_data(l_i,l_j)
	  get_dense_blockZ = .true.
	else
	  a_block(io+1,jo+1) = 0.0_dp
	end if
      else
	call coords_global_to_local(a_scalapack_info_obj, g_i, g_j, l_i, l_j)
	if (l_i > 0 .and. l_j > 0) then
	  a_block(io+1,jo+1) = a_data(l_i,l_j)
	  get_dense_blockZ = .true.
	else
	  a_block(io+1,jo+1) = 0.0_dp
	end if
      end if
    end do
    end do
  else
    get_dense_blockZ = .true.
    if (use_a_ij_transpose) then
      a_block(1:block_ni,1:block_nj) = transpose(a_data(b_l%dense_row_of_row(j):b_l%dense_row_of_row(j)+block_nj-1, &
				                        b_l%dense_row_of_row(i):b_l%dense_row_of_row(i)+block_ni-1))
    else
      a_block(1:block_ni,1:block_nj) = a_data(b_l%dense_row_of_row(i):b_l%dense_row_of_row(i)+block_ni-1, &
		                              b_l%dense_row_of_row(j):b_l%dense_row_of_row(j)+block_nj-1)
    end if
  end if

end function get_dense_blockZ

function RS_SparseMatrix_TraceMult_dden_dsp(a, b, w, a_T, b_T, diag_mask, offdiag_mask)
  type(MatrixD), intent(in) :: a
  type(RS_SparseMatrixD), intent(in) ::  b
  real(dp), intent(in), optional :: w(:)
  logical, intent(in), optional :: a_T, b_T
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)
  real(dp) :: RS_SparseMatrix_TraceMult_dden_dsp

  if (present(w)) then
    RS_SparseMatrix_TraceMult_dden_dsp = sum(partial_TraceMult(a, b, a_T, b_T, diag_mask, offdiag_mask)*w)
  else
    RS_SparseMatrix_TraceMult_dden_dsp = sum(partial_TraceMult(a, b, a_T, b_T, diag_mask, offdiag_mask))
  end if

end function RS_SparseMatrix_TraceMult_dden_dsp

function RS_SparseMatrix_TraceMult_dden_zsp(a, b, w, a_T, b_H, diag_mask, offdiag_mask)
  type(MatrixD), intent(in) :: a
  type(RS_SparseMatrixZ), intent(in) ::  b
  real(dp), intent(in), optional :: w(:)
  logical, intent(in), optional :: a_T, b_H
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: RS_SparseMatrix_TraceMult_dden_zsp

  if (present(w)) then
    RS_SparseMatrix_TraceMult_dden_zsp = sum(partial_TraceMult(a, b, a_T, b_H, diag_mask, offdiag_mask)*w)
  else
    RS_SparseMatrix_TraceMult_dden_zsp = sum(partial_TraceMult(a, b, a_T, b_H, diag_mask, offdiag_mask))
  end if

end function RS_SparseMatrix_TraceMult_dden_zsp

function RS_SparseMatrix_TraceMult_zden_dsp(a, b, w, a_H, b_T, diag_mask, offdiag_mask)
  type(MatrixZ), intent(in) :: a
  type(RS_SparseMatrixD), intent(in) ::  b
  real(dp), intent(in), optional :: w(:)
  logical, intent(in), optional :: a_H, b_T
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: RS_SparseMatrix_TraceMult_zden_dsp

  if (present(w)) then
    RS_SparseMatrix_TraceMult_zden_dsp = sum(partial_TraceMult(a, b, a_H, b_T, diag_mask, offdiag_mask)*w)
  else
    RS_SparseMatrix_TraceMult_zden_dsp = sum(partial_TraceMult(a, b, a_H, b_T, diag_mask, offdiag_mask))
  end if

end function RS_SparseMatrix_TraceMult_zden_dsp

function RS_SparseMatrix_TraceMult_dsp_zden(a, b, w, a_T, b_H, diag_mask, offdiag_mask)
  type(RS_SparseMatrixD), intent(in) ::  a
  type(MatrixZ), intent(in) :: b
  real(dp), intent(in), optional :: w(:)
  logical, intent(in), optional :: a_T, b_H
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: RS_SparseMatrix_TraceMult_dsp_zden

  if (present(w)) then
    RS_SparseMatrix_TraceMult_dsp_zden = sum(partial_TraceMult(a, b, a_T, b_H, diag_mask, offdiag_mask)*w)
  else
    RS_SparseMatrix_TraceMult_dsp_zden = sum(partial_TraceMult(a, b, a_T, b_H, diag_mask, offdiag_mask))
  end if

end function RS_SparseMatrix_TraceMult_dsp_zden

function RS_SparseMatrix_TraceMult_dsp_dden(a, b, w, a_T, b_T, diag_mask, offdiag_mask)
  type(RS_SparseMatrixD), intent(in) ::  a
  type(MatrixD), intent(in) :: b
  real(dp), intent(in), optional :: w(:)
  logical, intent(in), optional :: a_T, b_T
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: RS_SparseMatrix_TraceMult_dsp_dden

  if (present(w)) then
    RS_SparseMatrix_TraceMult_dsp_dden = sum(partial_TraceMult(a, b, a_T, b_T, diag_mask, offdiag_mask)*w)
  else
    RS_SparseMatrix_TraceMult_dsp_dden = sum(partial_TraceMult(a, b, a_T, b_T, diag_mask, offdiag_mask))
  end if

end function RS_SparseMatrix_TraceMult_dsp_dden

function RS_SparseMatrix_TraceMult_zsp_zden(a, b, w, a_H, b_H, diag_mask, offdiag_mask)
  type(RS_SparseMatrixZ), intent(in) ::  a
  type(MatrixZ), intent(in) :: b
  real(dp), intent(in), optional :: w(:)
  logical, intent(in), optional :: a_H, b_H
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: RS_SparseMatrix_TraceMult_zsp_zden

  if (present(w)) then
    RS_SparseMatrix_TraceMult_zsp_zden = sum(partial_TraceMult(a, b, a_H, b_H, diag_mask, offdiag_mask)*w)
  else
    RS_SparseMatrix_TraceMult_zsp_zden = sum(partial_TraceMult(a, b, a_H, b_H, diag_mask, offdiag_mask))
  end if

end function RS_SparseMatrix_TraceMult_zsp_zden

function RS_SparseMatrix_TraceMult_zden_zsp(a, b, w, a_H, b_H, diag_mask, offdiag_mask)
  type(MatrixZ), intent(in) :: a
  type(RS_SparseMatrixZ), intent(in) ::  b
  real(dp), intent(in), optional :: w(:)
  logical, intent(in), optional :: a_H, b_H
  logical, intent(in), optional :: diag_mask(:), offdiag_mask(:)
  complex(dp) :: RS_SparseMatrix_TraceMult_zden_zsp

  if (present(w)) then
    RS_SparseMatrix_TraceMult_zden_zsp = sum(partial_TraceMult(a, b, a_H, b_H, diag_mask, offdiag_mask)*w)
  else
    RS_SparseMatrix_TraceMult_zden_zsp = sum(partial_TraceMult(a, b, a_H, b_H, diag_mask, offdiag_mask))
  end if

end function RS_SparseMatrix_TraceMult_zden_zsp

subroutine RS_SparseMatrixD_multDiagRL_d(this, A, diag)
  type(RS_SparseMatrixD), intent(inout) :: this
  type(RS_SparseMatrixD), intent(in) :: A
  real(dp), intent(in) :: diag(:)

  integer i, jj, j, block_ni, block_nj, io, jo
  real(dp) fac

  do i=1, A%l%N
    block_ni = A%l%block_size(i)
    do jj=A%l%row_indices(i), A%l%row_indices(i+1)-1
      j = A%l%col(jj)
      block_nj = A%l%block_size(j)
      do jo=0, block_nj-1
      do io=0, block_ni-1
	fac = 0.5_dp*(diag(this%l%dense_row_of_row(i)+io)+diag(this%l%dense_row_of_row(j)+jo))
	this%data(this%l%data_ptrs(jj)+io+block_ni*jo) = fac*A%data(this%l%data_ptrs(jj)+io+block_ni*jo)
      end do
      end do
    end do
  end do

end subroutine RS_SparseMatrixD_multDiagRL_d

subroutine RS_SparseMatrixZ_multDiagRL_d(this, A, diag)
  type(RS_SparseMatrixZ), intent(inout) :: this
  type(RS_SparseMatrixZ), intent(in) :: A
  real(dp), intent(in) :: diag(:)

  integer i, jj, j, block_ni, block_nj, io, jo
  real(dp) fac

  do i=1, A%l%N
    block_ni = A%l%block_size(i)
    do jj=A%l%row_indices(i), A%l%row_indices(i+1)-1
      j = A%l%col(jj)
      block_nj = A%l%block_size(j)
      do jo=0, block_nj-1
      do io=0, block_ni-1
	fac = 0.5_dp*(diag(this%l%dense_row_of_row(i)+io)+diag(this%l%dense_row_of_row(j)+jo))
	this%data(this%l%data_ptrs(jj)+io+block_ni*jo) = fac*A%data(this%l%data_ptrs(jj)+io+block_ni*jo)
      end do
      end do
    end do
  end do

end subroutine RS_SparseMatrixZ_multDiagRL_d

subroutine check_sparse_layout(l)
  type(RS_SparseMatrixL), intent(in) :: l
  integer i, ji, j

  do i=1, l%N
    do ji=l%row_indices(i), l%row_indices(i+1)-1
      j = l%col(ji)

      call print("check_sparse " // i // " " // j)

      if (ji > l%row_indices(i)) then
	if (j <= l%col(ji-1)) then
	  call print("ERROR row " // i // " col " // j // " (ji="//ji//") out of order")
	  call print(l%col(l%row_indices(i):l%row_indices(i+1)-1))
	end if
      end if

      if (.not. any(l%col(l%row_indices(j):l%row_indices(j+1)-1) == i)) then
	call print ("ERROR row "//i//" col "//j//" has no transpose match")
	call print(l%col(l%row_indices(j):l%row_indices(j+1)-1))
      end if

    end do
  end do
end subroutine

end module RS_SparseMatrix_module

max_stars_repo_name: albapa/QUIP

max_issues_repo_name: albapa/QUIP

max_forks_repo_name: albapa/QUIP


text:
!********o*********o*********o*********o*********o*********o*********o**
! This program allows to players to play a game that is intended to
! be a crude model of basketball.  The idea was originally implemented
! on programmable calculators in the early days of that technology.
! The simulation uses fortran's pseudo-random number generator.
!
! It's purpose for this class is as a demonstration of a 'menu-driven'
! user interface.
!
! Inputs
!   A series of choices from the players.
!
! Outputs
!   A series of status updates shown to the players.
!---------------------
! Author:   Bob Harris
! Class:    CMPSC 201F
! Section:  1234
! Date Due: 29 Oct 2003
!********o*********o*********o*********o*********o*********o*********o**

program hoops1977

use, intrinsic :: iso_fortran_env, only : stdin=>input_unit
implicit none

! game configuration

integer, parameter :: gameSeconds     = 40*60   ! 40 minute game
integer, parameter :: overtimeSeconds =  5*60   ! 5 minute overtime(s)

real, parameter    :: stealChance     = 0.25    ! chance pass is stolen
real, parameter    :: reboundChance   = 0.10    ! chance of rebound

! current 'state' of the game

integer :: whosBall         ! which player has the ball (1 or 2)
integer :: score1, score2   ! current score for each player
integer :: secondsLeft      ! number of seconds remaining in the game
integer :: mins, secs       ! working variables to display clock

real    :: shotChance       ! chance of making current shot (0..1)
integer :: pct              ! working variables to display chance

! player commands

character :: command          ! the latest command from the current player

character, parameter :: pass = "p",  shoot = "s"

! other variables
integer :: ierr
real    :: r                ! a random number, 0.0 <= r < 1.0

logical :: takeTheShot      ! true  => take the shot we have
                            ! false => don't
logical :: lookForShot      ! true  => pass ball to get open shot
                            ! false => don't

call random_init(.false., .false.)

!initialize the score and time
score1 = 0
score2 = 0
secondsLeft = gameSeconds
! jump ball
call RANDOM_NUMBER(r)
if (r < 0.5) then
    whosBall = 1
else
    whosBall = 2
end if
print '(1X,A,I1,A)', "Player ", whosBall, " wins the tip"
lookForShot = .true.
call random_number(shotChance)

!----------
! continually loop, processing commands until players decide to quit
!----------

main : do

  if (secondsLeft <= 0) then

    if (score1 == score2) then
      print *,       "Overtime!"
      secondsLeft    = overTimeSeconds
    else
      call gameIsOver()
    endif

  else
    mins = secondsLeft / 60
    secs = secondsLeft - (mins*60)
    pct  = nint(100 * shotChance) ! round to nearest %

    print *,                  "--------------"
    print '(1X,A,I3)',        "Player 1 score: ", score1
    print '(1X,A,I3)',        "Player 2 score: ", score2
    print '(1X,A,I3,A,I2.2)', "Time Left:      ", mins,":",secs
    print '(1X,A,I3,A)',      "You are open for a ",pct,"% shot"
    print '(1X,A,I1,A)',      "Player *", whosBall, "*, what shall we do?"
    print '(3X,A)', "q: take my ball and go home"
    print '(3X,A1,A)', pass,     ": look for another shot"
    print '(3X,A1,A)', shoot,    ": take that shot"
    read(stdin,'(A1)', iostat=ierr) command
  end if

  if (command == "q" .or. IS_IOSTAT_END(ierr)) call gameIsOver()
  if(ierr /= 0) cycle main

  !----------
  ! perform the command
  !----------

  takeTheShot = .false.
  lookForShot = .false.

  select case (command)
    case default
      cycle main
    case (pass)
      lookForShot = .true.
    case (shoot)
      takeTheShot = .true.
  end select

  !----------
  ! do the simulation
  !----------

  ! take a shot, if we're supposed to

  if (takeTheShot) then
      ! see if the shot is made
    call RANDOM_NUMBER(r)
    if (r < shotChance) then
      ! shot was made, count it on the score board

      print *, "Nice shot"
      if (whosBall == 1) then
          score1 = score1 + 2
      else
          score2 = score2 + 2
      end if

      whosBall = 3 - whosBall ! give ball to other player
    else
      ! shot was missed

      print *, "What a brick!"

      ! see if shooter gets own rebound

      call RANDOM_NUMBER(r)
      if (r < reboundChance) then
          print *, "But you got your own rebound"
      else
          whosBall = 3 - whosBall ! give ball to other player
      end if

    end if

    ! look for the next shot

    lookForShot = .true.
  end if

  ! pass the ball around until someone has an open shot

  look: do while (lookForShot)
    ! run some time off the clock

    call RANDOM_NUMBER(r)
    secondsLeft = secondsLeft - INT(12 + 33*r)
    if (secondsLeft <= 0) exit look

    ! see if the ball was stolen

    call RANDOM_NUMBER(r)
    if (r < stealChance) then
      whosBall = 3 - whosBall ! give ball to other player
      print '(1X,A,I1,A)', "Player ", whosBall, " steals the ball"
      cycle look ! go try for an open shot for this player
    end if

    ! figure out chance of making this shot

    call RANDOM_NUMBER(shotChance)
    lookForShot = .false.

  end do look

end do main

contains

subroutine gameIsOver()

print *, "(sound of buzzer)"
print '(1X,A,I3)', "Player 1 score: ", score1
print '(1X,A,I3)', "Player 2 score: ", score2
if (score1 > score2) then
  stop "Congratulations player 1"
elseif (score1 < score2) then
  stop "Congratulations player 2"
elseif (score1 == 0) then
  stop "No one scored!"
else
  stop "It's a tie!"
end if

end subroutine gameIsOver

end program hoops1977

max_stars_repo_name: fortran-gaming/basketball

max_issues_repo_name: fortran-gaming/basketball

max_forks_repo_name: fortran-gaming/basketball


text:
These are archived reviews of Nuggets from 2005. For more recent comments, please visit the Nugget entry.

20050726 14:47:34 nbsp Nugget is fine, but the selection of fresh food doesnt give me a woody. Users/RocksandDirt

20050812 19:23:34 nbsp Im very skeptical of the price savings scoreboard at Nugget stores. It doesnt explain what the numbers mean, and it just doesnt seem possible that Nugget could beat Safeway and Albertsons on price across the board. Produce and dairy prices are pretty competitive, and prices for many specialty items seem more than fair, but just about everything else (e.g. cereal, homewares, pasta sauce) seems to be priced much higher than at other stores. Users/RickEle

According to their website, the Price Challenge is that you bring in $100 of groceries from another store, then they scan them in and see if it would have been cheaper at Nugget.  If it isnt, they pay double the difference plus $10.  The scoreboard presumably tracks wins versus losses on that challenge.  Maybe the savings on produce and dairy make up for cereal and pasta?  Users/RyanSandler

20050813 10:53:43 nbsp The raw sausages they sell in a variety of different flavors are excellent when BBQd.  Of particular note in this regard is the pineapple teriyaki sausage. Users/AlexPomeranz

20050815 18:46:55 nbsp They raised the rice bowl price in the last few months. Used to be 2.99  which was a great deal.  Now its $3.99, except for Happy Hour 112  $2.99 during random days/weeks. Users/EdWins ES

20051004 21:47:02 nbsp While I have no information about prices at this store, back when I intensively did compare prices between a health food store and supermarkets I discovered that its critical which items you are comparing. The health food items were considerably cheaper in the health food store, and the unhealth food items were cheaper at the supermarkets. I assume this has to do with supplier volume. Users/SteveDavison

20051022 12:15:46 nbsp HORSESHIT They raised the price of their sandwich by a full DOLLAR (theyre now $6.50). Users/MattJojola

20051026 11:32:22 nbsp The panini sandwiches are tiny for $6.99, and not particularly wonderful. Users/ApolloStumpy

20051112 01:09:22 nbsp Always a wonderfully rare little oddity bit of food to bring home and try. Nugget wasnt that expensive really, and they always carried the strange ingredients I needed for more exotic cooking. Friendly staff, lots of cheese (woo!), and a fantastic selection natural or exotic foods. If theres one thing I miss about Davis, its Nugget. Users/TarynHornen

20051210 21:41:46 nbsp If you like brussels sprouts, try their brussels sprouts salad with pancetta and pine nuts. Its even better heated! Users/JohnnaSarkis

20051221 04:27:57 nbsp I think the Nugget is too expensive for me. That said, I happen to like their Chinese food. Late nights, at like, 9:30, theyll clean out the bins and put the remaining food into big containers to sell off for $3.75 each. Theyre pretty good, I think. SS

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!> Tests the Harmonic Bath module.
PROGRAM test_harmonic_bath
    USE harmonic_bath
    USE kinds
    USE input_output
    USE unit_conversions
    IMPLICIT NONE
    REAL(dp) :: beta
    REAL(dp), ALLOCATABLE :: test(:)
    REAL(dp), ALLOCATABLE :: x(:),p(:),xt(:),pt(:)
    nstate = 2 ; temperature = 77.0_dp
    beta = 1.0_dp / ( temperature * convert('kelvin','au_energy') )

    CALL initialize_bath
    ALLOCATE( test(nbath*nosc), x(nstate), p(nstate), xt(nstate),pt(nstate) )
    test=0.0_dp ; x=0.0_dp ; p=0.0_dp ; xt=0.0_dp ; pt=0.0_dp

    CALL sample_thermal_wigner(x_bath,p_bath,beta)

    IF (x_bath(1) == 0) THEN
        WRITE(*,*) 'Thermal bath not being sampled correctly!'
        STOP 1
    END IF

    test = bilinear_harmonic_force_pldm(x_bath,x,p,xt,pt)
    IF (test(1) == 0) THEN
        WRITE(*,*) 'Bilinear harmonic force pldm returning incorrect value'
        STOP 1
    END IF

    test = bilinear_harmonic_force_twa(x_bath,x,p)
    IF (test(1) == 0) THEN
        WRITE(*,*) 'Bilinear harmonic force twa returning incorrect value'
        STOP 1
    END IF

    CALL finalize_bath
    DEALLOCATE( test, x, p, xt, pt )

END PROGRAM test_harmonic_bath


text:
cc Copyright (C) 2010: Vladimir Rokhlin
cc Copyright (C) 2009-2012: Leslie Greengard and Zydrunas Gimbutas
cc Contact: greengard@cims.nyu.edu
cc
cc This software is being released under a modified FreeBSD license
cc (see COPYING in home directory).
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c    $Date: 2010-01-05 12:57:52 -0500 (Tue, 05 Jan 2010) $
c    $Revision: 782 $
c
c
c        Tensor product triangle quadrature rules for smooth functions.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
c
c
        subroutine triagauc(n,vert1,vert2,vert3,rnodes,
     1      weights,ifinit,w)
        implicit real *8 (a-h,o-z)
        dimension w(1),vert1(1),
     1      vert3(1),vert2(1),rnodes(1),weights(1)
c
c       this subroutine constructs a tensor product Gaussian
c       quadrature formula on the triangle in the plane
c       specified by its vertices. Note that the one-dimensional
c       Gaussian quadratures are lined up parallel to the side of
c       the triangle opposite the vertex vert1; the total number
c       of nodes and weights to be created is n**2. Actually, this
c       is simply a memory manager for the subroutine triagau0,
c       which does all the work.
c
c                              input parameters:
c
c  n - the number of nodes in each direction. note that the
c       total number of nodes to be created is n**2, and the
c       algebraic order of the quadrature is (n-1)*2
c  vert1,vert2,vert3 - the vertices of the triangle on which the
c       quadrature rule is to be constructed
c
c                              output parameters:
c
c  rnodes - the array of n**2 nodes in the plane (all inside the
c       user-supplied traiangle)
c  weights - the quadrature weights corresponding to the nodes rnodes
c
c                              work arrays:
c
c  w - must be at least 4*n+5 real *8 locations
c
c
c        . . . allocate memory in the work array
c
        it=1
        lt=n+1
c
        iwhts=it+lt
        lwhts=n+2
c
        itvert=iwhts+lwhts
        ltvert=n+2
c
        iwhtsver=itvert+ltvert
        lwhtsver=n+1
c
c       construct the quadrature formula on the triangle
c
        call triagau0(n,vert1,vert2,vert3,rnodes,weights,
     1      ifinit,w(it),w(iwhts),w(itvert),w(iwhtsver) )
c
        return
        end
c
c
c
c
c
        subroutine triagau0(n,vert1,vert2,vert3,rnodes,weights,
     1      ifinit,t,whts,tvert,whtsvert)
        implicit real *8 (a-h,o-z)
        dimension w(12),z1(2),z2(2),z3(2),vert1(1),
     1      vert3(1),vert2(1),rnodes(2,n,n),weights(n,n)
c
        dimension t(1),whts(1),tvert(1),whtsvert(1)
c
c       this subroutine constructs a tensor product Gaussian
c       quadrature formula on the triangle in the plane
c       specified by its vertices. Note that the one-dimensional
c       Gaussian quadratures are lined up parallel to the side of
c       the triangle opposite the vertex vert1; the total number
c       of nodes and weights to be created is n**2.
c
c                              input parameters:
c
c  n - the number of nodes in each direction. note that the
c       total number of nodes to be created is n**2, and the
c       algebraic order of the quadrature is (n-1)*2
c  vert1,vert2,vert3 - the vertices of the triangle on which the
c       quadrature rule is to be constructed
c
c                              output parameters:
c
c  rnodes - the array of n**2 nodes in the plane (all inside the
c       user-supplied traiangle)
c  weights - the quadrature weights corresponding to the nodes rnodes
c
c                              work arrays:
c
c  t,whts,tvert,whtsvert - must be n+1 real *8 locations each
c
c
c       . . . if required, construct the Gaussian nodes and
c             weights on the interval [-1,1]
c
        if(ifinit .eq. 1) then
        ifwhts=1
        call legewhts(n,t,whts,ifwhts)
        endif
c
c        construct affine transformation putting one side
c       of the user-specified triangle on the real axis. The side
c       to be put on the real axis is the one opposite the vertex vert1.
c
        call trianini(vert1,vert2,vert3,w)
c
c       construct the discretization of the version of the triangle
c       moved so that its side is on the real axis
c
        call trianfor(w,vert1,z1)
        call trianfor(w,vert2,z2)
        call trianfor(w,vert3,z3)
c
c       . . . construct the gaussian discretization of the
c             hight of the triangle
c
        b=z1(2)
        u=b/2
        v=b/2
        d=dabs(b/2)
        do 1200 i=1,n
        tvert(i)=u*t(i)+v
        whtsvert(i)=whts(i)*d
 1200 continue
c
c        find the equations of the two non-horizontal sides of
c        the shifted triangle
c
        call trialine(z1(1),z2(1),z1(2),z2(2),a12,b12,c12)
c
        call trialine(z1(1),z3(1),z1(2),z3(2),a13,b13,c13)
c
c        one after another, construct the horizontal discretizations
c
        do 2000 i=1,n
c
c       find the current this horizontal section
c
        y=tvert(i)
        x12=-(b12*y+c12)/a12
        x13=-(b13*y+c13)/a13

cccc        call prin2('x12=*',x12,1)
cccc        call prin2('x13=*',x13,1)
c
        d=dabs(x12-x13)/2
c
        u=(x13-x12)/2
        v=(x13+x12)/2
c
cccc        call prin2('second u=*',u,1)
cccc        call prin2('second v=*',v,1)
c
        do 1400 j=1,n
c
        thor=u*t(j)+v
c
        rnodes(1,j,i)=thor
        rnodes(2,j,i)=y
c
        weights(j,i)=d*whts(j)*whtsvert(i)
c
 1400 continue
 2000 continue
c
c       now, move the nodes in the plane back to the original triangle
c
        do 2400 i=1,n
        do 2200 j=1,n
c
        call trianbak(w,rnodes(1,j,i),z3)
        rnodes(1,j,i)=z3(1)
        rnodes(2,j,i)=z3(2)
 2200 continue
 2400 continue
c
        return
        end
c
c
c
c
c
        subroutine trialine(x1,x2,y1,y2,a,b,c)
        implicit real *8 (a-h,o-z)
c
c        construct the equation of the line passing through the
c        two points (x1,y1), (x2,y2)
c
        dx=x2-x1
        dy=y2-y1
c
        if(dabs(dx) .lt. dabs(dy)) goto 1200
c
        b=1
        a=-dy/dx
        c=-(a*x1+b*y1)
c
        goto 1400
c
 1200 continue
c
        a=1
        b=-dx/dy
        c=-(a*x1+b*y1)
c
 1400 continue
c
c       . . . normalize the coefficients
c
        d=a**2+b**2
        d=dsqrt(d)
c
        a=a/d
        b=b/d
        c=c/d
c
        return
        end
c
c
c
c
c
        subroutine trianini(vert1,vert2,vert3,w)
        implicit real *8 (a-h,o-z)
        dimension vert1(2),vert2(2),vert3(2),w(1),zin(2),zout(2)
        data ixshift/1/,iyshift/2/,iumat/3/,ivmat/7/
c
c       this subroutine constructs an affine transformation putting one side
c       of the user-specified triangle on the real axis. The side
c       to be put on the real axis is the one opposite the vertex vert1.
c       this subroutine also constructs the transformation inverse to the
c       above one. The actual transformations are applied to particular
c       user-specified points by the entries trianfor and trianbak of this
c       subroutine
c
c                              input parameters:
c
c  vert1,vert2,vert3 - the vertices of the triangle one of whose
c       sides is to be put on the x axis
c
c                              output parameters:
c
c  w - the real *8 array of length 10 containing the transformations;
c       it is to be used by the entries trianfor, trianbak (see below)
c
        call triambld(vert1(1),vert1(2),vert2(1),vert2(2),vert3(1),
     1      vert3(2),w(ixshift),w(iyshift),w(iumat),w(ivmat) )
c
        return
c
c
c
c
        entry trianfor(w,zin,zout)
c
c       this entry applies to the user-specified point zin \in R^2
c       the first of the transformations constructed by the entry
c       trianini (see above), obtaining the point zout \in R^2.
c
c
        call triarotf(zin(1),zin(2),w(ixshift),w(iyshift),
     1      w(iumat),zout(1),zout(2) )
c
cccc        call prin2('in trianfor, zout=*',zout,2)

        return
c
c
c
c
        entry trianbak(w,zin,zout)
c
c       this entry applies to the user-specified point zin \in R^2
c       the second of the transformations constructed by the entry
c       trianini (see above), obtaining the point zout \in R^2.
c
        call triarotb(zin(1),zin(2),w(ixshift),w(iyshift),
     1      w(ivmat),zout(1),zout(2) )
c
        return
        end
c
c
c
c
c
        subroutine triarotf(x,y,xshift,yshift,umat,xout,yout)
        implicit real *8 (a-h,o-z)
        dimension umat(2,2)
c
c        apply the shift
c
        xx=x+xshift
        yy=y+yshift

cccc         call prin2('in triarotf, xx=*',xx,1)
cccc         call prin2('in triarotf, yy=*',yy,1)

c
c        apply the rotation
c
        xout=umat(1,1)*xx+umat(1,2)*yy
        yout=umat(2,1)*xx+umat(2,2)*yy
        return
        end
c
c
c
c
c
        subroutine triarotb(x,y,xshift,yshift,vmat,xout,yout)
        implicit real *8 (a-h,o-z)
        dimension vmat(2,2)

c
c        apply the rotation
c
        xx=vmat(1,1)*x+vmat(1,2)*y
        yy=vmat(2,1)*x+vmat(2,2)*y
c
cccc         call prin2('in triarotb, xx=*',xx,1)
cccc         call prin2('in triarotb, yy=*',yy,1)
c
c        apply the shift
c
        xout=xx-xshift
        yout=yy-yshift
        return
        end
c
c
c
c
c
        subroutine triambld(x1,y1,x2,y2,x3,y3,
     1      xshift,yshift,umat,vmat)
        implicit real *8 (a-h,o-z)
        dimension umat(2,2),vmat(2,2)
c
c        construct a motion in the plane that will put the side of
c        the triangle opposite the vertex (x1,y1) on th2 x axis
c
c        . . . the translation
c
        xshift=-x2
        yshift=-y2
c
        xx1=x1+xshift
        xx2=x2+xshift
        xx3=x3+xshift
c
        yy1=y1+yshift
        yy2=y2+yshift
        yy3=y3+yshift
c
c        . . . rotation
c
        u=(x3-x2)
        v=y3-y2
        d=dsqrt(u**2+v**2)
        u=u/d
        v=v/d
c
        umat(1,1)=u
        umat(1,2)=v
        umat(2,1)=-v
        umat(2,2)=u
c
        vmat(1,1)=umat(1,1)
        vmat(2,2)=umat(2,2)
        vmat(1,2)=umat(2,1)
        vmat(2,1)=umat(1,2)
c
        return
        end
c
c
c
c
c
        subroutine triaarrm(x,y,n)
        implicit real *8 (a-h,o-z)
        dimension x(1),y(1)
c
        do 1200 i=1,n
        y(i)=x(i)
 1200 continue
        return
        end
c
c
c
c
c
        subroutine ctriaarrm(x,y,n)
        implicit real *8 (a-h,o-z)
        complex *16 x(1),y(1)
c
        do 1200 i=1,n
        y(i)=x(i)
 1200 continue
        return
        end
c
c
c
c
c
        subroutine triaplot(iw,z)
        implicit real *8 (a-h,o-z)
        dimension z(2,3)
c
 1200 format(2x,e11.5,2x,e11.5)
c
        write(iw,1200) z(1,1),z(2,1)
        write(iw,1200) z(1,2),z(2,2)
        write(iw,1200) z(1,3),z(2,3)
        write(iw,1200) z(1,1),z(2,1)
c
 1400 format(2x)
        write(iw,1400)

        return
        end
c
c
c
c
c
        subroutine triadiv(z1,z2,z3,zout)
        implicit real *8 (a-h,o-z)
        dimension z1(2),z2(2),z3(2),zout(2,3,4),w12(2),
     1      w13(2),w23(2)
c
c        this subroutine subdivides the input triangle
c        with the nodes z1, z2, z3 into four triangles
c
c        . . . construct the middles of the sides
c
        w12(1)=(z1(1)+z2(1))/2
        w12(2)=(z1(2)+z2(2))/2
c
        w13(1)=(z1(1)+z3(1))/2
        w13(2)=(z1(2)+z3(2))/2
c
        w23(1)=(z2(1)+z3(1))/2
        w23(2)=(z2(2)+z3(2))/2
c
c       . . . construct the new triangles
c
        zout(1,1,1)=w12(1)
        zout(2,1,1)=w12(2)
c
        zout(1,2,1)=z1(1)
        zout(2,2,1)=z1(2)
c
        zout(1,3,1)=w13(1)
        zout(2,3,1)=w13(2)
c
c
        zout(1,1,2)=w13(1)
        zout(2,1,2)=w13(2)
c
        zout(1,2,2)=z3(1)
        zout(2,2,2)=z3(2)
c
        zout(1,3,2)=w23(1)
        zout(2,3,2)=w23(2)
c
c
        zout(1,1,3)=w23(1)
        zout(2,1,3)=w23(2)
c
        zout(1,2,3)=z2(1)
        zout(2,2,3)=z2(2)
c
        zout(1,3,3)=w12(1)
        zout(2,3,3)=w12(2)
c
c
        zout(1,1,4)=w23(1)
        zout(2,1,4)=w23(2)
c
        zout(1,2,4)=w12(1)
        zout(2,2,4)=w12(2)
c
        zout(1,3,4)=w13(1)
        zout(2,3,4)=w13(2)
c
        return
        end
c
c
c
c
c


text:
module error
  !! Development tool used to add and print error messages, as well as handle logic flow in the event an error is thrown.

  implicit none

  public :: add_error_message, print_error_list_to_shell, terminate_with_failure

  private

  character(len=1024) :: error_list(1:100)
  integer :: err_i = 0

contains

  subroutine add_error_message(err_s)

    character(len=*), intent(in) :: err_s

    err_i = err_i + 1
    error_list(err_i) = "Error: " // trim(err_s)

  end subroutine add_error_message

  subroutine print_error_list_to_shell()

    integer :: i, l, u

    l = lbound(error_list,1)
    u = ubound(error_list,1)

    do i = l, u
       write(6,*) error_list(i)
    end do

  end subroutine print_error_list_to_shell

  subroutine terminate_with_failure()
    error stop 1
  end subroutine terminate_with_failure

end module error

max_stars_repo_name: ArmchairPhysicist/peano

max_issues_repo_name: ArmchairPhysicist/peano

max_forks_repo_name: ArmchairPhysicist/peano


text:
      subroutine tceigen(ca,cw,ceig,n,ndim)
      implicit none
      integer*4 ,intent(in):: n,ndim
      integer*4 ,parameter ::itmax=30
      integer*4 i,j,k,i1,i2,it
      complex*16 ,intent(inout):: ca(ndim,n)
      complex*16 ,intent(out):: cw(n,n),ceig(n)
      complex*16 cs,cc,ccs,ccc,
     $     cpk,cr,cp,cu,cd,ctr,cdet,csqr,c1,c2
      real*8 a,w,xc,xs,xp
      cw=(0.d0,0.d0)
      do i=1,n
        cw(i,i)=(1.d0,0.d0)
      enddo
      do i=1,n-1
        cp=ca(i+1,i)
        if(imag(cp) .ne. 0.d0)then
          a=abs(cp)
          cu=conjg(cp)/a
          do j=i+1,n
            ca(i+1,j)=ca(i+1,j)*cu
          enddo
          cu=conjg(cu)
          ca(1:n,i+1)=ca(1:n,i+1)*cu
          cw(1:n,i+1)=cw(1:n,i+1)*cu
        else
          a=dble(cp)
        endif
        xp=a
        do j=i+2,n
          a=abs(dcmplx(xp,abs(ca(j,i))))
          if(a .ne. 0.d0)then
            xc=xp/a
            cs=ca(j,i)/a
            ccs=conjg(cs)
            do k=i+1,n
              cpk=ca(i+1,k)
              ca(i+1,k)= cpk*xc +ca(j,k)*ccs
              ca(j  ,k)=-cpk*cs +ca(j,k)*xc
            enddo
            do k=1,n
              cpk=ca(k,i+1)
              ca(k,i+1)= cpk*xc +ca(k,j)*cs
              ca(k,j  )=-cpk*ccs+ca(k,j)*xc
              cpk=cw(k,i+1)
              cw(k,i+1)= cpk*xc +cw(k,j)*cs
              cw(k,j  )=-cpk*ccs+cw(k,j)*xc
            enddo
            xp=a
          endif
        enddo
        ca(i+1,i)=xp
      enddo
      i1=1
      i2=n
 1    if(i2 .eq. 1)then
        go to 3000
      endif
      if(i1 .eq. i2)then
        i2=i1-1
        i1=1
        go to 1
      endif
      it=0
 10   do i=i2-1,i1,-1
        w=abs(ca(i,i))+abs(ca(i+1,i+1))
        if(abs(dble(ca(i+1,i)))+w .eq. w)then
          ca(i+1,i)=0.d0
          if(i .eq. i2-1)then
            i2=i2-1
          else
            i1=i+1
          endif
          go to 1
        endif
      enddo
      ctr=ca(i2-1,i2-1)+ca(i2,i2)
      cdet=2.d0*(ca(i2-1,i2-1)*ca(i2,i2)-
     $     dble(ca(i2,i2-1))*ca(i2-1,i2))
      csqr=sqrt(ctr**2-2.d0*cdet)
      c1=ctr+csqr
      c2=ctr-csqr
      if(abs(c1) .gt. abs(c2))then
        cr=cdet/c1
      else
        cr=cdet/c2
      endif
      cp=ca(i1,i1)-cr
      a=abs(dcmplx(abs(cp),dble(ca(i1+1,i1))))
      cc=cp/a
      xs=dble(ca(i1+1,i1))/a
      ccc=conjg(cc)
      do k=i1,n
        cpk=ca(i1,k)
        ca(i1  ,k)= cpk*ccc+ca(i1+1,k)*xs
        ca(i1+1,k)=-cpk*xs +ca(i1+1,k)*cc
      enddo
      do k=1,i1+1
        cpk=ca(k,i1)
        ca(k,i1  )= cpk*cc +ca(k,i1+1)*xs
        ca(k,i1+1)=-cpk*xs +ca(k,i1+1)*ccc
      enddo
      if(i1+2 .le. i2)then
        ca(i1+2,i1  )= dble(ca(i1+2,i1+1))*xs
        ca(i1+2,i1+1)= dble(ca(i1+2,i1+1))*ccc
      endif
      do k=1,n
        cpk=cw(k,i1)
        cw(k,i1  )= cpk*cc +cw(k,i1+1)*xs
        cw(k,i1+1)=-cpk*xs +cw(k,i1+1)*ccc
      enddo
      do i=i1,i2-2
        cp=ca(i+1,i)
        a=abs(dcmplx(abs(cp),dble(ca(i+2,i))))
        if(a .ne. 0.d0)then
          cc=cp/a
          xs=dble(ca(i+2,i))/a
          ccc=conjg(cc)
          do k=i+1,n
            cpk=ca(i+1,k)
            ca(i+1,k)= cpk*ccc+ca(i+2,k)*xs
            ca(i+2,k)=-cpk*xs +ca(i+2,k)*cc
          enddo
          ca(i+1,i)=a
          do k=1,i+2
            cpk=ca(k,i+1)
            ca(k,i+1)= cpk*cc +ca(k,i+2)*xs
            ca(k,i+2)=-cpk*xs +ca(k,i+2)*ccc
          enddo
          if(i+3 .le. i2)then
            ca(i+3,i+1)= dble(ca(i+3,i+2))*xs
            ca(i+3,i+2)= dble(ca(i+3,i+2))*ccc
          endif
          do k=1,n
            cpk=cw(k,i+1)
            cw(k,i+1)= cpk*cc +cw(k,i+2)*xs
            cw(k,i+2)=-cpk*xs +cw(k,i+2)*ccc
          enddo
        endif
      enddo
      if(imag(ca(i2,i2-1)) .ne. 0.d0)then
        a=abs(ca(i2,i2-1))
        cu=conjg(ca(i2,i2-1))/a
        ca(i2,i2-1)=a
c        do i=i2+1,n
        ca(i2,i2+1:n)=ca(i2,i2+1:n)*cu
c        enddo
        cu=conjg(cu)
c        do i=1,i2-1
        ca(1:i2-1,i2)=ca(1:i2-1,i2)*cu
c        enddo
c        do i=1,n
        cw(:,i2)=cw(:,i2)*cu
c        enddo
      endif
      it=it+1
      if(it .gt. itmax)then
        write(*,*)'TCEIGEN Convergence fail.'
        ca(i2,i2-1)=0.d0
      endif
      go to 10
 3000 do i=2,n
        do j=i-1,1,-1
          cd=ca(j,j)-ca(i,i)
          if(cd .ne. 0.d0)then
            cu=ca(j,i)/cd
c            do k=i,n
            ca(j,i:n)=ca(j,i:n)+cu*ca(i,i:n)
c            enddo
c            do k=1,j-1
            ca(1:j-1,i)=ca(1:j-1,i)-cu*ca(1:j-1,j)
c            enddo
            ca(j,i)=0.d0
c            do k=1,n
            cw(:,i)=cw(:,i)-cu*cw(:,j)
c            enddo
          endif
        enddo
      enddo
      do i=1,n
        ceig(i)=ca(i,i)
      enddo
      ca(1:n,:)=cw
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD