test_hahn1.c

开放gsl矩阵运算

const size_t hahn1_N = 236;
const size_t hahn1_P = 7;

/* double hahn1_x0[7] = { 10, -1, 0.05, -0.00001, -0.05, 0.001, -0.000001 }; */

double hahn1_x0[7] = { 1, -0.1, 0.005, -0.000001, -0.005, 0.0001, -0.0000001}; 

double hahn1_x[7] = {
1.0776351733E+00,
-1.2269296921E-01,
4.0863750610E-03,
-1.4262662514E-06,
-5.7609940901E-03,
2.4053735503E-04,
-1.2314450199E-07
};

double hahn1_sumsq = 1.5324382854E+00;

double hahn1_sigma[7] = {
  1.7070154742E-01,
  1.2000289189E-02,
  2.2508314937E-04,
  2.7578037666E-07,
  2.4712888219E-04,
  1.0449373768E-05,
  1.3027335327E-08
};

double hahn1_F1[236] = {
        .591E0,
       1.547E0,
       2.902E0,
       2.894E0,
       4.703E0,
       6.307E0,
       7.03E0 ,
       7.898E0,
       9.470E0,
       9.484E0,
      10.072E0,
      10.163E0,
      11.615E0,
      12.005E0,
      12.478E0,
      12.982E0,
      12.970E0,
      13.926E0,
      14.452E0,
      14.404E0,
      15.190E0,
      15.550E0,
      15.528E0,
      15.499E0,
      16.131E0,
      16.438E0,
      16.387E0,
      16.549E0,
      16.872E0,
      16.830E0,
      16.926E0,
      16.907E0,
      16.966E0,
      17.060E0,
      17.122E0,
      17.311E0,
      17.355E0,
      17.668E0,
      17.767E0,
      17.803E0,
      17.765E0,
      17.768E0,
      17.736E0,
      17.858E0,
      17.877E0,
      17.912E0,
      18.046E0,
      18.085E0,
      18.291E0,
      18.357E0,
      18.426E0,
      18.584E0,
      18.610E0,
      18.870E0,
      18.795E0,
      19.111E0,
        .367E0,
        .796E0,
       0.892E0,
       1.903E0,
       2.150E0,
       3.697E0,
       5.870E0,
       6.421E0,
       7.422E0,
       9.944E0,
      11.023E0,
      11.87E0 ,
      12.786E0,
      14.067E0,
      13.974E0,
      14.462E0,
      14.464E0,
      15.381E0,
      15.483E0,
      15.59E0 ,
      16.075E0,
      16.347E0,
      16.181E0,
      16.915E0,
      17.003E0,
      16.978E0,
      17.756E0,
      17.808E0,
      17.868E0,
      18.481E0,
      18.486E0,
      19.090E0,
      16.062E0,
      16.337E0,
      16.345E0,
      16.388E0,
      17.159E0,
      17.116E0,
      17.164E0,
      17.123E0,
      17.979E0,
      17.974E0,
      18.007E0,
      17.993E0,
      18.523E0,
      18.669E0,
      18.617E0,
      19.371E0,
      19.330E0,
       0.080E0,
       0.248E0,
       1.089E0,
       1.418E0,
       2.278E0,
       3.624E0,
       4.574E0,
       5.556E0,
       7.267E0,
       7.695E0,
       9.136E0,
       9.959E0,
       9.957E0,
      11.600E0,
      13.138E0,
      13.564E0,
      13.871E0,
      13.994E0,
      14.947E0,
      15.473E0,
      15.379E0,
      15.455E0,
      15.908E0,
      16.114E0,
      17.071E0,
      17.135E0,
      17.282E0,
      17.368E0,
      17.483E0,
      17.764E0,
      18.185E0,
      18.271E0,
      18.236E0,
      18.237E0,
      18.523E0,
      18.627E0,
      18.665E0,
      19.086E0,
       0.214E0,
       0.943E0,
       1.429E0,
       2.241E0,
       2.951E0,
       3.782E0,
       4.757E0,
       5.602E0,
       7.169E0,
       8.920E0,
      10.055E0,
      12.035E0,
      12.861E0,
      13.436E0,
      14.167E0,
      14.755E0,
      15.168E0,
      15.651E0,
      15.746E0,
      16.216E0,
      16.445E0,
      16.965E0,
      17.121E0,
      17.206E0,
      17.250E0,
      17.339E0,
      17.793E0,
      18.123E0,
      18.49E0 ,
      18.566E0,
      18.645E0,
      18.706E0,
      18.924E0,
      19.1E0  ,
       0.375E0,
       0.471E0,
       1.504E0,
       2.204E0,
       2.813E0,
       4.765E0,
       9.835E0,
      10.040E0,
      11.946E0,
      12.596E0,
      13.303E0,
      13.922E0,
      14.440E0,
      14.951E0,
      15.627E0,
      15.639E0,
      15.814E0,
      16.315E0,
      16.334E0,
      16.430E0,
      16.423E0,
      17.024E0,
      17.009E0,
      17.165E0,
      17.134E0,
      17.349E0,
      17.576E0,
      17.848E0,
      18.090E0,
      18.276E0,
      18.404E0,
      18.519E0,
      19.133E0,
      19.074E0,
      19.239E0,
      19.280E0,
      19.101E0,
      19.398E0,
      19.252E0,
      19.89E0 ,
      20.007E0,
      19.929E0,
      19.268E0,
      19.324E0,
      20.049E0,
      20.107E0,
      20.062E0,
      20.065E0,
      19.286E0,
      19.972E0,
      20.088E0,
      20.743E0,
      20.83E0 ,
      20.935E0,
      21.035E0,
      20.93E0 ,
      21.074E0,
      21.085E0,
      20.935E0
};


