msci Code

#include <stdio.h>

#include <math.h>

#include "gsl/gsl_linalg.h"
#include "gsl/gsl_poly.h"

double nearest_of_three(double value,
			double x1,
			double x2,
			double x3) {

  double diff1, diff2, diff3;
  double result;

  diff1=fabs(value-x1);
  diff2=fabs(value-x2);
  diff3=fabs(value-x3);

  if (diff1 > diff2) {
    if (diff2 > diff3) result = x3; else result = x2;
  } else {
    if (diff1 > diff3) result = x3; else result = x1;
  }

  return result;

}

//minimization function which brackets the root
//and then approximates root by fitting a third-
//order polynomial

double supernewton(void (*funcd) (double, void *, double *, double *),
		void *params,	//fucntion parameters
		double x1,	//first bracket
		double x2,	//second bracket
		double xtol,	//desired tolerance in x direction
		double ytol,	//desired tolerance in y direction
		long maxiter,	//maximum number of iterations
		long &err,	//error code--positive for success
		double &y1,	//to avoid re-calculating
		double &dydx1,	//these items
		double y2,
		double dydx2)
{

  //function evaluations at the brackets:
  /*
  double y1;
  double dydx1;
  double y2;
  double dydx2;
  */

  //polynomial coefficients:
  double a, b, c, d;

  double x0;		//solution

  double y0;		//solution value for y
  double dydx0;		//solution value for dydx

  //intermediates in calculation:
  double x1_2, x1_3, x2_2, x2_3;
  long nroot;

  //roots of polynomial:
  double x0_1, x0_2, x0_3;

  //number of iterations:
  long i;

  //x and y error:
  double xerr, yerr;

  //approximate slope of the function:
  double m;

  gsl_vector *yp;
  gsl_matrix *A;

  int gsl_status;	//error state from GSL calls

  yp=gsl_vector_alloc(4);
  A=gsl_matrix_alloc(4, 4);

  err=0;
  i=0;
  do {
/*
    printf("Brackets: [%f, %f]\n", x1, x2);
    printf("          [%f, %f]\n", y1, y2);
    printf("          [%f, %f]\n", dydx1, dydx2);
*/

    i++;

    //populate the matrix and solution vector:
    x1_2=x1*x1;
    x1_3=x1_2*x1;
    x2_2=x2*x2;
    x2_3=x2_2*x2;

    gsl_matrix_set(A, 0, 0, 1);
    gsl_matrix_set(A, 0, 1, x1);
    gsl_matrix_set(A, 0, 2, x1_2);
    gsl_matrix_set(A, 0, 3, x1_3);

    gsl_vector_set(yp, 0, y1);

    gsl_matrix_set(A, 1, 0, 1);
    gsl_matrix_set(A, 1, 1, x2);
    gsl_matrix_set(A, 1, 2, x2_2);
    gsl_matrix_set(A, 1, 3, x2_3);

    gsl_vector_set(yp, 1, y2);

    gsl_matrix_set(A, 2, 0, 0);
    gsl_matrix_set(A, 2, 1, 1);
    gsl_matrix_set(A, 2, 2, 2*x1);
    gsl_matrix_set(A, 2, 3, 3*x1_2);

    gsl_vector_set(yp, 2, dydx1);

    gsl_matrix_set(A, 3, 0, 0);
    gsl_matrix_set(A, 3, 1, 1);
    gsl_matrix_set(A, 3, 2, 2*x2);
    gsl_matrix_set(A, 3, 3, 3*x2_2);

    gsl_vector_set(yp, 3, dydx2);

    //solve using Householder transformations:
    gsl_status=gsl_linalg_HH_svx(A, yp);

    if (gsl_status != 0) {
      fprintf(stderr, "Supernewton: gsl_linalg_HH_svx returned the following error:\n");
      fprintf(stderr, "\"%s\" (%d)\n", gsl_strerror(gsl_status), gsl_status);
      err=-i;
      break;
    }

    a=gsl_vector_get(yp, 0);
    b=gsl_vector_get(yp, 1);
    c=gsl_vector_get(yp, 2);
    d=gsl_vector_get(yp, 3);

    //printf("Polynomial coeffs.:%f, %f, %f, %f\n", a, b, c, d);

    //solve the cubic:
    nroot=gsl_poly_solve_cubic(c/d, b/d, a/d, &x0_1, &x0_2, &x0_3);

    if (nroot==1) {
      //if there is only one root, we use that one:
      x0=x0_1;
      //printf("Root: %f\n", x0);
    } else {
      //if there are three,
      //take the root that's closest to the centre:
      //printf("Roots: %f, %f, %f\n", x0_1, x0_2, x0_3);
      x0=(x1+x2)/2;
      x0=nearest_of_three(x0, x0_1, x0_2, x0_3);
    }

    //evaluate the function at the approximated root
    (*funcd) (x0, params, &y0, &dydx0);

    m=(y2-y1)/(x2-x1);
    yerr=fabs(y0);
    //approximate the x error from the y error and approx. slope of func.:
    //xerr=fabs(y0/m);
    xerr=fabs(x1-x2);

    //test for convergence:
    //(*note* this convergence test is designed for speed, NOT accuracy)
    if (xerr < xtol || yerr < ytol) {
      err=i;
      break;
    }

    if (i > maxiter) {
      fprintf(stderr, "Maximum number of iterations exceeded in supernewton\n");
      err = -maxiter;
      break;
      printf("Brackets: [%f, %f]\n", x1, x2);
      printf("Polynomial coeffs.:%f, %f, %f, %f\n", a, b, c, d);
      printf("Roots: %f, %f, %f\n", x0_1, x0_2, x0_3);
    }

    //rebracket the "true" root:
    if (y0*y1 > 0) {
      x1=x0;
      y1=y0;
      dydx1=dydx0;
    } else {
      x2=x0;
      y2=y0;
      dydx2=dydx0;
    }

  } while(1);

  //place the function values in the fourth and third last parameters:
  y1=y0;
  dydx1=dydx0;

  //printf("Supernewton: %d iterations required to reach convergence\n", i);

  //clean up:
  gsl_matrix_free(A);
  gsl_vector_free(yp);

  return x0;

}