| /* integration/qpsrt.c |
| * |
| * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 2 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 ; |
| } |
| |
| |