gsl_multifit__lmset.c 1.3 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758
  1. static int
  2. set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale)
  3. {
  4. lmder_state_t *state = (lmder_state_t *) vstate;
  5. gsl_matrix *r = state->r;
  6. gsl_vector *tau = state->tau;
  7. gsl_vector *diag = state->diag;
  8. gsl_vector *work1 = state->work1;
  9. gsl_permutation *perm = state->perm;
  10. int signum;
  11. GSL_MULTIFIT_FN_EVAL_F_DF (fdf, x, f, J);
  12. state->par = 0;
  13. state->iter = 1;
  14. state->fnorm = enorm (f);
  15. gsl_vector_set_all (dx, 0.0);
  16. /* store column norms in diag */
  17. if (scale)
  18. {
  19. compute_diag (J, diag);
  20. }
  21. else
  22. {
  23. gsl_vector_set_all (diag, 1.0);
  24. }
  25. /* set delta to 100 |D x| or to 100 if |D x| is zero */
  26. state->xnorm = scaled_enorm (diag, x);
  27. state->delta = compute_delta (diag, x);
  28. /* Factorize J into QR decomposition */
  29. gsl_matrix_memcpy (r, J);
  30. gsl_linalg_QRPT_decomp (r, tau, perm, &signum, work1);
  31. gsl_vector_set_zero (state->rptdx);
  32. gsl_vector_set_zero (state->w);
  33. /* Zero the trial vector, as in the alloc function */
  34. gsl_vector_set_zero (state->f_trial);
  35. #ifdef DEBUG
  36. printf("r = "); gsl_matrix_fprintf(stdout, r, "%g");
  37. printf("perm = "); gsl_permutation_fprintf(stdout, perm, "%d");
  38. printf("tau = "); gsl_vector_fprintf(stdout, tau, "%g");
  39. #endif
  40. return GSL_SUCCESS;
  41. }