htswanalysis/MACS/lib/gsl/gsl-1.11/integration/qpsrt.c

   1 /* integration/qpsrt.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 static inline void
  21 qpsrt (gsl_integration_workspace * workspace);
  22
  23 static inline
  24 void qpsrt (gsl_integration_workspace * workspace)
  25 {
  26   const size_t last = workspace->size - 1;
  27   const size_t limit = workspace->limit;
  28
  29   double * elist = workspace->elist;
  30   size_t * order = workspace->order;
  31
  32   double errmax ;
  33   double errmin ;
  34   int i, k, top;
  35
  36   size_t i_nrmax = workspace->nrmax;
  37   size_t i_maxerr = order[i_nrmax] ;
  38
  39   /* Check whether the list contains more than two error estimates */
  40
  41   if (last < 2)
  42     {
  43       order[0] = 0 ;
  44       order[1] = 1 ;
  45       workspace->i = i_maxerr ;
  46       return ;
  47     }
  48
  49   errmax = elist[i_maxerr] ;
  50
  51   /* This part of the routine is only executed if, due to a difficult
  52      integrand, subdivision increased the error estimate. In the normal
  53      case the insert procedure should start after the nrmax-th largest
  54      error estimate. */
  55
  56   while (i_nrmax > 0 && errmax > elist[order[i_nrmax - 1]])
  57     {
  58       order[i_nrmax] = order[i_nrmax - 1] ;
  59       i_nrmax-- ;
  60     }
  61
  62   /* Compute the number of elements in the list to be maintained in
  63      descending order. This number depends on the number of
  64      subdivisions still allowed. */
  65
  66   if(last < (limit/2 + 2))
  67     {
  68       top = last ;
  69     }
  70   else
  71     {
  72       top = limit - last + 1;
  73     }
  74
  75   /* Insert errmax by traversing the list top-down, starting
  76      comparison from the element elist(order(i_nrmax+1)). */
  77
  78   i = i_nrmax + 1 ;
  79
  80   /* The order of the tests in the following line is important to
  81      prevent a segmentation fault */
  82
  83   while (i < top && errmax < elist[order[i]])
  84     {
  85       order[i-1] = order[i] ;
  86       i++ ;
  87     }
  88
  89   order[i-1] = i_maxerr ;
  90
  91   /* Insert errmin by traversing the list bottom-up */
  92
  93   errmin = elist[last] ;
  94
  95   k = top - 1 ;
  96
  97   while (k > i - 2 && errmin >= elist[order[k]])
  98     {
  99       order[k+1] = order[k] ;
 100       k-- ;
 101     }
 102
 103   order[k+1] = last ;
 104
 105   /* Set i_max and e_max */
 106
 107   i_maxerr = order[i_nrmax] ;
 108
 109   workspace->i = i_maxerr ;
 110   workspace->nrmax = i_nrmax ;
 111 }
 112
 113