double hahn1_F0[236] = {
    24.41E0,
    34.82E0,
    44.09E0,
    45.07E0,
    54.98E0,
    65.51E0,
    70.53E0,
    75.70E0,
    89.57E0,
    91.14E0,
    96.40E0,
    97.19E0,
   114.26E0,
   120.25E0,
   127.08E0,
   133.55E0,
   133.61E0,
   158.67E0,
   172.74E0,
   171.31E0,
   202.14E0,
   220.55E0,
   221.05E0,
   221.39E0,
   250.99E0,
   268.99E0,
   271.80E0,
   271.97E0,
   321.31E0,
   321.69E0,
   330.14E0,
   333.03E0,
   333.47E0,
   340.77E0,
   345.65E0,
   373.11E0,
   373.79E0,
   411.82E0,
   419.51E0,
   421.59E0,
   422.02E0,
   422.47E0,
   422.61E0,
   441.75E0,
   447.41E0,
   448.7E0 ,
   472.89E0,
   476.69E0,
   522.47E0,
   522.62E0,
   524.43E0,
   546.75E0,
   549.53E0,
   575.29E0,
   576.00E0,
   625.55E0,
    20.15E0,
    28.78E0,
    29.57E0,
    37.41E0,
    39.12E0,
    50.24E0,
    61.38E0,
    66.25E0,
    73.42E0,
    95.52E0,
   107.32E0,
   122.04E0,
   134.03E0,
   163.19E0,
   163.48E0,
   175.70E0,
   179.86E0,
   211.27E0,
   217.78E0,
   219.14E0,
   262.52E0,
   268.01E0,
   268.62E0,
   336.25E0,
   337.23E0,
   339.33E0,
   427.38E0,
   428.58E0,
   432.68E0,
   528.99E0,
   531.08E0,
   628.34E0,
   253.24E0,
   273.13E0,
   273.66E0,
   282.10E0,
   346.62E0,
   347.19E0,
   348.78E0,
   351.18E0,
   450.10E0,
   450.35E0,
   451.92E0,
   455.56E0,
   552.22E0,
   553.56E0,
   555.74E0,
   652.59E0,
   656.20E0,
    14.13E0,
    20.41E0,
    31.30E0,
    33.84E0,
    39.70E0,
    48.83E0,
    54.50E0,
    60.41E0,
    72.77E0,
    75.25E0,
    86.84E0,
    94.88E0,
    96.40E0,
   117.37E0,
   139.08E0,
   147.73E0,
   158.63E0,
   161.84E0,
   192.11E0,
   206.76E0,
   209.07E0,
   213.32E0,
   226.44E0,
   237.12E0,
   330.90E0,
   358.72E0,
   370.77E0,
   372.72E0,
   396.24E0,
   416.59E0,
   484.02E0,
   495.47E0,
   514.78E0,
   515.65E0,
   519.47E0,
   544.47E0,
   560.11E0,
   620.77E0,
    18.97E0,
    28.93E0,
    33.91E0,
    40.03E0,
    44.66E0,
    49.87E0,
    55.16E0,
    60.90E0,
    72.08E0,
    85.15E0,
    97.06E0,
   119.63E0,
   133.27E0,
   143.84E0,
   161.91E0,
   180.67E0,
   198.44E0,
   226.86E0,
   229.65E0,
   258.27E0,
   273.77E0,
   339.15E0,
   350.13E0,
   362.75E0,
   371.03E0,
   393.32E0,
   448.53E0,
   473.78E0,
   511.12E0,
   524.70E0,
   548.75E0,
   551.64E0,
   574.02E0,
   623.86E0,
    21.46E0,
    24.33E0,
    33.43E0,
    39.22E0,
    44.18E0,
    55.02E0,
    94.33E0,
    96.44E0,
   118.82E0,
   128.48E0,
   141.94E0,
   156.92E0,
   171.65E0,
   190.00E0,
   223.26E0,
   223.88E0,
   231.50E0,
   265.05E0,
   269.44E0,
   271.78E0,
   273.46E0,
   334.61E0,
   339.79E0,
   349.52E0,
   358.18E0,
   377.98E0,
   394.77E0,
   429.66E0,
   468.22E0,
   487.27E0,
   519.54E0,
   523.03E0,
   612.99E0,
   638.59E0,
   641.36E0,
   622.05E0,
   631.50E0,
   663.97E0,
   646.9E0 ,
   748.29E0,
   749.21E0,
   750.14E0,
   647.04E0,
   646.89E0,
   746.9E0 ,
   748.43E0,
   747.35E0,
   749.27E0,
   647.61E0,
   747.78E0,
   750.51E0,
   851.37E0,
   845.97E0,
   847.54E0,
   849.93E0,
   851.61E0,
   849.75E0,
   850.98E0,
   848.23E0
};


