VMult-function in the MatrixFree-Framework

[Maxi Miller]

In my project I already managed to replace the system_matrix in the solve()-call by a custom LinearOperator, but still I had to retain the matrix itself (for the preconditioner, for example). Thus I am now trying to use the MatrixFree-framework for that.

Thus, some questions arose:

• When using the LinearOperator, the function vmult() is used when solving the system with a GMRES-solver (or similar). Does the same hold up for the MatrixFree-framework?
• If the above is correct: In each GMRES-iteration I have to calculate the residual of my equation, once for the current solution and once for the current solution + a small step. At the moment I calculate the residual by summing the terms into my residual vector, and then distributing the values while setting the boundary condition to 0. If I understand it correctly, that happens as soon as I call phi.distribute_local_to_global(dst), with phi a FEEvaluation-element. This means that I first have to assemble the residual vector, then distribute it, and then can use the vmult-function. Is that correct, or do I have to use another approach?

Thanks!

[Daniel Arndt]

• When using the LinearOperator, the function vmult() is used when solving the system with a GMRES-solver (or similar). Does the same hold up for the MatrixFree-framework?

Yes.

• If the above is correct: In each GMRES-iteration I have to calculate the residual of my equation, once for the current solution and once for the current solution + a small step. At the moment I calculate the residual by summing the terms into my residual vector, and then distributing the values while setting the boundary condition to 0. If I understand it correctly, that happens as soon as I call phi.distribute_local_to_global(dst), with phi a FEEvaluation-element.

Yes.

• This means that I first have to assemble the residual vector, then distribute it, and then can use the vmult-function. Is that correct, or do I have to use another approach?

Isn't assembling and distributing already what `vmult` is supposed to do? Can you elaborate on what you are doing in each GMRES iteration?

Best,

Daniel

[Maxi Miller]

My aim is to implement the jacobian approximation, i.e. (if F(u) = 0 is the function to solve) I wanted to implement

dst = (F(current_solution + epsilon * src) - F(current_solution))/epsilon.

In my LinearOperator, my code is

LinearOperator<LinearAlgebraTrilinos::MPI::Vector> jacobian_operator;
jacobian_operator.vmult = [&](LinearAlgebraTrilinos::MPI::Vector &dst, const LinearAlgebraTrilinos::MPI::Vector &src){
    local_src = src;
    extended_src = src;
    local_solution = present_solution;

    double epsilon = 1e-6;
    if(local_src.l2_norm() != 0){
        epsilon = sqrt((1 + local_solution.l2_norm()) * 1e-16) / local_src.l2_norm();
    };
    forward_eps_solution = present_solution;
    backward_eps_solution = present_solution;

    extended_src *= epsilon;
    forward_eps_solution += extended_src;
    backward_eps_solution -= extended_src;

    compute_residual(backward_eps_solution, cur_residual);
    compute_residual(forward_eps_solution, eps_residual);
    eps_residual -= cur_residual;
    eps_residual /= (2 * epsilon);
    dst.reinit(locally_owned_dofs,
   mpi_communicator);
    dst = eps_residual;
};

Is the same/similar possible for the matrix-free-approach?

Thanks!

[Daniel Arndt]

My aim is to implement the jacobian approximation, i.e. (if F(u) = 0 is the function to solve) I wanted to implement

dst = (F(current_solution + epsilon * src) - F(current_solution))/epsilon.

In my LinearOperator, my code is

LinearOperator<LinearAlgebraTrilinos::MPI::Vector> jacobian_operator;
jacobian_operator.vmult = [&](LinearAlgebraTrilinos::MPI::Vector &dst, const LinearAlgebraTrilinos::MPI::Vector &src){
    local_src = src;
    extended_src = src;
    local_solution = present_solution;

    double epsilon = 1e-6;
    if(local_src.l2_norm() != 0){
        epsilon = sqrt((1 + local_solution.l2_norm()) * 1e-16) / local_src.l2_norm();
    };
    forward_eps_solution = present_solution;
    backward_eps_solution = present_solution;

    extended_src *= epsilon;
    forward_eps_solution += extended_src;
    backward_eps_solution -= extended_src;

    compute_residual(backward_eps_solution, cur_residual);
    compute_residual(forward_eps_solution, eps_residual);
    eps_residual -= cur_residual;
    eps_residual /= (2 * epsilon);
    dst.reinit(locally_owned_dofs,
   mpi_communicator);
    dst = eps_residual;
};

Is the same/similar possible for the matrix-free-approach?

Sure, you can always write a wrapper around the actual vmult call (which I assume is happening in compute_residual()).

Best,

Daniel

[Maxi Miller]

No, that code is the vmult-call I am using in my LinearOperator (and is used directly in my solve()-call). The calculate_residual()-function is similar to the function in step-15, with the difference, that I do not provide alpha, but the vector, and return the calculated residual value. Thus, as far as I could understand, I have to do the following in the matrix-free cell loop:

• Calculate residual of the current solution
• Calculate residual of the current solution + epsilon * src
• Calculate dst, and return it

On the other hand, would it be sufficient to unroll the residual-calculation, i.e. if the residual is F(u) = nabla^2 u + f, to write (in the cell-loop): (nabla^2(u + eps*src) + f - nabla^2u - f)/epsilon? Or would that lead to wrong results?

