src/spline.c

   1 /* spline.c -- Do spline interpolation. */
   2
   3 /******************************************************************************
   4         David L. Fox
   5         2140 Braun Dr.
   6         Golden, CO 80401
   7
   8 This program has been placed in the public domain by its author.
   9
  10 Version Date            Change
  11 1.1     17 Dec 1985     Modify getdata() to realloc() more memory if needed.
  12 1.0     14 May 1985
  13
  14 spline [options] [file]
  15
  16 Spline reads pairs of numbers from file (or the standard input, if file is missing).
  17 The pairs are interpreted as abscissas and ordinates of a function.  The output of
  18 spline consists of similar pairs generated from the input by interpolation with
  19 cubic splines. (See R. W. Hamming, Numerical Methods for Scientists and Engineers,
  20 2nd ed., pp. 349ff.)  There are sufficiently many points in the output to appear
  21 smooth when plotted (e.g., by graph(1)).  The output points are approximately evenly
  22 spaced and include the input points.
  23
  24 There may be one or more options of the form: -o [argument [arg2] ].
  25
  26 The available options are:
  27
  28 -a      Generate abscissa values automatically.  The input consists of a list of
  29         ordinates.  The generated abscissas are evenly spaced with spacing given by
  30         the argument, or 1 if the next argument is not a number.
  31
  32 -f      Set the first derivative of the spline function at the left and right end
  33         points to the first and second arguments following -f.  If only one numerical
  34         argument follows -f then that value is used for both left and right end points.
  35
  36 -k      The argument of the k option is the value of the constant used to calculate
  37         the boundary values by y''[0] = ky''[1] and y''[n] = ky''[n-1].  The default
  38         value is k = 0.
  39
  40 -m      The argument gives the maximum number of input data points.  The default is 1000.
  41         If the input contains more than this many points additional memory will be
  42         allocated.  This option may be used to slightly increase efficiency for large
  43         data sets or reduce the amount of memory required for small data sets.
  44
  45 -n      The number of output points is given by the argument.  The default value is 100.
  46
  47 -p      The splines used for interpolation are forced to be periodic, i.e. y'[0] = y'[n].
  48         The first and last ordinates should be equal.
  49
  50 -s      Set the second derivative of the spline function at the left and right end
  51         points to the first and second arguments following -s.  If only one numerical
  52         argument follows -s then that value is used for both left and right end points.
  53
  54 -x      The argument (and arg2, if present) are the lower (and upper) limits on the
  55         abscissa values in the output.  If the x option is not specified these limits
  56         are calculated from the data.  Automatic abscissas start at the lower limit
  57         (default 0).
  58
  59 The data need not be monotonic in x but all the x values must be distinct.
  60 Non-numeric data in the input is ignored.
  61 ******************************************************************************/
  62
  63 #include        <a:stdio.h>
  64 #include        <ctype.h>
  65
  66 /* The constant DOUBLE is a compile time switch.
  67         If #define DOUBLE appears here double pecision variables are
  68         used to store all data and parameters.
  69         Otherwise, float variables are used for the data.
  70         For most smoothing and and interpolation applications single
  71         precision is more than adequate.
  72         Double precision is used to solve the system of linear equations
  73         in either case.
  74 */
  75 /* #define DOUBLE       */
  76
  77 #ifdef  DOUBLE
  78 #define double  real;           /* Type used for data storage. */
  79 #define IFMT    "%F"            /* Input format. */
  80 #define OFMT    "%18.14g %18.14g\n"     /* Output format. */
  81 #else
  82 #define float   real;           /* Type used for data storage. */
  83 #define IFMT    "%f"            /* Input format. */
  84 #define OFMT    "%8.5g %8.5g\n" /* Output format. */
  85 #endif
  86
  87 /* Numerical constants: These may be machine and/or precision dependent. */
  88 #define HUGE    (1.e38)         /* Largest floating point number. */
  89 #define EPS     1.e-5           /* Test for zero when solving linear equations. */
  90
  91 /* Default parameters */
  92 #define NPTS    1000
  93 #define NINT    100
  94 #define DEFSTEP 1.      /* Default step size for automatic abcissas. */
  95 #define DEFBVK  0.      /* Default boundary value constant. */
  96
  97 /* Boundary condition types. */
  98 #define EXTRAP  0       /* Extrapolate second derivative:
  99                            y''(0) = k * y''(1), y''(n) = k * y''(n-1) */
 100 #define FDERIV  1       /* Fixed first derivatives y'(0) and y'(n). */
 101 #define SDERIV  2       /* Fixed second derivatives y''(0) and y''(n). */
 102 #define PERIOD  3       /* Periodic:  derivatives equal at end points. */
 103
 104 /* Token types for command line processing. */
 105 #define OPTION  1
 106 #define NUMBER  2
 107 #define OTHER   3
 108
 109 /* Define error numbers. */
 110 #define MEMERR  1
 111 #define NODATA  2
 112 #define BADOPT  4
 113 #define BADFILE 5
 114 #define XTRAARG 6
 115 #define XTRAPTS 7
 116 #define SINGERR 8
 117 #define DUPX    9
 118 #define BADBC   10
 119 #define RANGERR 11
 120
 121 /* Constants and flags are global. */
 122 int aflag = FALSE;      /* Automatic abscissa flag. */
 123 real step = DEFSTEP;    /* Spacing for automatic abscissas. */
 124 int bdcode = EXTRAP;    /* Type of boundary conditions:
 125                         0 extrapolate f'' with coeficient k.
 126                         1 first derivatives given
 127                         2 second derivatives given
 128                         3 periodic */
 129 real leftd = 0.,        /* Boundary values of derivatives. */
 130    rightd = 0.;
 131 real k = DEFBVK;        /* Boundary value constant. */
 132 int rflag = 0;          /* 0: take range from data, 1: min. given, 2: min. & max. given. */
 133 real xmin = HUGE, xmax; /* Output range. */
 134 unsigned nint = NINT;   /* Number of intervals in output. */
 135 unsigned maxpts = NPTS; /* Maximun number of points. */
 136 unsigned nknots = 0;    /* Number of input points. */
 137 int sflag = FALSE;      /* If TRUE data must be sorted. */
 138 char *datafile;         /* Name of data file */
 139
 140 int xcompare();         /* Function to compare x values of two data points. */
 141
 142 struct pair {
 143         real x, y;
 144         } *data;        /* Points to list of data points. */
 145
 146 struct bandm {
 147         double diag, upper, rhs;
 148         } *eqns;        /* Points to elements of band matrix used to calculate
 149                         coefficients of the splines. */
 150
 151 main(argc, argv)
 152 int argc;
 153 char **argv;
 154 {
 155         setup(argc, argv);
 156         if (nknots == 1) {      /* Cannot interpolate with one data point.  */
 157                 printf(OFMT,data->x, data->y);
 158                 exit(0);
 159         }
 160         if (sflag)      /* Sort data if needed. */
 161                 qsort(data, nknots, sizeof(struct pair), xcompare);
 162         calcspline();
 163         interpolate();
 164 }
 165
 166 /* xcompare -- Compare abcissas of two data points (for qsort). */
 167 xcompare(arg1, arg2)
 168 struct pair *arg1, *arg2;
 169 {
 170         if (arg1->x > arg2->x)
 171                 return(1);
 172         else if (arg1->x < arg2->x)
 173                 return(-1);
 174         else
 175                 return(0);
 176 }
 177
 178 /* error -- Print error message and abort fatal errors. */
 179 error(errno, auxmsg)
 180 int errno;
 181 char *auxmsg;
 182 {       char *msg;
 183         int fatal, usemsg;
 184         static char *usage =
 185         "usage: spline [options] [file]\noptions:\n",
 186         *options = "-a spacing\n-k const\n-n intervals\n-m points\n-p\n-x xmin xmax\n";
 187
 188         fatal = TRUE;   /* Default is fatal error. */
 189         usemsg = FALSE; /* Default no usage message. */
 190         fprintf(stderr, "spline: ");
 191         switch(errno) {
 192         case MEMERR:
 193                 msg = "not enough memory for %u data points\n";
 194                 break;
 195         case NODATA:
 196                 msg = "data file %s empty\n";
 197                 break;
 198         case BADOPT:
 199                 msg = "unknown option: %c\n";
 200                 usemsg = TRUE;
 201                 break;
 202         case BADFILE:
 203                 msg = "cannot open file: %s\n";
 204                 break;
 205         case XTRAARG:
 206                 msg = "extra argument ignored: %s\n";
 207                 fatal = FALSE;
 208                 usemsg = TRUE;
 209                 break;
 210         case XTRAPTS:
 211                 fatal = FALSE;
 212                 msg = "%s";
 213                 break;
 214         case DUPX:
 215                 msg = "duplicate abcissa value: %s\n";
 216                 break;
 217         case RANGERR:
 218                 msg = "xmax < xmin not allowed %s\n";
 219                 break;
 220         /* The following errors "can't happen." */
 221         /* If they occur some sort of bug is indicated. */
 222         case SINGERR:
 223                 msg = "singular matrix encountered %s\n";
 224                 break;
 225         case BADBC:
 226                 msg = "internal error: bad boundary value code %d\n";
 227                 break;
 228         default:
 229                 fprintf(stderr, "unknown error number: %d\n", errno);
 230                 exit(1);
 231         }
 232         fprintf(stderr, msg, auxmsg);
 233         if (usemsg)
 234                 fprintf(stderr,"%s%s", usage, options);
 235         if (fatal)
 236                 exit(1);
 237 }
 238
 239 /* setup -- Initalize all constants and read data. */
 240 setup(argc, argv)
 241 int argc;
 242 char **argv;
 243 {       char *malloc();
 244         FILE fdinp,             /* Source of input. */
 245                 doarg();
 246
 247         fdinp = doarg(argc, argv);      /* Process command line arguments. */
 248
 249         /* Allocate memory for data and band matrix of coefficients. */
 250         if ((data = malloc(maxpts*sizeof(struct pair))) == NULL) {
 251                 error(MEMERR, (char *)maxpts);
 252         }
 253
 254         getdata(fdinp);         /* Read data from fdinp. */
 255         if (fdinp != stdin)
 256                 fclose(fdinp);  /* Close input data file. */
 257         if (nknots == 0) {
 258                 error(NODATA, datafile);
 259         }
 260         /* Allocate memory for calculation of spline coefficients. */
 261         if ((eqns = malloc((nknots+1)*sizeof(struct bandm))) == NULL) {
 262                 error(MEMERR, (char *)nknots);
 263         }
 264 }
 265
 266 /* doarg -- Process arguments. */
 267 FILE
 268 doarg(argc, argv)
 269 int argc;
 270 char **argv;
 271 {       int i, type;
 272         double atof();
 273         char *s, str[15];
 274         FILE fdinp;
 275
 276         s = argv[i=1];
 277         type = gettok(&i, &s, argv, str);
 278         do {
 279                 if (type == OPTION) {
 280                         switch(*str) {
 281                         case 'a':       /* Automatic abscissa. */
 282                                 aflag = TRUE;
 283                                 rflag = rflag < 2 ? 1 : rflag;
 284                                 if (xmin == HUGE)       /* Initialize xmin, if needed. */
 285                                         xmin = 0.;
 286                                 if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
 287                                         if ((step = atof(str)) <= 0.)
 288                                                 error(RANGERR,"");
 289                                         type = gettok(&i, &s, argv, str);
 290                                 }
 291                                 break;
 292                         case 'f':       /* Fix first derivative at boundaries. */
 293                                 bdcode = FDERIV;
 294                                 if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
 295                                         leftd = atof(str);
 296                                         if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
 297                                                 rightd = atof(str);
 298                                                 type = gettok(&i, &s, argv, str);
 299                                         }
 300                                         else {
 301                                                 rightd = leftd;
 302                                         }
 303                                 }
 304                                 break;
 305                         case 'k':       /* Set boundary value constant. */
 306                                 bdcode = EXTRAP;
 307                                 if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
 308                                         k = atof(str);
 309                                         type = gettok(&i, &s, argv, str);
 310                                 }
 311                                 break;
 312                         case 'm':       /* Set number of intervals in output. */
 313                                 if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
 314                                         maxpts = (unsigned)atoi(str);
 315                                         type = gettok(&i, &s, argv, str);
 316                                 }
 317                                 break;
 318                         case 'n':       /* Set number of intervals in output. */
 319                                 if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
 320                                         nint = (unsigned)atoi(str);
 321                                         type = gettok(&i, &s, argv, str);
 322                                 }
 323                                 break;
 324                         case 'p':       /* Require periodic interpolation function. */
 325                                 bdcode = PERIOD;
 326                                 type = gettok(&i, &s, argv, str);
 327                                 break;
 328                         case 's':       /* Fix second derivative at boundaries. */
 329                                 bdcode = SDERIV;
 330                                 if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
 331                                         leftd = atof(str);
 332                                         if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
 333                                                 rightd = atof(str);
 334                                                 type = gettok(&i, &s, argv, str);
 335                                         }
 336                                         else {
 337                                                 rightd = leftd;
 338                                         }
 339                                 }
 340                                 break;
 341                         case 'x':       /* Set range of x values. */
 342                                 rflag = 1;
 343                                 if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
 344                                         xmin = atof(str);
 345                                         if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
 346                                                 xmax = atof(str);
 347                                                 rflag = 2;
 348                                                 type = gettok(&i, &s, argv, str);
 349                                                 if (xmax <= xmin)
 350                                                         error(RANGERR, "");
 351                                         }
 352                                 }
 353                                 break;
 354                         default:
 355                                 error(BADOPT, (char *)argv[i][1]);
 356                         }
 357                 }
 358                 else {
 359                         if (argc > i) {
 360                                 datafile = argv[i];
 361                                 if ((fdinp = fopen(argv[i++], "r")) == NULL) {
 362                                         error(BADFILE, argv[i-1]);
 363                                 }
 364                                 if (argc > i)
 365                                         error(XTRAARG, argv[i]);
 366                         }
 367                         else
 368                                 fdinp = stdin;
 369                         break;
 370                 }
 371         } while (i < argc);
 372         return fdinp;
 373 }
 374
 375 /* gettok -- Get one token from command line, return type. */
 376 gettok(indexp, locp, argv, str)
 377 int *indexp;    /* Pointer to index in argv array. */
 378 char **locp;    /* Pointer to current location in argv[*indexp]. */
 379 char **argv;
 380 char *str;
 381 {       char *s;
 382         char *strcpy(), *strchr();
 383         int type;
 384
 385         s = *locp;
 386         while (isspace(*s) || *s == ',')
 387                 ++s;
 388         if (*s == '\0')         /* Look at next element in argv. */
 389                 s = argv[++*indexp];
 390         if (*s == '-' && isalpha(s[1])) {
 391                 /* Found an option. */
 392                 *str = *++s;
 393                 str[1] = '\0';
 394                 ++s;
 395                 type = OPTION;
 396         }
 397         else if (is_number(s)) {
 398                 while (is_number(s))
 399                         *str++ = *s++;
 400                 *str = '\0';
 401                 type = NUMBER;
 402         }
 403         else {
 404                 strcpy(str, s);
 405                 s = strchr(s, '\0');
 406                 type = OTHER;
 407         }
 408         *locp = s;
 409         return(type);
 410 }
 411
 412 /* is_number -- Return TRUE if argument is the ASCII representation of a number. */
 413 is_number(string)
 414 char *string;
 415 {
 416         if (isdigit(*string) ||
 417            *string == '.'    ||
 418            *string == '+'   ||
 419            (*string == '-' && (isdigit(string[1]) || string[1] == '.')))
 420                 return(TRUE);
 421         else
 422                 return(FALSE);
 423 }
 424
 425 /* getdata -- Read data points from fdinp. */
 426 getdata(fdinp)
 427 FILE fdinp;
 428 {       int n, i;
 429         real lastx, min, max;
 430         struct pair *dp;
 431         char msg[60], *realloc();
 432
 433         nknots = 0;
 434         lastx = -HUGE;
 435         dp = data;      /* Point to head of list of data. */
 436         min = HUGE;
 437         max = -HUGE;
 438         do {
 439                 if (aflag) {    /* Generate abcissa.  */
 440                         dp->x = xmin + nknots*step;
 441                         n = 0;
 442                 }
 443                 else {          /* Read abcissa. */
 444                         while ((n = (fscanf(fdinp,IFMT,&dp->x))) == 0)
 445                                 ;       /* Skip non-numeric input. */
 446                 }
 447                 if (n == 1) {
 448                         if (min > dp->x)
 449                                 min = dp->x;
 450                         if (max < dp->x)
 451                                 max = dp->x;
 452                         if (lastx > dp->x) {            /* Check for monotonicity. */
 453                                 sflag = TRUE;   /* Data must be sorted. */
 454                         }
 455                         lastx = dp->x;
 456                 }
 457                 /* Read ordinate. */
 458                 while ((n = (fscanf(fdinp, IFMT, &dp->y))) == 0)
 459                         ;       /* Skip non-numeric input. */
 460                 if (++nknots >= maxpts) {               /* Too many points, allocate more memory. */
 461                         if ((data = realloc(data, (maxpts *= 2)*sizeof(struct pair))) == NULL) {
 462                                 error(MEMERR, (char *)maxpts);
 463                         }
 464                         dp = data + nknots;
 465                         sprintf(msg, "more than %d input points, more memory allocated\n",
 466                            maxpts/2);
 467                         error(XTRAPTS,msg);
 468                 }
 469                 else
 470                         ++dp;
 471         } while (n == 1);
 472         --nknots;
 473         if (aflag) {    /* Compute maximum ordinate. */
 474                 max = xmin + (nknots-1)*step;
 475         }
 476         if (rflag < 2) {
 477                 xmax = max;
 478                 if (rflag < 1)
 479                         xmin = min;
 480         }
 481 }
 482
 483 /* calcspline -- Calculate coeficients of spline functions. */
 484 calcspline()
 485 {
 486         calccoef();
 487         solvband();
 488 }
 489
 490 /* calccoef -- Calculate coefficients of linear equations determining spline functions. */
 491 calccoef()
 492 {       int i, j;
 493         struct bandm *ptr, tmp1, tmp2;
 494         double h, h1;
 495         char str[10];
 496
 497         ptr = eqns;
 498         /* Initialize first row of matrix. */
 499         if ((h1 = data[1].x - data[0].x) == 0.) {
 500                 sprintf(str, "%8.5g", data[0].x);
 501                 error(DUPX, str);
 502         }
 503         switch(bdcode) {        /* First equation depends on boundary conditions. */
 504         case EXTRAP:
 505                 ptr->upper = -k;
 506                 ptr->diag = 1.;
 507                 ptr->rhs = 0.;
 508                 break;
 509         case FDERIV:            /* Fixed first derivatives at ends. */
 510                 ptr->upper = 1.;
 511                 ptr->diag = 2.;
 512                 h = data[1].x - data[0].x;
 513                 ptr->rhs = 6.*((data[1].y - data[0].y)/(h*h) - leftd/h);
 514                 break;
 515         case SDERIV:
 516                 ptr->upper = 0.;
 517                 ptr->diag = 1.;
 518                 ptr->rhs = leftd;
 519                 break;
 520         case PERIOD:            /* Periodic splines. */
 521                 ptr->upper = 1.;
 522                 ptr->diag = 2.;
 523                 break;
 524         default:
 525                 error(BADBC, (char *) bdcode);
 526         }
 527         ++ptr;
 528
 529         /* Initialize rows 1 to n-1, sub-diagonal band is assumed to be 1. */
 530         for (i=1; i     < nknots-1; ++i, ++ptr) {
 531                 h = h1;
 532                 if ((h1 = data[i+1].x - data[i].x) == 0.) {
 533                         sprintf(str, "%8.5g", data[i].x);
 534                         error(DUPX, str);
 535                 }
 536                 ptr->diag = 2.*(h + h1)/h;
 537                 ptr->upper = h1/h;
 538                 ptr->rhs = 6.*((data[i+1].y-data[i].y)/(h*h1) -
 539                    (data[i].y - data[i-1].y)/(h*h));
 540         }
 541         /* Initialize last row. */
 542         switch(bdcode) {
 543         case EXTRAP:            /* Standard case. */
 544                 ptr->diag = 1.;
 545                 ptr->upper = -k;        /* Upper holds actual sub-diagonal value. */
 546                 ptr->rhs = 0.;
 547                 break;
 548         case FDERIV:            /* Fixed first derivatives at ends. */
 549                 ptr->upper = 1.;
 550                 ptr->diag = 2.;
 551                 h = data[nknots-1].x - data[nknots-2].x;
 552                 ptr->rhs = 6.*((data[nknots-2].y - data[nknots-1].y)/(h*h) + rightd/h);
 553                 break;
 554         case SDERIV:
 555                 ptr->upper = 0.;
 556                 ptr->diag = 1.;
 557                 ptr->rhs = rightd;
 558                 break;
 559         case PERIOD:    /* Use periodic boundary conditions. */
 560                 /* First and last row are not in tridiagonal form. */
 561                 h = data[1].x - data[0].x;
 562                 h1 = data[nknots-1].x - data[nknots-2].x;
 563                 ptr->diag = 1.;
 564                 ptr->upper = 0.;
 565                 tmp1.diag = -1.;
 566                 tmp1.upper = 0;
 567                 tmp1.rhs = 0.;
 568                 tmp2.diag = 2.*h1/h;
 569                 tmp2.upper = h1/h;
 570                 tmp2.rhs = 6.*((data[1].y - data[0].y)/(h*h) -
 571                            (data[nknots-1].y - data[nknots-2].y)/(h1*h));
 572                 /* Transform periodic boundary equations to tri-diagonal form. */
 573                 for (i = 1; i < nknots - 1; ++i) {
 574                         tmp1.upper /= tmp1.diag;
 575                         tmp1.rhs /= tmp1.diag;
 576                         ptr->diag /= tmp1.diag;
 577                         ptr->upper /= tmp1.diag;
 578                         tmp1.diag = tmp1.upper - eqns[i].diag;
 579                         tmp1.upper = -eqns[i].upper;
 580                         tmp1.rhs -= eqns[i].rhs;
 581                         tmp2.upper /= tmp2.diag;
 582                         tmp2.rhs /= tmp2.diag;
 583                         eqns->diag /= tmp2.diag;
 584                         eqns->upper /= tmp2.diag;
 585                         tmp2.diag = tmp2.upper - eqns[nknots-1-i].diag/eqns[nknots-1-i].upper;
 586                         tmp2.upper = -1./eqns[nknots-1-i].upper;
 587                         tmp2.rhs -= eqns[nknots-1-i].rhs/eqns[nknots-1-i].upper;
 588                 }
 589                 /* Add in remaining terms of boundary condition equation. */
 590                 ptr->upper += tmp1.diag;
 591                 ptr->diag += tmp1.upper;
 592                 ptr->rhs = tmp1.rhs;
 593                 eqns->diag += tmp2.upper;
 594                 eqns->upper += tmp2.diag;
 595                 eqns->rhs = tmp2.rhs;
 596                 break;
 597         default:
 598                 error(BADBC, (char *) bdcode);
 599         }
 600 }
 601
 602 /* solvband -- Solve band matrix for spline functions. */
 603 solvband()
 604 {       int i, flag;
 605         struct bandm *ptr;
 606         double k1;
 607         double fabs();
 608         int fcompare();
 609
 610         ptr = eqns;
 611         flag = FALSE;
 612         /* Make a pass to triangularize matrix. */
 613         for (i=1; i < nknots - 1; ++i, ++ptr) {
 614                 if (fabs(ptr->diag) < EPS) {
 615                 /* Near zero on diagonal, pivot. */
 616                         if (fabs(ptr->upper) < EPS)
 617                                 error(SINGERR, "");
 618                         flag = TRUE;
 619                         ptr->diag = i;          /* Keep row index in diag.
 620                                                 Actual value of diag is always 1. */
 621                         if (i == nknots - 2) {
 622                                 flag = 2;
 623                                 /* Exchange next to last and last rows. */
 624                                 k1 = ptr->rhs/ptr->upper;
 625                                 if (fabs((ptr+1)->upper) < EPS)
 626                                         error(SINGERR, "");
 627                                 ptr->rhs = (ptr+1)->rhs/(ptr+1)->upper;
 628                                 ptr->upper = (ptr+1)->diag/(ptr+1)->upper;
 629                                 (ptr+1)->upper = 0.;
 630                                 (ptr+1)->rhs = k1;
 631                         }
 632                         else {
 633                                 ptr->rhs = (ptr+1)->rhs - (k1 = ptr->rhs/ptr->upper)*(ptr+1)->diag;
 634                                 ptr->upper = (ptr+1)->upper;    /* This isn't super-diagonal element
 635                                                                 but rather one to its right.
 636                                                                 Must watch for this when
 637                                                                 back substituting. */
 638                                 (++ptr)->diag = i-1;
 639                                 ++i;
 640                                 ptr->upper = 0;
 641                                 ptr->rhs = k1;
 642                                 (ptr+1)->rhs -= ptr->rhs;
 643                         }
 644                 }
 645                 else {
 646                         ptr->upper /= ptr->diag;
 647                         ptr->rhs /= ptr->diag;
 648                         ptr->diag = i-1;                /* Used to reorder solution if needed. */
 649                         (ptr+1)->diag -= ptr->upper;
 650                         (ptr+1)->rhs -= ptr->rhs;
 651                 }
 652         }
 653         /* Last row is a special case. */
 654         if (flag != 2) {
 655                 /* If flag == 2 last row is already computed. */
 656                 ptr->upper /= ptr->diag;
 657                 ptr->rhs /= ptr->diag;
 658                 ptr->diag = ptr - eqns;
 659                 ++ptr;
 660                 ptr->diag -= (ptr-1)->upper * ptr->upper;
 661                 ptr->rhs -= (ptr-1)->rhs * ptr->upper;
 662                 ptr->rhs /= ptr->diag;
 663                 ptr->diag = ptr - eqns;
 664         }
 665
 666         /* Now make a pass back substituting for solution. */
 667         --ptr;
 668         for ( ; ptr >= eqns; --ptr) {
 669                 if ((int)ptr->diag != ptr - eqns) {
 670                         /* This row and one above have been exchanged in a pivot. */
 671                         --ptr;
 672                         ptr->rhs -= ptr->upper * (ptr+2)->rhs;
 673                 }
 674                 else
 675                         ptr->rhs -= ptr->upper * (ptr+1)->rhs;
 676         }
 677         if (flag) {     /* Undo reordering done by pivoting. */
 678                 qsort(eqns, nknots, sizeof(struct bandm), fcompare);
 679         }
 680 }
 681
 682 /* fcompare -- Compare two floating point numbers. */
 683 fcompare(arg1, arg2)
 684 real *arg1, *arg2;
 685 {
 686         if (arg1 > arg2)
 687                 return(1);
 688         else if (arg1 < arg2)
 689                 return(-1);
 690         else
 691                 return(0);
 692 }
 693
 694 /* interpolate -- Do spline interpolation. */
 695 interpolate()
 696 {       int i;
 697         struct pair *dp;
 698         struct bandm *ep;
 699         real h, xi, yi, hi, xu, xl, limit;
 700
 701         h = (xmax - xmin)/nint;
 702         ep = eqns;
 703         dp = data + 1;
 704         for (xi = xmin; xi < xmax + 0.25*h; xi += h) {
 705                 while (dp->x < xi && dp < data + nknots - 1) {
 706                         ++dp;   /* Skip to correct interval. */
 707                         ++ep;
 708                 }
 709                 if (dp < data + nknots - 1 && dp->x < xmax)
 710                         limit = dp->x;
 711                 else
 712                         limit = xmax;
 713                 for ( ; xi < limit - 0.25*h; xi += h) {
 714                         /* Do interpolation. */
 715                         hi = dp->x - (dp-1)->x;
 716                         xu = dp->x - xi;
 717                         xl = xi - (dp-1)->x;
 718                         yi = ((ep+1)->rhs*xl*xl/(6.*hi) + dp->y/hi - (ep+1)->rhs*hi/6.)*xl +
 719                              (ep->rhs*xu*xu/(6.*hi) + (dp-1)->y/hi - ep->rhs*hi/6.)*xu;
 720                         printf(OFMT, xi, yi);
 721                 }
 722                 if (limit != dp->x) {   /* Interpolate. */
 723                         hi = dp->x - (dp-1)->x;
 724                         xu = dp->x - xmax;
 725                         xl = xmax - (dp-1)->x;
 726                         yi = ((ep+1)->rhs*xl*xl/(6.*hi) + dp->y/hi - (ep+1)->rhs*hi/6.)*xl +
 727                              (ep->rhs*xu*xu/(6.*hi) + (dp-1)->y/hi - ep->rhs*hi/6.)*xu;
 728                         printf(OFMT, xmax, yi);
 729                 }
 730                 else {          /* Print knot. */
 731                         printf(OFMT, dp->x, dp->y);
 732                         xi = dp->x;
 733                 }
 734         }
 735 }