aprod_ez Subroutine

    subroutine aprod_ez ( me, mode, m, n, x, y )

    implicit none

    class(lsqr_solver_ez),intent(inout) :: me
    integer,intent(in) :: mode  !! * If `mode = 1`, compute `y = y + A*x`.
                                !!   `y` should be altered without changing x.
                                !! * If `mode = 2`, compute `x = x + A(transpose)*y`.
                                !!   `x` should be altered without changing `y`.
    integer,intent(in) :: m     !! number of rows in `A` matrix
    integer,intent(in) :: n     !! number of columns in `A` matrix
    real(wp),dimension(:),intent(inout) :: x  !! [n]
    real(wp),dimension(:),intent(inout) :: y  !! [m]

    integer :: i    !! counter
    integer :: r    !! row index
    integer :: c    !! column index

    if (m/=me%m .or. n/=me%n) error stop 'lsqr_solver_ez class not properly initialized'

    select case (mode)

    case(1)    ! y = y + A*x

        !   A    x   Ax
        !  ---   -   -
        !  X0X   X   X
        !  0X0 * X = X
        !  00X   X   X
        !  00X       X

        ! A*x:
        if (.not. allocated(me%Ax)) allocate(me%Ax(me%m))
        me%Ax = zero
        do i = 1, me%num_nonzero_elements
            r = me%irow(i)
            c = me%icol(i)
            me%Ax(r) = me%Ax(r) + me%a(i)*x(c)
        end do

        y = y + me%Ax

    case(2)   ! x = x + A(transpose)*y

        !   A     Y   ATy
        !  ---    -   -
        !  X000   Y   X
        !  0X00 * Y = X
        !  X0XX   Y   X
        !         Y

        ! A(transpose)*y
        if (.not. allocated(me%Aty)) allocate(me%Aty(me%n))
        me%Aty = zero
        do i = 1, me%num_nonzero_elements
            r = me%irow(i)
            c = me%icol(i)
            me%Aty(c) = me%Aty(c) + me%a(i)*y(r)
        end do

        x = x + me%Aty

    case default
       error stop 'invalid mode in aprod_ez'
    end select

    end subroutine aprod_ez