Create a LinearOperator with a function instead of a matrix

[Maxi Miller]

I tried to follow the examples for setting up a LinearOperator (such as shown in Step-22 and Step-20), but instead of a matrix I wanted to provide my own function which should serve as vmult-function (i.e. a matrix-free linear operator). Is that possible? My initial test was the following code:

  class approximator : public Subscriptor
    {
    public:
       approximator(std::function<void(const Vector<double>&,
                                                    Vector<double>&)> residual_function,
                            const Vector<double> &evaluation_point) :
            residual_function(residual_function)
        {
            u = evaluation_point;
            u.add(1);
            epsilon_val = u.l1_norm() * 1e-6;
            residual_function(u, residual_0);
        };
   
        void vmult(Vector<double> &dst, const Vector<double> &src) const;
   
    private:
        std::function<void(const Vector<double>&,
                           Vector<double>&)> residual_function;
        double epsilon_val;
        Vector<double> u, residual_0, residual_1;
    };


const auto op_M = linear_operator(jacobian_approximation(std::bind(&MinimalSurfaceProblem<dim>::compute_residual,
                                                                       this,
                                                                       std::placeholders::_1,
                                                                       std::placeholders::_2),
                                                             present_solution));

but that resulted in

error: use of deleted function ‘dealii::LinearOperator<Range, Domain, Payload> dealii::linear_operator(Matrix&&) [with Range = dealii::Vector<double>; Domain = dealii::Vector<double>; Payload = dealii::internal::LinearOperatorImplementation::EmptyPayload; Matrix = Step15::jacobian_approximation; <template-parameter-1-5> = void]’
  389 |                                                              present_solution));

Are there other ways?

Thanks!

[Matthias Maier]

Hi,

the following example (which is a simplified version of test
lac/linear_operator_01) should get you started. In short: Simply create
an empty LinearOperator object (with the appropriate template
parameters) and populate the corresponding std::function objects.

Best,
Matthias



#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/vector_memory.templates.h>

using namespace dealii;

struct LeftVector {
typedef double value_type;
value_type value;

LeftVector &operator=(value_type new_value)
{
value = new_value;
return *this;
}
LeftVector &operator*=(value_type scale)
{
value *= scale;
return *this;
}
LeftVector &operator/=(value_type scale)
{
value /= scale;
return *this;
}
LeftVector &operator+=(const LeftVector &u)
{
value += u.value;
return *this;
}
int size() const
{
return 1;
}
std::size_t memory_consumption() const
{
return 1;
}
};

struct RightVector {
typedef double value_type;
value_type value;

RightVector &operator=(value_type new_value)
{
value = new_value;
return *this;
}
RightVector &operator*=(value_type scale)
{
value *= scale;
return *this;
}
RightVector &operator/=(value_type scale)
{
value /= scale;
return *this;
}
RightVector &operator+=(const RightVector &u)
{
value += u.value;
return *this;
}
int size() const
{
return 1;
}
std::size_t memory_consumption() const
{
return 1;
}
};

int main()
{
LinearOperator<LeftVector, RightVector> multiply2;

multiply2.vmult = [](LeftVector &v, const RightVector &u) {
v.value = 2 * u.value;
};
multiply2.vmult_add = [](LeftVector &v, const RightVector &u) {
v.value += 2 * u.value;
};
multiply2.Tvmult = [](RightVector &v, const LeftVector &u) {
v.value = 2 * u.value;
};
multiply2.Tvmult_add = [](RightVector &v, const LeftVector &u) {
v.value += 2 * u.value;
};
multiply2.reinit_range_vector = [](LeftVector &, bool omit_zeroing_values) {
// do nothing
};
multiply2.reinit_domain_vector = [](RightVector &, bool omit_zeroing_values) {
// do nothing
};

// Small test:

RightVector u = {4.};
LeftVector v = {0.};

multiply2.vmult(v, u);
std::cout << "2 * " << u.value << " = " << v.value << std::endl;

multiply2.vmult_add(v, u);
std::cout << "... + 2 * " << u.value << " = " << v.value << std::endl;

return 0;
}

[Matthias Maier]

A short comment:

On Tue, Jul 16, 2019, at 11:18 CDT, Matthias Maier <tam...@43-1.org> wrote:

> struct LeftVector {
> };

> struct RightVector {
> };

These two classes are of course just decoration (showing the minimal
interface a vector has to possess). There is usually no need to define
custom Vector clases.

[Maxi Miller]

Is it then possible to use this LinearOperator in a preconditioner? The code works without one (i.e. PreconditionIdentity()), but the Iterations are increasing quite fast, thus negating all speedup gained from the LinearOperator.

I tried to find something in the examples, but only found data about using matrices as input for preconditioners, not pure LinearOperators containing only a vmult-function.

[Wolfgang Bangerth]

On 7/19/19 9:18 AM, 'Maxi Miller' via deal.II User Group wrote:
> Is it then possible to use this LinearOperator in a preconditioner? The code
> works without one (i.e. PreconditionIdentity()), but the Iterations are
> increasing quite fast, thus negating all speedup gained from the LinearOperator.
> I tried to find something in the examples, but only found data about using
> matrices as input for preconditioners, not pure LinearOperators containing
> only a vmult-function.

