Hi, I posted a question about troubles I originally had on http://forum.kde.org/viewtopic.php?f=74&t=96197 .

HybridNonLinear solver with numerical evaluation of the Jacobian uses Givens rotations to update the Jacobian. Givens rotations are evaluated in r1updt.h, which receives std::vector<JacobiRotation<Scalar> > &v_givens, &w_givens; these vectors have the size of number of unknowns (n), but their values are unitialized. Loops in r1updt update elements of v_givens resp. w_givens, but only if corresponding value v[i] resp. w[i] is non-zero. If the value is zero, the element of v_givens resp. w_givens is not updated (keeping unitialized garbage in it):

 for (j=n-2; j>=0; --j) {
    w[j] = 0.;
    if (v[j] != 0.) {
       /* ... */
       v_givens[j]=givens;
       /* ... */
    }
    /* BUG: what if v[j]==0 ?!! */
 }

*All* values of v_givens resp. w_givens are used in r1mpyq (without testing whether they are valid or not), effectively propagating garbage into the Jacobian.

My solution was to assign zero rotation to those elements, i.e.

 for (j=n-2; j>=0; --j) {
    w[j] = 0.;
    if (v[j] != 0.) {
       /* ... */
    }
    else v_givens[j]=JacobiRotation<Scalar>(1,0);  // <-- added
 }
 
I would like to have someone knowing details of the algorithm to say whether such solution is correct (it works for me).