# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "<stdin>"
module matmult
! matrix multiply and its adjoint
  use adj_mod
  implicit none
  real ,private, dimension (:,:), pointer :: bb
  contains
  subroutine matmult_init ( bb_in )
    real , dimension (:,:), pointer :: bb_in
    bb => bb_in
  end subroutine
  function matmult_lop ( adj, add, x, y) result(stat)
    integer :: stat
    logical,intent(in) :: adj,add
    real,dimension(:) :: x,y
    call adjnull (adj,add,x,y )
    call matmult_lop2(adj,add,x,y )
    stat=0
  end function
  subroutine matmult_lop2(adj,add,x,y)
    logical,intent(in) :: adj,add
    real, dimension (:) :: x
    real, dimension (:) :: y
    integer ix, iy
    do ix= 1, size(x)
      do iy= 1, size(y)
        if ( adj) then
          x(ix) = x(ix) + bb(iy,ix) * y(iy)
        else
          y(iy) = y(iy) + bb(iy,ix) * x(ix)
        end if
      end do
    end do
  end subroutine
  subroutine matmult_close()
  end subroutine
end module