123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114 |
- /* integration/qpsrt.c
- *
- * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 3 of the License, or (at
- * your option) any later version.
- *
- * This program is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- */
- static inline void
- qpsrt (gsl_integration_workspace * workspace);
- static inline
- void qpsrt (gsl_integration_workspace * workspace)
- {
- const size_t last = workspace->size - 1;
- const size_t limit = workspace->limit;
- double * elist = workspace->elist;
- size_t * order = workspace->order;
- double errmax ;
- double errmin ;
- int i, k, top;
- size_t i_nrmax = workspace->nrmax;
- size_t i_maxerr = order[i_nrmax] ;
-
- /* Check whether the list contains more than two error estimates */
- if (last < 2)
- {
- order[0] = 0 ;
- order[1] = 1 ;
- workspace->i = i_maxerr ;
- return ;
- }
- errmax = elist[i_maxerr] ;
- /* This part of the routine is only executed if, due to a difficult
- integrand, subdivision increased the error estimate. In the normal
- case the insert procedure should start after the nrmax-th largest
- error estimate. */
- while (i_nrmax > 0 && errmax > elist[order[i_nrmax - 1]])
- {
- order[i_nrmax] = order[i_nrmax - 1] ;
- i_nrmax-- ;
- }
- /* Compute the number of elements in the list to be maintained in
- descending order. This number depends on the number of
- subdivisions still allowed. */
-
- if(last < (limit/2 + 2))
- {
- top = last ;
- }
- else
- {
- top = limit - last + 1;
- }
-
- /* Insert errmax by traversing the list top-down, starting
- comparison from the element elist(order(i_nrmax+1)). */
-
- i = i_nrmax + 1 ;
-
- /* The order of the tests in the following line is important to
- prevent a segmentation fault */
- while (i < top && errmax < elist[order[i]])
- {
- order[i-1] = order[i] ;
- i++ ;
- }
-
- order[i-1] = i_maxerr ;
-
- /* Insert errmin by traversing the list bottom-up */
-
- errmin = elist[last] ;
-
- k = top - 1 ;
-
- while (k > i - 2 && errmin >= elist[order[k]])
- {
- order[k+1] = order[k] ;
- k-- ;
- }
-
- order[k+1] = last ;
- /* Set i_max and e_max */
- i_maxerr = order[i_nrmax] ;
-
- workspace->i = i_maxerr ;
- workspace->nrmax = i_nrmax ;
- }
|