lin_op.hpp

//          Copyright Gunter Winkler 2004 - 2007.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

// authors: Gunter Winkler <guwi17 at gmx dot de>
//          Konstantin Kutzkow


#ifndef __INCLUDE_LIN_OP_HPP__
#define __INCLUDE_LIN_OP_HPP__


template < class MATRIX >
class LinOp
{
public:
    template < class VEC >
    VEC operator()(const VEC & x) const {
        VEC y( image_size() );
        mult(x,y);
        return y;
    }

    template < class VEC1, class VEC2 >
    void mult(const VEC1 & x, VEC2 & y) const {
        y = prod(A, x);
    }
    
    template < class VEC1, class VEC2, class VEC3 >
    void residuum(const VEC1 & x, const VEC2 & b, VEC3 & y) const {
  //orig:      y = b - prod(A, x);
        y = b; 
        axpy_prod(A,-x,y,false);
    }
    
    std::size_t image_size() const {
        return A.size1();
    }
    
    std::size_t preimage_size() const {
        return A.size2();
    }

    LinOp(MATRIX const & matrix) : A(matrix) { }
        
private:
    MATRIX const & A;
};


#endif