Multiple Jacobian calls with kinsol and Picard iterations

[mottelet]

Hello,

I noticed that using kinsol Picard iterations triggers a Jacobian evaluation at every iteration, although the Jacobian is supposed to be constant.

What am I doing wrong in the following trivial example (solve x*x-2=0) ? Here is the output of the program, showing that the jacobian function is called at each step though it shouldn't (IMHO):

jac jac jac jac jac jac jac jac jac jac jac jac jac jac y=1.414213

Thanks for your help

Stéphane Mottelet

#include <stdio.h>
#include <kinsol/kinsol.h>             /* access to KINSOL func., consts. */
#include <nvector/nvector_serial.h>    /* access to serial N_Vector       */
#include <sunlinsol/sunlinsol_dense.h>  /* access to band SUNLinearSolver  */
#include <sunmatrix/sunmatrix_dense.h>  /* access to dense SUNMatrix        */
#include <sundials/sundials_types.h>   /* defs. of realtype, sunindextype */

int func(N_Vector U, N_Vector F, void *user_data)
{
  double u = NV_Ith_S(U,0);
  NV_Ith_S(F,0) = u*u-2;
  return(0);
}
int jac(N_Vector u, N_Vector f, SUNMatrix J,
               void *user_data, N_Vector tmp1, N_Vector tmp2)
{
  SM_ELEMENT_D(J,0,0) = 2;
  printf("jac\n");
  return(0);
}

int main()
{
  N_Vector y, scale;
  void *kmem;
  SUNMatrix J;
  SUNLinearSolver LS;
  int flag;

  y = N_VNew_Serial(1);
  scale = N_VNew_Serial(1);
  N_VConst(1.0, y);
  N_VConst(1.0,scale);

  kmem = KINCreate();
  flag = KINInit(kmem, func, y);

  J = SUNDenseMatrix(1, 1);
  LS = SUNLinSol_Dense(y, J);
  KINSetLinearSolver(kmem, LS, J);
  KINSetJacFn(kmem, jac);

  flag = KINSol(kmem,y,KIN_PICARD,scale,scale);

  printf("flag=%d,y=%f\n",flag,NV_Ith_S(y,0));
}

[jandrej]

If your function is f(x) = x^2-2 then df(x)/dx = 2.0*x. Your Jacobian returns df(x)/dx = 2.0.

[mottelet]

If your function is f(x) = x^2-2 then df(x)/dx = 2.0*x. Your Jacobian returns df(x)/dx = 2.0.

I know how to compute a Jacobian, please see the Kinsol documentation about Picard method were the rhs is spliited in a linear part and a nonlinear part (https://sundials.readthedocs.io/en/latest/kinsol/Mathematics_link.html?highlight=picard#basic-fixed-point-iteration) :

For Picard iteration, as implemented in kinsol, we consider a special form of the nonlinear function F , such that F (u) = Lu − N (u), where L is a constant nonsingular matrix and N is (in general) nonlinear. Then the fixed-point function G is defined as G(u) = u − L−1F (u). The Picard iteration is given by:...,

[gardner48]

Currently the Picard iteration does not utilize the same linear system reuse controls that are present when using the modified Newton iteration. As such, the Picard iteration evaluates the linear system each iteration. We'll look at incorporating the reuse logic into the Picard iteration as part of a future release.

[mottelet]

OK thanks, so let us the issue open ok ?

[gardner48]

Yes, let's keep this issue open.