#$=head1 NAME
#$
#$matmult - matrix multiplication
#$
#$=head1 SYNOPSIS
#$
#$C<call matmult_init(bb)>
#$
#$C<ierr=matmult_lop(adj,add,model,data)>
#$
#$=head1 PARAMETERS
#$
#$=over 4
#$
#$=item  bb  C<real(:,:)>  
#$
#$       Matrix
#$
#$=item adj,add,model,data  
#$
#$      Standard operator interface
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$Matrix multiplication
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
#$
module matmult {                              # matrix multiply and its adjoint
real, dimension (:,:), pointer :: bb
#%  _init( bb)
#%  _lop(  x, y) 
integer ix, iy
do ix= 1, size(x) {
do iy= 1, size(y) {
        if( adj)
		x(ix) = x(ix) + bb(iy,ix) * y(iy)
        else	
		y(iy) = y(iy) + bb(iy,ix) * x(ix)
        }}
}