/gsl/trunk

To get this branch, use:
bzr branch /lh/gsl/trunk
2193 by bjg
added gpl headers
1
/* interpolation/cspline.c
2
 * 
3648 by bjg
updated copyright notices
3
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2004 Gerard Jungman
2193 by bjg
added gpl headers
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
4212 by bjg
updated to gpl v3
7
 * the Free Software Foundation; either version 3 of the License, or (at
2193 by bjg
added gpl headers
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
3825 by bjg
updated FSF address in copyright notices
17
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
2193 by bjg
added gpl headers
18
 */
19
850 by jungman
work in progress
20
/* Author:  G. Jungman
21
 */
1484 by bjg
added #include <config.h> to all top-level source files
22
#include <config.h>
23
#include <stdlib.h>
1997 by jungman
added gsl/ and fixed build to use it
24
#include <gsl/gsl_errno.h>
25
#include <gsl/gsl_linalg.h>
26
#include <gsl/gsl_vector.h>
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
27
#include "integ_eval.h"
3009 by bjg
made #includes consistent to simplify build
28
#include <gsl/gsl_interp.h>
850 by jungman
work in progress
29
1460 by bjg
indentification
30
typedef struct
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
31
{
32
  double * c;
33
  double * g;
34
  double * diag;
35
  double * offdiag;
36
} cspline_state_t;
850 by jungman
work in progress
37
38
39
/* common initialization */
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
40
static void *
41
cspline_alloc (size_t size)
850 by jungman
work in progress
42
{
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
43
  cspline_state_t * state = (cspline_state_t *) malloc (sizeof (cspline_state_t));
44
45
  if (state == NULL)
46
    {
47
      GSL_ERROR_NULL("failed to allocate space for state", GSL_ENOMEM);
48
    }
49
  
50
  state->c = (double *) malloc (size * sizeof (double));
51
  
52
  if (state->c == NULL)
53
    {
54
      free (state);
55
      GSL_ERROR_NULL("failed to allocate space for c", GSL_ENOMEM);
56
    }
57
58
  state->g = (double *) malloc (size * sizeof (double));
59
  
60
  if (state->g == NULL)
61
    {
62
      free (state->c);
63
      free (state);
64
      GSL_ERROR_NULL("failed to allocate space for g", GSL_ENOMEM);
65
    }
66
67
  state->diag = (double *) malloc (size * sizeof (double));
68
  
69
  if (state->diag == NULL)
70
    {
71
      free (state->g);
72
      free (state->c);
73
      free (state);
74
      GSL_ERROR_NULL("failed to allocate space for diag", GSL_ENOMEM);
75
    }
76
77
  state->offdiag = (double *) malloc (size * sizeof (double));
78
  
79
  if (state->offdiag == NULL)
80
    {
81
      free (state->diag);
82
      free (state->g);
83
      free (state->c);
84
      free (state);
85
      GSL_ERROR_NULL("failed to allocate space for offdiag", GSL_ENOMEM);
86
    }
87
88
  return state;
850 by jungman
work in progress
89
}
90
91
880 by jungman
work in progress
92
/* natural spline calculation
93
 * see [Engeln-Mullges + Uhlig, p. 254]
94
 */
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
95
static int
96
cspline_init (void * vstate, const double xa[], const double ya[],
97
              size_t size)
