htswanalysis/MACS/lib/gsl/gsl-1.11/multiroots/test_funcs.c

   1 /* multiroots/test_funcs.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <math.h>
  22 #include <gsl/gsl_vector.h>
  23 #include <gsl/gsl_matrix.h>
  24 #include <gsl/gsl_multiroots.h>
  25
  26 #include "test_funcs.h"
  27
  28 /* For information on testing see the following paper,
  29
  30    J.J More, B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained
  31    Optimization Software", ACM Transactions on Mathematical Software,
  32    Vol 7, No 1, (1981) p 17-41
  33
  34    */
  35
  36 /* Rosenbrock Function */
  37
  38 gsl_multiroot_function_fdf rosenbrock =
  39 {&rosenbrock_f,
  40  &rosenbrock_df,
  41  &rosenbrock_fdf,
  42  2, 0};
  43
  44 void
  45 rosenbrock_initpt (gsl_vector * x)
  46 {
  47   gsl_vector_set (x, 0, -1.2);
  48   gsl_vector_set (x, 1, 1.0);
  49 }
  50
  51 int
  52 rosenbrock_f (const gsl_vector * x, void *params, gsl_vector * f)
  53 {
  54   double x0 = gsl_vector_get (x, 0);
  55   double x1 = gsl_vector_get (x, 1);
  56
  57   double y0 = 1 - x0;
  58   double y1 = 10 * (x1 - x0 * x0);
  59
  60   gsl_vector_set (f, 0, y0);
  61   gsl_vector_set (f, 1, y1);
  62
  63   params = 0;                   /* avoid warning about unused parameters */
  64
  65   return GSL_SUCCESS;
  66 }
  67
  68 int
  69 rosenbrock_df (const gsl_vector * x, void *params, gsl_matrix * df)
  70 {
  71   double x0 = gsl_vector_get (x, 0);
  72
  73   double df00 = -1;
  74   double df01 = 0;
  75   double df10 = -20 * x0;
  76   double df11 = 10;
  77
  78   gsl_matrix_set (df, 0, 0, df00);
  79   gsl_matrix_set (df, 0, 1, df01);
  80   gsl_matrix_set (df, 1, 0, df10);
  81   gsl_matrix_set (df, 1, 1, df11);
  82
  83   params = 0;                   /* avoid warning about unused parameters */
  84
  85   return GSL_SUCCESS;
  86 }
  87
  88 int
  89 rosenbrock_fdf (const gsl_vector * x, void *params,
  90                 gsl_vector * f, gsl_matrix * df)
  91 {
  92   rosenbrock_f (x, params, f);
  93   rosenbrock_df (x, params, df);
  94
  95   return GSL_SUCCESS;
  96 }
  97
  98
  99 /* Freudenstein and Roth function */
 100
 101 gsl_multiroot_function_fdf roth =
 102 {&roth_f,
 103  &roth_df,
 104  &roth_fdf,
 105  2, 0};
 106
 107 void
 108 roth_initpt (gsl_vector * x)
 109 {
 110   gsl_vector_set (x, 0, 4.5);  /* changed from the value in the paper */
 111   gsl_vector_set (x, 1, 3.5);  /* otherwise the problem is too hard */
 112 }
 113
 114 int
 115 roth_f (const gsl_vector * x, void *params, gsl_vector * f)
 116 {
 117   double x0 = gsl_vector_get (x, 0);
 118   double x1 = gsl_vector_get (x, 1);
 119
 120   double y0 = -13.0 + x0 + ((5.0 - x1)*x1 - 2.0)*x1;
 121   double y1 = -29.0 + x0 + ((x1 + 1.0)*x1 - 14.0)*x1;
 122
 123   gsl_vector_set (f, 0, y0);
 124   gsl_vector_set (f, 1, y1);
 125
 126   params = 0;                   /* avoid warning about unused parameters */
 127
 128   return GSL_SUCCESS;
 129 }
 130
 131 int
 132 roth_df (const gsl_vector * x, void *params, gsl_matrix * df)
 133 {
 134   double x1 = gsl_vector_get (x, 1);
 135
 136   double df00 = 1;
 137   double df01 = -3 * x1 * x1 + 10 * x1 - 2;
 138   double df10 = 1;
 139   double df11 = 3 * x1 * x1 + 2 * x1 - 14;
 140
 141   gsl_matrix_set (df, 0, 0, df00);
 142   gsl_matrix_set (df, 0, 1, df01);
 143   gsl_matrix_set (df, 1, 0, df10);
 144   gsl_matrix_set (df, 1, 1, df11);
 145
 146   params = 0;                   /* avoid warning about unused parameters */
 147
 148   return GSL_SUCCESS;
 149 }
 150
 151 int
 152 roth_fdf (const gsl_vector * x, void *params,
 153                 gsl_vector * f, gsl_matrix * df)
 154 {
 155   roth_f (x, params, f);
 156   roth_df (x, params, df);
 157
 158   return GSL_SUCCESS;
 159 }
 160
 161
 162
 163 /* Powell badly scaled function */
 164
 165 gsl_multiroot_function_fdf powellscal =
 166 {&powellscal_f,
 167  &powellscal_df,
 168  &powellscal_fdf,
 169  2, 0};
 170
 171 void
 172 powellscal_initpt (gsl_vector * x)
 173 {
 174   gsl_vector_set (x, 0, 0.0);
 175   gsl_vector_set (x, 1, 1.0);
 176 }
 177
 178 int
 179 powellscal_f (const gsl_vector * x, void *params, gsl_vector * f)
 180 {
 181   double x0 = gsl_vector_get (x, 0);
 182   double x1 = gsl_vector_get (x, 1);
 183
 184   double y0 = 10000.0 * x0 * x1 - 1.0;
 185   double y1 = exp (-x0) + exp (-x1) - 1.0001;
 186
 187   gsl_vector_set (f, 0, y0);
 188   gsl_vector_set (f, 1, y1);
 189
 190   params = 0;                   /* avoid warning about unused parameters */
 191
 192   return GSL_SUCCESS;
 193 }
 194
 195 int
 196 powellscal_df (const gsl_vector * x, void *params, gsl_matrix * df)
 197 {
 198   double x0 = gsl_vector_get (x, 0);
 199   double x1 = gsl_vector_get (x, 1);
 200
 201   double df00 = 10000.0 * x1, df01 = 10000.0 * x0;
 202   double df10 = -exp (-x0), df11 = -exp (-x1);
 203
 204   gsl_matrix_set (df, 0, 0, df00);
 205   gsl_matrix_set (df, 0, 1, df01);
 206   gsl_matrix_set (df, 1, 0, df10);
 207   gsl_matrix_set (df, 1, 1, df11);
 208
 209   params = 0;                   /* avoid warning about unused parameters */
 210
 211   return GSL_SUCCESS;
 212 }
 213
 214 int
 215 powellscal_fdf (const gsl_vector * x, void *params,
 216                   gsl_vector * f, gsl_matrix * df)
 217 {
 218   powellscal_f (x, params, f);
 219   powellscal_df (x, params, df);
 220
 221   return GSL_SUCCESS;
 222 }
 223
 224
 225 /* Brown badly scaled function */
 226
 227 gsl_multiroot_function_fdf brownscal =
 228 {&brownscal_f,
 229  &brownscal_df,
 230  &brownscal_fdf,
 231  2, 0};
 232
 233 void
 234 brownscal_initpt (gsl_vector * x)
 235 {
 236   gsl_vector_set (x, 0, 1.0);
 237   gsl_vector_set (x, 1, 1.0);
 238 }
 239
 240 int
 241 brownscal_f (const gsl_vector * x, void *params, gsl_vector * f)
 242 {
 243   double x0 = gsl_vector_get (x, 0);
 244   double x1 = gsl_vector_get (x, 1);
 245
 246   double y0 = x0 - 1e6;
 247   double y1 = x0 * x1 - 2;
 248
 249   gsl_vector_set (f, 0, y0);
 250   gsl_vector_set (f, 1, y1);
 251
 252   params = 0;                   /* avoid warning about unused parameters */
 253
 254   return GSL_SUCCESS;
 255 }
 256
 257 int
 258 brownscal_df (const gsl_vector * x, void *params, gsl_matrix * df)
 259 {
 260   double x0 = gsl_vector_get (x, 0);
 261   double x1 = gsl_vector_get (x, 1);
 262
 263   double df00 = 1.0, df01 = 0.0;
 264   double df10 = x1, df11 = x0;
 265
 266   gsl_matrix_set (df, 0, 0, df00);
 267   gsl_matrix_set (df, 0, 1, df01);
 268   gsl_matrix_set (df, 1, 0, df10);
 269   gsl_matrix_set (df, 1, 1, df11);
 270
 271   params = 0;                   /* avoid warning about unused parameters */
 272
 273   return GSL_SUCCESS;
 274 }
 275
 276 int
 277 brownscal_fdf (const gsl_vector * x, void *params,
 278                   gsl_vector * f, gsl_matrix * df)
 279 {
 280   brownscal_f (x, params, f);
 281   brownscal_df (x, params, df);
 282
 283   return GSL_SUCCESS;
 284 }
 285
 286
 287 /* Powell Singular Function */
 288
 289 gsl_multiroot_function_fdf powellsing =
 290 {&powellsing_f,
 291  &powellsing_df,
 292  &powellsing_fdf,
 293  4, 0};
 294
 295 void
 296 powellsing_initpt (gsl_vector * x)
 297 {
 298   gsl_vector_set (x, 0, 3.0);
 299   gsl_vector_set (x, 1, -1.0);
 300   gsl_vector_set (x, 2, 0.0);
 301   gsl_vector_set (x, 3, 1.0);
 302 }
 303
 304 int
 305 powellsing_f (const gsl_vector * x, void *params, gsl_vector * f)
 306 {
 307   double x0 = gsl_vector_get (x, 0);
 308   double x1 = gsl_vector_get (x, 1);
 309   double x2 = gsl_vector_get (x, 2);
 310   double x3 = gsl_vector_get (x, 3);
 311
 312   double y0 = x0 + 10 * x1;
 313   double y1 = sqrt (5.0) * (x2 - x3);
 314   double y2 = pow (x1 - 2 * x2, 2.0);
 315   double y3 = sqrt (10.0) * pow (x0 - x3, 2.0);
 316
 317   gsl_vector_set (f, 0, y0);
 318   gsl_vector_set (f, 1, y1);
 319   gsl_vector_set (f, 2, y2);
 320   gsl_vector_set (f, 3, y3);
 321
 322   params = 0;                   /* avoid warning about unused parameters */
 323
 324   return GSL_SUCCESS;
 325 }
 326
 327 int
 328 powellsing_df (const gsl_vector * x, void *params, gsl_matrix * df)
 329 {
 330   double x0 = gsl_vector_get (x, 0);
 331   double x1 = gsl_vector_get (x, 1);
 332   double x2 = gsl_vector_get (x, 2);
 333   double x3 = gsl_vector_get (x, 3);
 334
 335   double df00 = 1, df01 = 10, df02 = 0, df03 = 0;
 336   double df10 = 0, df11 = 0, df12 = sqrt (5.0), df13 = -df12;
 337   double df20 = 0, df21 = 2 * (x1 - 2 * x2), df22 = -2 * df21, df23 = 0;
 338   double df30 = 2 * sqrt (10.0) * (x0 - x3), df31 = 0, df32 = 0, df33 = -df30;
 339
 340   gsl_matrix_set (df, 0, 0, df00);
 341   gsl_matrix_set (df, 0, 1, df01);
 342   gsl_matrix_set (df, 0, 2, df02);
 343   gsl_matrix_set (df, 0, 3, df03);
 344
 345   gsl_matrix_set (df, 1, 0, df10);
 346   gsl_matrix_set (df, 1, 1, df11);
 347   gsl_matrix_set (df, 1, 2, df12);
 348   gsl_matrix_set (df, 1, 3, df13);
 349
 350   gsl_matrix_set (df, 2, 0, df20);
 351   gsl_matrix_set (df, 2, 1, df21);
 352   gsl_matrix_set (df, 2, 2, df22);
 353   gsl_matrix_set (df, 2, 3, df23);
 354
 355   gsl_matrix_set (df, 3, 0, df30);
 356   gsl_matrix_set (df, 3, 1, df31);
 357   gsl_matrix_set (df, 3, 2, df32);
 358   gsl_matrix_set (df, 3, 3, df33);
 359
 360   params = 0;                   /* avoid warning about unused parameters */
 361
 362   return GSL_SUCCESS;
 363 }
 364
 365 int
 366 powellsing_fdf (const gsl_vector * x, void *params,
 367                     gsl_vector * f, gsl_matrix * df)
 368 {
 369   powellsing_f (x, params, f);
 370   powellsing_df (x, params, df);
 371
 372   return GSL_SUCCESS;
 373 }
 374
 375
 376 /* Wood function */
 377
 378 gsl_multiroot_function_fdf wood =
 379 {&wood_f,
 380  &wood_df,
 381  &wood_fdf,
 382  4, 0};
 383
 384 void
 385 wood_initpt (gsl_vector * x)
 386 {
 387   gsl_vector_set (x, 0, -3.0);
 388   gsl_vector_set (x, 1, -1.0);
 389   gsl_vector_set (x, 2, -3.0);
 390   gsl_vector_set (x, 3, -1.0);
 391 }
 392
 393 int
 394 wood_f (const gsl_vector * x, void *params, gsl_vector * f)
 395 {
 396   double x0 = gsl_vector_get (x, 0);
 397   double x1 = gsl_vector_get (x, 1);
 398   double x2 = gsl_vector_get (x, 2);
 399   double x3 = gsl_vector_get (x, 3);
 400
 401   double t1 = x1 - x0 * x0;
 402   double t2 = x3 - x2 * x2;
 403
 404   double y0 = -200.0 * x0 * t1 - (1 - x0);
 405   double y1 = 200.0 * t1 + 20.2 * (x1 - 1) + 19.8 * (x3 - 1);
 406   double y2 = -180.0 * x2 * t2 - (1 - x2);
 407   double y3 = 180.0 * t2 + 20.2 * (x3 - 1) + 19.8 * (x1 - 1);
 408
 409   gsl_vector_set (f, 0, y0);
 410   gsl_vector_set (f, 1, y1);
 411   gsl_vector_set (f, 2, y2);
 412   gsl_vector_set (f, 3, y3);
 413
 414   params = 0;                   /* avoid warning about unused parameters */
 415
 416   return GSL_SUCCESS;
 417 }
 418
 419 int
 420 wood_df (const gsl_vector * x, void *params, gsl_matrix * df)
 421 {
 422   double x0 = gsl_vector_get (x, 0);
 423   double x1 = gsl_vector_get (x, 1);
 424   double x2 = gsl_vector_get (x, 2);
 425   double x3 = gsl_vector_get (x, 3);
 426
 427   double t1 = x1 - 3 * x0 * x0;
 428   double t2 = x3 - 3 * x2 * x2;
 429
 430   double df00 = -200.0 * t1 + 1, df01 = -200.0 * x0, df02 = 0, df03 = 0;
 431   double df10 = -400.0*x0, df11 = 200.0 + 20.2, df12 = 0, df13 = 19.8;
 432   double df20 = 0, df21 = 0, df22 = -180.0 * t2 + 1, df23 = -180.0 * x2;
 433   double df30 = 0, df31 = 19.8, df32 = -2 * 180 * x2, df33 = 180.0 + 20.2;
 434
 435   gsl_matrix_set (df, 0, 0, df00);
 436   gsl_matrix_set (df, 0, 1, df01);
 437   gsl_matrix_set (df, 0, 2, df02);
 438   gsl_matrix_set (df, 0, 3, df03);
 439
 440   gsl_matrix_set (df, 1, 0, df10);
 441   gsl_matrix_set (df, 1, 1, df11);
 442   gsl_matrix_set (df, 1, 2, df12);
 443   gsl_matrix_set (df, 1, 3, df13);
 444
 445   gsl_matrix_set (df, 2, 0, df20);
 446   gsl_matrix_set (df, 2, 1, df21);
 447   gsl_matrix_set (df, 2, 2, df22);
 448   gsl_matrix_set (df, 2, 3, df23);
 449
 450   gsl_matrix_set (df, 3, 0, df30);
 451   gsl_matrix_set (df, 3, 1, df31);
 452   gsl_matrix_set (df, 3, 2, df32);
 453   gsl_matrix_set (df, 3, 3, df33);
 454
 455   params = 0;                   /* avoid warning about unused parameters */
 456
 457   return GSL_SUCCESS;
 458 }
 459
 460 int
 461 wood_fdf (const gsl_vector * x, void *params,
 462                     gsl_vector * f, gsl_matrix * df)
 463 {
 464   wood_f (x, params, f);
 465   wood_df (x, params, df);
 466
 467   return GSL_SUCCESS;
 468 }
 469
 470
 471 /* Helical Valley Function */
 472
 473 gsl_multiroot_function_fdf helical =
 474 {&helical_f,
 475  &helical_df,
 476  &helical_fdf,
 477  3, 0};
 478
 479 void
 480 helical_initpt (gsl_vector * x)
 481 {
 482   gsl_vector_set (x, 0, -1.0);
 483   gsl_vector_set (x, 1, 0.0);
 484   gsl_vector_set (x, 2, 0.0);
 485 }
 486
 487 int
 488 helical_f (const gsl_vector * x, void *params, gsl_vector * f)
 489 {
 490   double x0 = gsl_vector_get (x, 0);
 491   double x1 = gsl_vector_get (x, 1);
 492   double x2 = gsl_vector_get (x, 2);
 493
 494   double t1, t2;
 495   double y0, y1, y2;
 496
 497   if (x0 > 0)
 498     {
 499       t1 = atan(x1/x0) / (2.0 * M_PI);
 500     }
 501   else if (x0 < 0)
 502     {
 503       t1 = 0.5 + atan(x1/x0) / (2.0 * M_PI);
 504     }
 505   else
 506     {
 507       t1 = 0.25 * (x1 > 0 ? +1 : -1);
 508     }
 509
 510   t2 = sqrt(x0*x0 + x1*x1) ;
 511
 512   y0 = 10 * (x2 - 10 * t1);
 513   y1 = 10 * (t2 - 1);
 514   y2 = x2 ;
 515
 516   gsl_vector_set (f, 0, y0);
 517   gsl_vector_set (f, 1, y1);
 518   gsl_vector_set (f, 2, y2);
 519
 520   params = 0;                   /* avoid warning about unused parameters */
 521
 522   return GSL_SUCCESS;
 523 }
 524
 525 int
 526 helical_df (const gsl_vector * x, void *params, gsl_matrix * df)
 527 {
 528   double x0 = gsl_vector_get (x, 0);
 529   double x1 = gsl_vector_get (x, 1);
 530
 531   double t = x0 * x0 + x1 * x1 ;
 532   double t1 = 2 * M_PI * t ;
 533   double t2 = sqrt(t) ;
 534
 535   double df00 = 100*x1/t1, df01 = -100.0 * x0/t1, df02 = 10.0;
 536   double df10 = 10*x0/t2, df11 = 10*x1/t2, df12 = 0;
 537   double df20 = 0, df21 = 0, df22 = 1.0;
 538
 539   gsl_matrix_set (df, 0, 0, df00);
 540   gsl_matrix_set (df, 0, 1, df01);
 541   gsl_matrix_set (df, 0, 2, df02);
 542
 543   gsl_matrix_set (df, 1, 0, df10);
 544   gsl_matrix_set (df, 1, 1, df11);
 545   gsl_matrix_set (df, 1, 2, df12);
 546
 547   gsl_matrix_set (df, 2, 0, df20);
 548   gsl_matrix_set (df, 2, 1, df21);
 549   gsl_matrix_set (df, 2, 2, df22);
 550
 551   params = 0;                   /* avoid warning about unused parameters */
 552
 553   return GSL_SUCCESS;
 554 }
 555
 556 int
 557 helical_fdf (const gsl_vector * x, void *params,
 558                     gsl_vector * f, gsl_matrix * df)
 559 {
 560   helical_f (x, params, f);
 561   helical_df (x, params, df);
 562
 563   return GSL_SUCCESS;
 564 }
 565
 566
 567 /* Discrete Boundary Value Function */
 568
 569 #define N 10
 570
 571 gsl_multiroot_function_fdf dbv =
 572 {&dbv_f,
 573  &dbv_df,
 574  &dbv_fdf,
 575  N, 0};
 576
 577 void
 578 dbv_initpt (gsl_vector * x)
 579 {
 580   size_t i;
 581   double h = 1.0 / (N + 1.0);
 582
 583   for (i = 0; i < N; i++)
 584     {
 585       double t = (i + 1) * h;
 586       double z = t * (t - 1);
 587       gsl_vector_set (x, i, z);
 588     }
 589 }
 590
 591 int
 592 dbv_f (const gsl_vector * x, void *params, gsl_vector * f)
 593 {
 594   size_t i;
 595
 596   double h = 1.0 / (N + 1.0);
 597
 598   for (i = 0; i < N; i++)
 599     {
 600       double z, ti = (i + 1) * h;
 601       double xi = 0, xim1 = 0, xip1 = 0;
 602
 603       xi = gsl_vector_get (x, i);
 604
 605       if (i > 0)
 606         xim1 = gsl_vector_get (x, i - 1);
 607
 608       if (i < N - 1)
 609         xip1 = gsl_vector_get (x, i + 1);
 610
 611       z = 2 * xi - xim1 - xip1 + h * h * pow(xi + ti + 1, 3.0) / 2.0;
 612
 613       gsl_vector_set (f, i, z);
 614
 615     }
 616
 617   params = 0;                   /* avoid warning about unused parameters */
 618
 619   return GSL_SUCCESS;
 620 }
 621
 622 int
 623 dbv_df (const gsl_vector * x, void *params, gsl_matrix * df)
 624 {
 625   size_t i, j;
 626
 627   double h = 1.0 / (N + 1.0);
 628
 629   for (i = 0; i < N; i++)
 630     for (j = 0; j < N; j++)
 631       gsl_matrix_set (df, i, j, 0.0);
 632
 633   for (i = 0; i < N; i++)
 634     {
 635       double dz_dxi, ti = (i + 1) * h;
 636
 637       double xi = gsl_vector_get (x, i);
 638
 639       dz_dxi = 2.0 + (3.0 / 2.0) * h * h * pow(xi + ti + 1, 2.0) ;
 640
 641       gsl_matrix_set (df, i, i, dz_dxi);
 642
 643       if (i > 0)
 644         gsl_matrix_set (df, i, i-1, -1.0);
 645
 646       if (i < N - 1)
 647         gsl_matrix_set (df, i, i+1, -1.0);
 648
 649     }
 650
 651   params = 0;                   /* avoid warning about unused parameters */
 652
 653   return GSL_SUCCESS;
 654 }
 655
 656 int
 657 dbv_fdf (const gsl_vector * x, void *params,
 658                     gsl_vector * f, gsl_matrix * df)
 659 {
 660   dbv_f (x, params, f);
 661   dbv_df (x, params, df);
 662
 663   return GSL_SUCCESS;
 664 }
 665
 666 /* Trigonometric Function */
 667
 668 gsl_multiroot_function_fdf trig =
 669 {&trig_f,
 670  &trig_df,
 671  &trig_fdf,
 672  N, 0};
 673
 674 void
 675 trig_initpt (gsl_vector * x)
 676 {
 677   size_t i;
 678
 679   for (i = 0; i < N; i++)       /* choose an initial point which converges */
 680     {
 681       gsl_vector_set (x, i, 0.05);
 682     }
 683 }
 684
 685 int
 686 trig_f (const gsl_vector * x, void *params, gsl_vector * f)
 687 {
 688   size_t i;
 689   double sum = 0;
 690
 691   for (i = 0; i < N; i++)
 692     {
 693       sum += cos(gsl_vector_get(x,i));
 694     }
 695
 696   for (i = 0; i < N; i++)
 697     {
 698       double xi = gsl_vector_get (x,i);
 699       double z = N - sum + (i + 1) * (1 - cos(xi)) - sin(xi);
 700
 701       gsl_vector_set (f, i, z);
 702     }
 703
 704   params = 0;                   /* avoid warning about unused parameters */
 705
 706   return GSL_SUCCESS;
 707 }
 708
 709 int
 710 trig_df (const gsl_vector * x, void *params, gsl_matrix * df)
 711 {
 712   size_t i, j;
 713
 714   for (i = 0; i < N; i++)
 715     {
 716       for (j = 0; j < N; j++)
 717         {
 718           double dz;
 719           double xi = gsl_vector_get(x, i);
 720           double xj = gsl_vector_get(x, j);
 721
 722           if (j == i)
 723             dz = sin(xi) + (i + 1) * sin(xi) - cos(xi);
 724           else
 725             dz = sin(xj);
 726
 727           gsl_matrix_set(df, i, j, dz);
 728         }
 729     }
 730
 731   params = 0;                   /* avoid warning about unused parameters */
 732
 733   return GSL_SUCCESS;
 734 }
 735
 736 int
 737 trig_fdf (const gsl_vector * x, void *params,
 738                     gsl_vector * f, gsl_matrix * df)
 739 {
 740   trig_f (x, params, f);
 741   trig_df (x, params, df);
 742
 743   return GSL_SUCCESS;
 744 }