fv_two_phase_hess_fun.m

function [res, precond] = fv_two_phase_hess_fun(params, jac, HessTemp, arg)

  Plength = params.Slength;
  if isempty(HessTemp)
    HessTemp = jac(1:end-1,:)'*jac(1:end-1,:);
  end
  if isempty(arg)
    res = HessTemp;
  else
    Pscale = jac(end,end-Plength+1:end);
    PscaleP = Pscale(:) .* arg(end-Plength+1:end);
    res = HessTemp * arg;
    res(end-Plength+1:end) = res(end-Plength+1:end) + Pscale(:) .* PscaleP;

    resref = (jac'*jac) * arg;
    if norm(res - resref) > 1e-12
        keyboard;
    end
  end

  if nargout == 2
    [precond.L, precond.U] = ilu(HessTemp, struct('type', 'crout'));
  end

end