# matrix multiply and its adjoint
#
subroutine matmult( adj, bb, x,nx, y,ny)
integer ix, iy, adj, nx, ny
real bb(ny,nx), x(nx), y(ny)
if( adj == 0 )
do iy= 1, ny
y(iy) = 0.
else
do ix= 1, nx
x(ix) = 0.
do ix= 1, nx {
do iy= 1, ny {
if( adj == 0 )
y(iy) = y(iy) + bb(iy,ix) * x(ix)
else
x(ix) = x(ix) + bb(iy,ix) * y(iy)
}}
return; end