[Daniel Arndt]

Maxi,

No, that code is the vmult-call I am using in my LinearOperator (and is used directly in my solve()-call). The calculate_residual()-function is similar to the function in step-15, with the difference, that I do not provide alpha, but the vector, and return the calculated residual value.

 

So this is where your cell loop currently lives and that can be represented by a matrix-vector product. :-)

 

Thus, as far as I could understand, I have to do the following in the matrix-free cell loop:

• Calculate residual of the current solution
• Calculate residual of the current solution + epsilon * src
• Calculate dst, and return it

On the other hand, would it be sufficient to unroll the residual-calculation, i.e. if the residual is F(u) = nabla^2 u + f, to write (in the cell-loop): (nabla^2(u + eps*src) + f - nabla^2u - f)/epsilon? Or would that lead to wrong results?

 

It is sufficient to only change compute_diagonal() to use the MatrixFree framework. Passing multiple vectors at once might be a little bit more difficult to implement.

Currently, you residual seems to depend linearly on u. In that case, the difference does not depend on u at all.

Best,

Daniel

[Maxi Miller]

Am Mittwoch, 24. Juli 2019 00:39:14 UTC+2 schrieb Daniel Arndt:

Maxi,

No, that code is the vmult-call I am using in my LinearOperator (and is used directly in my solve()-call). The calculate_residual()-function is similar to the function in step-15, with the difference, that I do not provide alpha, but the vector, and return the calculated residual value.

 

So this is where your cell loop currently lives and that can be represented by a matrix-vector product. :-)

 

Thus, as far as I could understand, I have to do the following in the matrix-free cell loop:

• Calculate residual of the current solution
• Calculate residual of the current solution + epsilon * src
• Calculate dst, and return it

On the other hand, would it be sufficient to unroll the residual-calculation, i.e. if the residual is F(u) = nabla^2 u + f, to write (in the cell-loop): (nabla^2(u + eps*src) + f - nabla^2u - f)/epsilon? Or would that lead to wrong results?

 

It is sufficient to only change compute_diagonal() to use the MatrixFree framework. Passing multiple vectors at once might be a little bit more difficult to implement.

Could you elaborate on that further (changing the diagonal value)? 

Currently, you residual seems to depend linearly on u. In that case, the difference does not depend on u at all.

I know, the equation is linear anyway, but was intended to be used as simple example. For further tests the non-linear equations will be used.

Best,

Daniel

[Daniel Arndt]

On the other hand, would it be sufficient to unroll the residual-calculation, i.e. if the residual is F(u) = nabla^2 u + f, to write (in the cell-loop): (nabla^2(u + eps*src) + f - nabla^2u - f)/epsilon? Or would that lead to wrong results?

 

It is sufficient to only change compute_diagonal() to use the MatrixFree framework. Passing multiple vectors at once might be a little bit more difficult to implement.

Could you elaborate on that further (changing the diagonal value)? 

Oh, that was a typo. I meant that you only have to implement calculate_residual and can leave the rest of you LinearOperator as is.

Best,

Daniel

[Maxi Miller]

Is there then a way to use LinearOperator as Input for a preconditioner? Else I still have to form the full system matrix (which I would like to avoid).

[Daniel Arndt]

Is there then a way to use LinearOperator as Input for a preconditioner? Else I still have to form the full system matrix (which I would like to avoid).

Most precoditioners require access to individual matrx entries. If the class you derive from linear operator provides a suitable interface, you can use it in a preconditioner.

Building suitable preconditioners compatible with matrix-free methods is non-trivial. The default choice would be PreconditionChebyshev.

Best,

Daniel

[Maxi Miller]

My LinearOperator only provides a vmult-interface, nothing else, but for initializing the preconditioner I still need a sparse matrix, thus I can not use the LinearOperator, as far as I understand. Or is there a way to use it?

[Wolfgang Bangerth]

On 7/24/19 12:23 PM, 'Maxi Miller' via deal.II User Group wrote:
> My LinearOperator only provides a vmult-interface, nothing else, but for
> initializing the preconditioner I still need a sparse matrix, thus I can not
> use the LinearOperator, as far as I understand. Or is there a way to use it?

You just found the problem everyone faces. The LinearOperator interface is
useful if you only know how to express a linear operator through its action,
without knowing the entries of the matrix that represent this action. This is
all you need for the iterative solvers we have.

But to work well, these iterators all need preconditioners, and most of the
preconditioners we have require knowledge of the entries of the matrix --
i.e., there is no longer any advantage to using LinearOperator over building a
matrix. The ways you can work around this are Chebyshev iterations, for
example; step-20 also in great detail explores this issue. But the point is
that it's a difficult issue.

Best
W.


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

[Maxi Miller]

Thus I intended to use a GMG-preconditioner, with a Chebyshev-smoother. Usually the GMG-preconditioner does not require the main matrix either (yes, small matrices are still required, but they do not require the same amount of storage space), but here I run into the problem which is mentioned in a separate question, i.e. I can not get the gradients/values of my function on non-active cells while having to loop over both non-active and active cells for building the MG-matrices, if that makes sense?