1 // This code has been mechanically translated from the original FORTRAN 2 // code at http://netlib.org/quadpack. 3 4 /** Authors: Lars Tandle Kyllingstad 5 Copyright: Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved. 6 License: Boost License 1.0 7 */ 8 module scid.ports.quadpack.qpsrt; 9 10 11 import scid.core.fortran; 12 13 14 15 16 /// 17 void qpsrt(Real)(int limit, int last, ref int maxerr, ref Real ermax, 18 Real* elist_, int* iord_, ref int nrmax) 19 { 20 //***begin prologue dqpsrt 21 //***refer to dqage,dqagie,dqagpe,dqawse 22 //***routines called (none) 23 //***revision date 810101 (yymmdd) 24 //***keywords sequential sorting 25 //***author piessens,robert,appl. math. & progr. div. - k.u.leuven 26 // de doncker,elise,appl. math. & progr. div. - k.u.leuven 27 //***purpose this routine maintains the descending ordering in the 28 // list of the local error estimated resulting from the 29 // interval subdivision process. at each call two error 30 // estimates are inserted using the sequential search 31 // method, top-down for the largest error estimate and 32 // bottom-up for the smallest error estimate. 33 //***description 34 // 35 // ordering routine 36 // standard fortran subroutine 37 // double precision version 38 // 39 // parameters (meaning at output) 40 // limit - integer 41 // maximum number of error estimates the list 42 // can contain 43 // 44 // last - integer 45 // number of error estimates currently in the list 46 // 47 // maxerr - integer 48 // maxerr points to the nrmax-th largest error 49 // estimate currently in the list 50 // 51 // ermax - double precision 52 // nrmax-th largest error estimate 53 // ermax = elist(maxerr) 54 // 55 // elist - double precision 56 // vector of dimension last containing 57 // the error estimates 58 // 59 // iord - integer 60 // vector of dimension last, the first k elements 61 // of which contain pointers to the error 62 // estimates, such that 63 // elist(iord(1)),..., elist(iord(k)) 64 // form a decreasing sequence, with 65 // k = last if last.le.(limit/2+2), and 66 // k = limit+1-last otherwise 67 // 68 // nrmax - integer 69 // maxerr = iord(nrmax) 70 // 71 //***end prologue dqpsrt 72 // 73 Real errmax,errmin; 74 int i=1,ibeg=1,ido=1,isucc=1,j=1,jbnd=1,jupbn=1,k=1; 75 auto elist = dimension(elist_, last); 76 auto iord = dimension(iord_, last); 77 // 78 // check whether the list contains more than 79 // two error estimates. 80 // 81 //***first executable statement dqpsrt 82 if (last > 2) goto l10; 83 iord[1] = 1; 84 iord[2] = 2; 85 goto l90; 86 // 87 // this part of the routine is only executed if, due to a 88 // difficult integrand, subdivision increased the error 89 // estimate. in the normal case the insert procedure should 90 // start after the nrmax-th largest error estimate. 91 // 92 l10: errmax = elist[maxerr]; 93 if (nrmax == 1) goto l30; 94 ido = nrmax-1; 95 for (i=1; i<=ido; i++) { // end: l20 96 isucc = iord[nrmax-1]; 97 // ***jump out of do-loop 98 if (errmax <= elist[isucc]) goto l30; 99 iord[nrmax] = isucc; 100 nrmax = nrmax-1; 101 l20:;} 102 // 103 // compute the number of elements in the list to be maintained 104 // in descending order. this number depends on the number of 105 // subdivisions still allowed. 106 // 107 l30: jupbn = last; 108 if (last > (limit/2+2)) jupbn = limit+3-last; 109 errmin = elist[last]; 110 // 111 // insert errmax by traversing the list top-down, 112 // starting comparison from the element elist(iord(nrmax+1)). 113 // 114 jbnd = jupbn-1; 115 ibeg = nrmax+1; 116 if (ibeg > jbnd) goto l50; 117 for (i=ibeg; i<=jbnd; i++) { // end: l40 118 isucc = iord[i]; 119 // ***jump out of do-loop 120 if (errmax >= elist[isucc]) goto l60; 121 iord[i-1] = isucc; 122 l40:;} 123 l50: iord[jbnd] = maxerr; 124 iord[jupbn] = last; 125 goto l90; 126 // 127 // insert errmin by traversing the list bottom-up. 128 // 129 l60: iord[i-1] = maxerr; 130 k = jbnd; 131 for (j=i; j<=jbnd; j++) { // end: l70 132 isucc = iord[k]; 133 // ***jump out of do-loop 134 if (errmin < elist[isucc]) goto l80; 135 iord[k+1] = isucc; 136 k = k-1; 137 l70:;} 138 iord[i] = last; 139 goto l90; 140 l80: iord[k+1] = last; 141 // 142 // set maxerr and ermax. 143 // 144 l90: maxerr = iord[nrmax]; 145 ermax = elist[maxerr]; 146 return; 147 } 148 149 unittest 150 { 151 alias qpsrt!float fqpsrt; 152 alias qpsrt!double dqpsrt; 153 alias qpsrt!real rqpsrt; 154 }