[Bug-gsl] bug in compute_gradient_direction?

[Yajun Wang]

Hi buddies:

I want to use the Nonlinear Least Squares Fitting. I am working on the
sum of 3-dimensional formulas. However, the program quits quietly
whenever I have only 2 formulas to solve.

I looked deep into the code, and find it quits in function:
compute_gradient_direction (lmpar.c):

  const size_t n = r->size2;
  size_t i, j;
  for (j = 0; j < n; j++)
    {
     double sum = 0;
      for (i = 0; i <= j; i++)
        {
          sum += gsl_matrix_get (r, i, j) * gsl_vector_get (qtf, i);
        }

      {
        size_t pj = gsl_permutation_get (p, j);
        double dpj = gsl_vector_get (diag, pj);

        gsl_vector_set (g, j, sum / dpj);
      }
    }

During my debug, r has size (2,3), thus r->size1 == 2. size2 ==3. And
qtf is of size 2. (Sorry, I have no clue what are those matrices and
vectors, I get the size from gdb).

However, in the loops, the matrix r is assumed to be square(same row
with column, or more row than column at least). But it is not the case
during my program.

I am not sure how this happen. Is this code expected to work on the
square matrix? If not, how could it be fixed? If yes, what's wrong
with my following setting?

   gsl_multifit_function_fdf  f; //fdf functions.
    gsl_vector_view x_x = gsl_vector_view_array( x_init, p);

    f.fdf = &(J_fdf);
    f.f = &(J_f);
    f.df = &(J_df);
    f.n = N; //error happens with N==2.
    f.p = p; //p is set to 3 all the time.
    f.params = &params;

    T = gsl_multifit_fdfsolver_lmsder;
    s = gsl_multifit_fdfsolver_alloc (T, N, p);
    gsl_multifit_fdfsolver_set(s, &f, &x_x.vector);
    int status;
    size_t iter = 0;
    if(N == 2)
        status = 0; //stop here, error.
    do
    {
        iter++;
        status = gsl_multifit_fdfsolver_iterate(s);
        if(status)
            break;
        status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
    }
    while(status == GSL_CONTINUE && iter < 100); //we at most iterate 100 times.
    gsl_multifit_fdfsolver_free(s);

Thank you very much.

regards,
yalding