int
hahn1_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double b[7];
  size_t i;

  for (i = 0; i < 7; i++)
    {
      b[i] = gsl_vector_get(x, i);
    }

  for (i = 0; i < 236; i++)
    {
      double x = hahn1_F0[i];
      double y = ((b[0] + x* (b[1]  + x * (b[2] + x * b[3])))
                  / (1 + x*(b[4]  + x *(b[5] + x*b[6]))));
      gsl_vector_set (f, i, hahn1_F1[i] - y);
    }

  return GSL_SUCCESS;
}

int
hahn1_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double b[7];
  size_t i;

  for (i = 0; i < 7; i++)
    {
      b[i] = gsl_vector_get(x, i);
    }

  for (i = 0; i < 236; i++)
    {
      double x = hahn1_F0[i];
      double u = (b[0] + x*(b[1] + x*(b[2] + x * b[3])));
      double v = (1 + x*(b[4] + x*(b[5] + x*b[6])));
      gsl_matrix_set (df, i, 0, -1/v);
      gsl_matrix_set (df, i, 1, -x/v);
      gsl_matrix_set (df, i, 2, -x*x/v);
      gsl_matrix_set (df, i, 3, -x*x*x/v);
      gsl_matrix_set (df, i, 4, x*u/(v*v));
      gsl_matrix_set (df, i, 5, x*x*u/(v*v));
      gsl_matrix_set (df, i, 6, x*x*x*u/(v*v));
    }

  return GSL_SUCCESS;
}

int
hahn1_fdf (const gsl_vector * x, void *params,
	   gsl_vector * f, gsl_matrix * df)
{
  hahn1_f (x, params, f);
  hahn1_df (x, params, df);

  return GSL_SUCCESS;
}