Yes, this is supposed to work!
Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Maxi Miller]

I.e. (based on step-37) I can write something like

LinearOperator jacobian_operator;
//vmult-operation is added later

 using SystemMatrixType = jacobian_operator;
 SystemMatrixType system_matrix;

and follow the example accordingly?

[Wolfgang Bangerth]

On 7/19/19 9:34 AM, 'Maxi Miller' via deal.II User Group wrote:
> I.e. (based on step-37) I can write something like
> |
> LinearOperator jacobian_operator;
> //vmult-operation is added later
>
> using SystemMatrixType = jacobian_operator;
> SystemMatrixType system_matrix;
> |
>
> and follow the example accordingly?

No. `jacobian_operator` here is the name of a local variable, but in a `using`
statement, the right hand side needs to be a *type*.

What I'm saying is that you should be able to write something of the sort
SolverCG<...> solver_cg;
solver_cg.solve (system_matrix, solution, rhs,
jacobian_operator);
where the last argument is any kind of type that provides a vmult something --
whether it's a matrix, a class with a vmult function, or a linear operator.

[Maxi Miller]

What I did was the replacement of system_matrix with a function, providing a vmult()-function. Now I would like to add a preconditioner based on that function to reduce the necessary GMRES-iterations. But until now most of the preconditioners require a matrix for initialization. Thus I wanted to replace that matrix with my LinearOperator, if that is possible?

[Wolfgang Bangerth]

Yes, as I already said, this is possible. Have you tried?

I think you should also be able to look at step-20 (and maybe 22) for examples.

[Maxi Miller]

Would that work together with the GMG-preconditioner? Will test it there as soon as my problem in my second question could be solved.

There I have the following code:

   

SolverControl        coarse_solver_control(1000, 1e-10, false, false);
    SolverGMRES<vector_t>   coarse_solver(coarse_solver_control);
    PreconditionIdentity id;
    MGCoarseGridIterativeSolver<vector_t,
            SolverGMRES<vector_t>,
            matrix_t,
            PreconditionIdentity>
            coarse_grid_solver(coarse_solver, coarse_matrix, id);

    using Smoother = LA::MPI::PreconditionJacobi;
    MGSmootherPrecondition<matrix_t, Smoother, vector_t> mg_smoother;
    mg_smoother.initialize(mg_matrices, Smoother::AdditionalData(0.5));
    mg_smoother.set_steps(2);

    mg::Matrix<vector_t> mg_matrix(mg_matrices);
    mg::Matrix<vector_t> mg_interface_up(mg_interface_matrices);
    mg::Matrix<vector_t> mg_interface_down(mg_interface_matrices);

    Multigrid<vector_t> mg(mg_matrix, coarse_grid_solver, mg_transfer, mg_smoother, mg_smoother);
    mg.set_edge_matrices(mg_interface_down, mg_interface_up);

    PreconditionMG<dim, vector_t, MGTransferPrebuilt<vector_t>> preconditioner(dof_handler,
                                                                               mg,
                                                                               mg_transfer);

and here I intent to replace all occurences of matrix_t with my custom LinearOperator (if that makes sense).

[Wolfgang Bangerth]

> Would that work together with the GMG-preconditioner?

You mean you want to use the linear operator as a smoother for the GMG? Or
where in the GMG do you want to use it?

> Will test it there as
> soon as my problem in my second question could be solved.
> There I have the following code:
> |
> SolverControl       coarse_solver_control(1000,1e-10,false,false);
> SolverGMRES<vector_t>  coarse_solver(coarse_solver_control);
> PreconditionIdentityid;
> MGCoarseGridIterativeSolver<vector_t,
> SolverGMRES<vector_t>,
>             matrix_t,
> PreconditionIdentity>
>             coarse_grid_solver(coarse_solver,coarse_matrix,id);
>
> usingSmoother=LA::MPI::PreconditionJacobi;
> MGSmootherPrecondition<matrix_t,Smoother,vector_t>mg_smoother;

>     mg_smoother.initialize(mg_matrices,Smoother::AdditionalData(0.5));

>     mg_smoother.set_steps(2);
>
>     mg::Matrix<vector_t>mg_matrix(mg_matrices);
>     mg::Matrix<vector_t>mg_interface_up(mg_interface_matrices);
>     mg::Matrix<vector_t>mg_interface_down(mg_interface_matrices);
>
> Multigrid<vector_t>mg(mg_matrix,coarse_grid_solver,mg_transfer,mg_smoother,mg_smoother);
>     mg.set_edge_matrices(mg_interface_down,mg_interface_up);
>

> PreconditionMG<dim,vector_t,MGTransferPrebuilt<vector_t>>preconditioner(dof_handler,

>
>  mg,
>
>  mg_transfer);
> |
>
> and here I intent to replace all occurences of matrix_t with my custom
> LinearOperator (if that makes sense).

I don't know the multigrid well enough, but I would expect that it really only
wants something that has a vmult.

My recommendation would be, if you don't know whether something would work, to
just try it out.