880 by jungman
work in progress
98
{
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
99
  cspline_state_t *state = (cspline_state_t *) vstate;
100
1459 by bjg
clean up for make strict, use size_t instead of unsigned int
101
  size_t i;
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
102
  size_t num_points = size;
3262 by jungman
fixed incorrect min_size setting for cspline
103
  size_t max_index = num_points - 1;  /* Engeln-Mullges + Uhlig "n" */
104
  size_t sys_size = max_index - 1;    /* linear system is sys_size x sys_size */
105
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
106
  state->c[0] = 0.0;
107
  state->c[max_index] = 0.0;
3262 by jungman
fixed incorrect min_size setting for cspline
108
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
109
  for (i = 0; i < sys_size; i++)
110
    {
111
      const double h_i   = xa[i + 1] - xa[i];
112
      const double h_ip1 = xa[i + 2] - xa[i + 1];
113
      const double ydiff_i   = ya[i + 1] - ya[i];
114
      const double ydiff_ip1 = ya[i + 2] - ya[i + 1];
3720 by bjg
handle equal abscissae in cspline
115
      const double g_i = (h_i != 0.0) ? 1.0 / h_i : 0.0;
116
      const double g_ip1 = (h_ip1 != 0.0) ? 1.0 / h_ip1 : 0.0;
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
117
      state->offdiag[i] = h_ip1;
118
      state->diag[i] = 2.0 * (h_ip1 + h_i);
3720 by bjg
handle equal abscissae in cspline
119
      state->g[i] = 3.0 * (ydiff_ip1 * g_ip1 -  ydiff_i * g_i);
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
120
    }
3262 by jungman
fixed incorrect min_size setting for cspline
121
3586 by bjg
clean up tridiagonal solvers (fixes BUG#4)
122
  if (sys_size == 1)
123
    {
124
      state->c[1] = state->g[0] / state->diag[0];
125
      return GSL_SUCCESS;
126
    }
127
  else
128
    {
129
      gsl_vector_view g_vec = gsl_vector_view_array(state->g, sys_size);
130
      gsl_vector_view diag_vec = gsl_vector_view_array(state->diag, sys_size);
131
      gsl_vector_view offdiag_vec = gsl_vector_view_array(state->offdiag, sys_size - 1);
132
      gsl_vector_view solution_vec = gsl_vector_view_array ((state->c) + 1, sys_size);
133
      
134
      int status = gsl_linalg_solve_symm_tridiag(&diag_vec.vector, 
135
                                                 &offdiag_vec.vector, 
136
                                                 &g_vec.vector, 
137
                                                 &solution_vec.vector);
138
      return status;
139
    }
880 by jungman
work in progress
140
}
141
142
143
/* periodic spline calculation
144
 * see [Engeln-Mullges + Uhlig, p. 256]
145
 */
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
146
static int
147
cspline_init_periodic (void * vstate, const double xa[], const double ya[],
148
                       size_t size)
880 by jungman
work in progress
149
{
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
150
  cspline_state_t *state = (cspline_state_t *) vstate;
151
1459 by bjg
clean up for make strict, use size_t instead of unsigned int
152
  size_t i;
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
153
  size_t num_points = size;
1981 by jungman
uses gsl_linalg tridiag solver now
154
  size_t max_index = num_points - 1;  /* Engeln-Mullges + Uhlig "n" */
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
155
  size_t sys_size = max_index;    /* linear system is sys_size x sys_size */
1981 by jungman
uses gsl_linalg tridiag solver now
156
157
  if (sys_size == 2) {
158
    /* solve 2x2 system */
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
159
    
1981 by jungman
uses gsl_linalg tridiag solver now
160
    const double h0 = xa[1] - xa[0];
161
    const double h1 = xa[2] - xa[1];
3891 by bjg
fix invalid cspline_periodic memory access for n=3, added more tests
162
1981 by jungman
uses gsl_linalg tridiag solver now
163
    const double A = 2.0*(h0 + h1);
164
    const double B = h0 + h1;
165
    double g[2];
166
    double det;
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
167
    
1981 by jungman
uses gsl_linalg tridiag solver now
168
    g[0] = 3.0 * ((ya[2] - ya[1]) / h1 - (ya[1] - ya[0]) / h0);
3891 by bjg
fix invalid cspline_periodic memory access for n=3, added more tests
169
    g[1] = 3.0 * ((ya[1] - ya[2]) / h0 - (ya[2] - ya[1]) / h1);
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
170
    
1981 by jungman
uses gsl_linalg tridiag solver now
171
    det = 3.0 * (h0 + h1) * (h0 + h1);
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
172
    state->c[1] = ( A * g[0] - B * g[1])/det;
173
    state->c[2] = (-B * g[0] + A * g[1])/det;
174
    state->c[0] = state->c[2];
175
    
176
    return GSL_SUCCESS;
177
  } else {
178
    
3132 by bjg
Fixed bug in gsl_interp_cspline_periodic which caused a
179
    for (i = 0; i < sys_size-1; i++) {
180
      const double h_i       = xa[i + 1] - xa[i];
181
      const double h_ip1     = xa[i + 2] - xa[i + 1];
182
      const double ydiff_i   = ya[i + 1] - ya[i];
183
      const double ydiff_ip1 = ya[i + 2] - ya[i + 1];
3720 by bjg
handle equal abscissae in cspline
184
      const double g_i = (h_i != 0.0) ? 1.0 / h_i : 0.0;
185
      const double g_ip1 = (h_ip1 != 0.0) ? 1.0 / h_ip1 : 0.0;
3132 by bjg
Fixed bug in gsl_interp_cspline_periodic which caused a
186
      state->offdiag[i] = h_ip1;
187
      state->diag[i] = 2.0 * (h_ip1 + h_i);
3720 by bjg
handle equal abscissae in cspline
188
      state->g[i] = 3.0 * (ydiff_ip1 * g_ip1 - ydiff_i * g_i);
3132 by bjg
Fixed bug in gsl_interp_cspline_periodic which caused a
189
    }
190
191
    i = sys_size - 1;
192
193
    {
194
      const double h_i       = xa[i + 1] - xa[i];
195
      const double h_ip1     = xa[1] - xa[0];
196
      const double ydiff_i   = ya[i + 1] - ya[i];
197
      const double ydiff_ip1 = ya[1] - ya[0];
3720 by bjg
handle equal abscissae in cspline
198
      const double g_i = (h_i != 0.0) ? 1.0 / h_i : 0.0;
199
      const double g_ip1 = (h_ip1 != 0.0) ? 1.0 / h_ip1 : 0.0;
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
200
      state->offdiag[i] = h_ip1;
201
      state->diag[i] = 2.0 * (h_ip1 + h_i);
3720 by bjg
handle equal abscissae in cspline
202
      state->g[i] = 3.0 * (ydiff_ip1 * g_ip1 - ydiff_i * g_i);
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
203
    }
204
    
205
    {
2732 by bjg
fixed missing const
206
      gsl_vector_view g_vec = gsl_vector_view_array(state->g, sys_size);
207
      gsl_vector_view diag_vec = gsl_vector_view_array(state->diag, sys_size);
208
      gsl_vector_view offdiag_vec = gsl_vector_view_array(state->offdiag, sys_size);
209
      gsl_vector_view solution_vec = gsl_vector_view_array ((state->c) + 1, sys_size);
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
210
      
2732 by bjg
fixed missing const
211
      int status = gsl_linalg_solve_symm_cyc_tridiag(&diag_vec.vector, 
212
                                                     &offdiag_vec.vector, 
213
                                                     &g_vec.vector, 
214
                                                     &solution_vec.vector);
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
215
      state->c[0] = state->c[max_index];
216
      
217
      return status;
218
    }
1981 by jungman
uses gsl_linalg tridiag solver now
219
  }
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
220
}
221
222
223
static
224
void
225
cspline_free (void * vstate)
226
{
227
  cspline_state_t *state = (cspline_state_t *) vstate;
228
  
229
  free (state->c);
230
  free (state->g);
231
  free (state->diag);
232
  free (state->offdiag);
233
  free (state);
234
}
235
236
/* function for common coefficient determination
237
 */
238
static inline void
239
coeff_calc (const double c_array[], double dy, double dx, size_t index,  
240
            double * b, double * c, double * d)
241
{
242
  const double c_i = c_array[index];
243
  const double c_ip1 = c_array[index + 1];
244
  *b = (dy / dx) - dx * (c_ip1 + 2.0 * c_i) / 3.0;
245
  *c = c_i;
246
  *d = (c_ip1 - c_i) / (3.0 * dx);
247
}
248
249
250
static
251
int
252
cspline_eval (const void * vstate,
253
              const double x_array[], const double y_array[], size_t size,
254
              double x,
255
              gsl_interp_accel * a,
256
              double *y)
257
{
2732 by bjg
fixed missing const
258
  const cspline_state_t *state = (const cspline_state_t *) vstate;
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
259
260
  double x_lo, x_hi;
261
  double dx;
262
  size_t index;
263
  
264
  if (a != 0)
265
    {
266
      index = gsl_interp_accel_find (a, x_array, size, x);
267
    }
268
  else
269
    {
270
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
271
    }
272
  
273
  /* evaluate */
274
  x_hi = x_array[index + 1];
275
  x_lo = x_array[index];
276
  dx = x_hi - x_lo;
277
  if (dx > 0.0)
278
    {
279
      const double y_lo = y_array[index];
280
      const double y_hi = y_array[index + 1];
281
      const double dy = y_hi - y_lo;
282
      double delx = x - x_lo;
283
      double b_i, c_i, d_i; 
284
      coeff_calc(state->c, dy, dx, index,  &b_i, &c_i, &d_i);
285
      *y = y_lo + delx * (b_i + delx * (c_i + delx * d_i));
286
      return GSL_SUCCESS;
287
    }
288
  else
289
    {
290
      *y = 0.0;
291
      return GSL_EINVAL;
292
    }
293
}
294
295
296
static
297
int
298
cspline_eval_deriv (const void * vstate,
299
                    const double x_array[], const double y_array[], size_t size,
300
                    double x,
301
                    gsl_interp_accel * a,
302
                    double *dydx)
303
{
2732 by bjg
fixed missing const
304
  const cspline_state_t *state = (const cspline_state_t *) vstate;
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
305
306
  double x_lo, x_hi;
307
  double dx;
308
  size_t index;
309
  
310
  if (a != 0)
311
    {
312
      index = gsl_interp_accel_find (a, x_array, size, x);
313
    }
314
  else
315
    {
316
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
317
    }
318
  
319
  /* evaluate */
320
  x_hi = x_array[index + 1];
321
  x_lo = x_array[index];
322
  dx = x_hi - x_lo;
323
  if (dx > 0.0)
324
    {
325
      const double y_lo = y_array[index];
326
      const double y_hi = y_array[index + 1];
327
      const double dy = y_hi - y_lo;
328
      double delx = x - x_lo;
329
      double b_i, c_i, d_i; 
330
      coeff_calc(state->c, dy, dx, index,  &b_i, &c_i, &d_i);
331
      *dydx = b_i + delx * (2.0 * c_i + 3.0 * d_i * delx);
332
      return GSL_SUCCESS;
333
    }
334
  else
335
    {
336
      *dydx = 0.0;
4704 by Brian Gough
improve interpolation error handling (Evgeny Kurbatov <EvgenyKurbatov@yandex.ru>)
337
      return GSL_EINVAL;
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
338
    }
339
}
340
341
342
static
343
int
344
cspline_eval_deriv2 (const void * vstate,
345
                     const double x_array[], const double y_array[], size_t size,
346
                     double x,
347
                     gsl_interp_accel * a,
348
                     double * y_pp)
349
{
2732 by bjg
fixed missing const
350
  const cspline_state_t *state = (const cspline_state_t *) vstate;
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
351
352
  double x_lo, x_hi;
353
  double dx;
354
  size_t index;
355
  
356
  if (a != 0)
357
    {
358
      index = gsl_interp_accel_find (a, x_array, size, x);
359
    }
360
  else
361
    {
362
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
363
    }
364
  
365
  /* evaluate */
366
  x_hi = x_array[index + 1];
367
  x_lo = x_array[index];
368
  dx = x_hi - x_lo;
369
  if (dx > 0.0)
370
    {
371
      const double y_lo = y_array[index];
372
      const double y_hi = y_array[index + 1];
373
      const double dy = y_hi - y_lo;
374
      double delx = x - x_lo;
375
      double b_i, c_i, d_i;
376
      coeff_calc(state->c, dy, dx, index,  &b_i, &c_i, &d_i);
377
      *y_pp = 2.0 * c_i + 6.0 * d_i * delx;
378
      return GSL_SUCCESS;
379
    }
380
  else
381
    {
382
      *y_pp = 0.0;
4704 by Brian Gough
improve interpolation error handling (Evgeny Kurbatov <EvgenyKurbatov@yandex.ru>)
383
      return GSL_EINVAL;
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
384
    }
385
}
386
387
388
static
389
int
390
cspline_eval_integ (const void * vstate,
391
                    const double x_array[], const double y_array[], size_t size,
392
                    gsl_interp_accel * acc,
393
                    double a, double b,
394
                    double * result)
395
{
2732 by bjg
fixed missing const
396
  const cspline_state_t *state = (const cspline_state_t *) vstate;
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
397
398
  size_t i, index_a, index_b;
399
  
400
  if (acc != 0)
401
    {
402
      index_a = gsl_interp_accel_find (acc, x_array, size, a);
403
      index_b = gsl_interp_accel_find (acc, x_array, size, b);
404
    }
405
  else
406
    {
407
      index_a = gsl_interp_bsearch (x_array, a, 0, size - 1);
408
      index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
409
    }
410
411
  *result = 0.0;
412
  
413
  /* interior intervals */
414
  for(i=index_a; i<=index_b; i++) {
415
    const double x_hi = x_array[i + 1];
416
    const double x_lo = x_array[i];
417
    const double y_lo = y_array[i];
418
    const double y_hi = y_array[i + 1];
419
    const double dx = x_hi - x_lo;
420
    const double dy = y_hi - y_lo;
421
    if(dx != 0.0) {
422
      double b_i, c_i, d_i; 
423
      coeff_calc(state->c, dy, dx, i,  &b_i, &c_i, &d_i);
424
      
425
      if (i == index_a || i == index_b)
426
        {
427
          double x1 = (i == index_a) ? a : x_lo;
428
          double x2 = (i == index_b) ? b : x_hi;
429
          *result += integ_eval(y_lo, b_i, c_i, d_i, x_lo, x1, x2);
430
        }
431
      else
432
        {
433
          *result += dx * (y_lo + dx*(0.5*b_i + dx*(c_i/3.0 + 0.25*d_i*dx)));
434
        }
1460 by bjg
indentification
435
    }
1981 by jungman
uses gsl_linalg tridiag solver now
436
    else {
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
437
      *result = 0.0;
4704 by Brian Gough
improve interpolation error handling (Evgeny Kurbatov <EvgenyKurbatov@yandex.ru>)
438
      return GSL_EINVAL;
1460 by bjg
indentification
439
    }
1981 by jungman
uses gsl_linalg tridiag solver now
440
  }
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
441
  
442
  return GSL_SUCCESS;
443
}
444
445
static const gsl_interp_type cspline_type = 
446
{
447
  "cspline", 
3262 by jungman
fixed incorrect min_size setting for cspline
448
  3,
2709 by bjg
standardized to gsl conventions, added higher-level 'spline' interface
449
  &cspline_alloc,
450
  &cspline_init,
451
  &cspline_eval,
452
  &cspline_eval_deriv,
453
  &cspline_eval_deriv2,
454
  &cspline_eval_integ,
455
  &cspline_free
456
};
457
458
const gsl_interp_type * gsl_interp_cspline = &cspline_type;
459
460
static const gsl_interp_type cspline_periodic_type = 
461
{
462
  "cspline-periodic", 
463
  2,
464
  &cspline_alloc,
465
  &cspline_init_periodic,
466
  &cspline_eval,
467
  &cspline_eval_deriv,
468
  &cspline_eval_deriv2,
469
  &cspline_eval_integ,
470
  &cspline_free
471
};
472
473
const gsl_interp_type * gsl_interp_cspline_periodic = &cspline_periodic_type;
474
475