[petsc-users] Using matrix-free method for SNES problems with an explicit Jacobian

[Justin Chang]

Okay that makes sense. Reason for asking this is because I am solving the
following convex optimization problem for a linear diffusion problem:

min 1/2 x^T * A * x + x^T * g
s.t. 0 <= x <= 1

Using Tao's BLMVM solver, the objective function and gradient vector reads
as follows:

f = 0.5*x^T*A*x + x^T*(r - A*x_0)
g = A*(x - x_0) + r

Where A, x_0, and r are the jacobian, initial guess, and residual
respectively. To compute these, I have TaoSetObjectiveAndGradientRoutine()
invoke the following function:

#undef __FUNCT__
#define __FUNCT__ "FormFunctionGradient"
PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g,
void *ctx)
{
  AppCtx *user = (AppCtx*)ctx;
  PetscScalar       xtHx;
  PetscErrorCode    ierr;

  PetscFunctionBegin;
  ierr = MatMult(user->A,x,g);CHKERRQ(ierr);
  ierr = VecDot(x,g,&xtHx);CHKERRQ(ierr);
  ierr = VecDot(x,user->f,f);CHKERRQ(ierr);
  *f += 0.5*xtHx;
  ierr = VecAXPY(g,1.0,user->f);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

where user->f = (r - A*x_0)  The timing results from -log_summary seem to
indicate that MatMult(user->A,x,g) is a bottleneck within TaoSolve() so I
was wondering if I could somehow make the computation of Vec g matrix-free
considering that user->A is constant. The problem proliferates if I decide
to use TRON. Or is this unavoidable (Jason any thoughts?)

Thanks,
Justin

On Wed, May 13, 2015 at 3:15 PM, Matthew Knepley <knepley at gmail.com> wrote:

> On Wed, May 13, 2015 at 2:54 PM, Justin Chang <jychang48 at gmail.com> wrote:
>
>> Hello everyone,
>>
>> From what I read online, one could use a matrix-free method if using the
>> conjugate gradient method for a symmetric and positive definite problem.
>> Some of the PETSc SNES examples like ex12 compute the Jacobian A
>> explicitly, but I was wondering if it's possible to convert this
>> problem into matrix free form. That is, to directly assemble the product of
>> A*x into vector form without first forming A and then invoking MatMult().
>>
>
> Yes, Jed advocates this all the time, but there is some infrastructure
> that you want for this
> to be really effective. For quadratic and higher orders, it makes sense to
> compute the action
> of your operator MF. This is done by defining something that looks a lot
> like your residual
> evaluation (element vec in and out) but computes the linearization.
> However, you still usually
> need a low order explicit matrix for preconditioning. I have not spent
> time hooking that all together.
>
>   Thanks,
>
>     Matt
>
>
>> Thanks,
>> Justin
>>
>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150513/521ef37d/attachment.html>