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.qelg; 9 10 11 import std.algorithm: min, max; 12 import std.math; 13 14 import scid.core.fortran; 15 16 17 18 19 /// 20 void qelg(Real)(ref int n, Real* epstab_, ref Real result, 21 ref Real abserr, Real* res3la_, ref int nres) 22 { 23 //***begin prologue dqelg 24 //***refer to dqagie,dqagoe,dqagpe,dqagse 25 //***routines called d1mach 26 //***revision date 830518 (yymmdd) 27 //***keywords epsilon algorithm, convergence acceleration, 28 // extrapolation 29 //***author piessens,robert,appl. math. & progr. div. - k.u.leuven 30 // de doncker,elise,appl. math & progr. div. - k.u.leuven 31 //***purpose the routine determines the limit of a given sequence of 32 // approximations, by means of the epsilon algorithm of 33 // p.wynn. an estimate of the absolute error is also given. 34 // the condensed epsilon table is computed. only those 35 // elements needed for the computation of the next diagonal 36 // are preserved. 37 //***description 38 // 39 // epsilon algorithm 40 // standard fortran subroutine 41 // double precision version 42 // 43 // parameters 44 // n - integer 45 // epstab(n) contains the new element in the 46 // first column of the epsilon table. 47 // 48 // epstab - double precision 49 // vector of dimension 52 containing the elements 50 // of the two lower diagonals of the triangular 51 // epsilon table. the elements are numbered 52 // starting at the right-hand corner of the 53 // triangle. 54 // 55 // result - double precision 56 // resulting approximation to the integral 57 // 58 // abserr - double precision 59 // estimate of the absolute error computed from 60 // result and the 3 previous results 61 // 62 // res3la - double precision 63 // vector of dimension 3 containing the last 3 64 // results 65 // 66 // nres - integer 67 // number of calls to the routine 68 // (should be zero at first call) 69 // 70 //***end prologue dqelg 71 // 72 Real delta1,delta2,delta3, 73 epmach,epsinf,error,err1,err2,err3,e0,e1,e1abs,e2,e3, 74 oflow,res,ss,tol1,tol2,tol3; 75 int i=1,ib=1,ib2=1,ie=1,indx=1,k1=1,k2=1,k3=1,limexp=1,newelm=1,num=1; 76 auto epstab = dimension(epstab_, 52); 77 auto res3la = dimension(res3la_, 3); 78 // 79 // list of major variables 80 // ----------------------- 81 // 82 // e0 - the 4 elements on which the computation of a new 83 // e1 element in the epsilon table is based 84 // e2 85 // e3 e0 86 // e3 e1 new 87 // e2 88 // newelm - number of elements to be computed in the new 89 // diagonal 90 // error - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2) 91 // result - the element in the new diagonal with least value 92 // of error 93 // 94 // machine dependent constants 95 // --------------------------- 96 // 97 // epmach is the largest relative spacing. 98 // oflow is the largest positive magnitude. 99 // limexp is the maximum number of elements the epsilon 100 // table can contain. if this number is reached, the upper 101 // diagonal of the epsilon table is deleted. 102 // 103 //***first executable statement dqelg 104 epmach = Real.epsilon; 105 oflow = Real.max; 106 nres = nres+1; 107 abserr = oflow; 108 result = epstab[n]; 109 if (n < 3) goto l100; 110 limexp = 50; 111 epstab[n+2] = epstab[n]; 112 newelm = (n-1)/2; 113 epstab[n] = oflow; 114 num = n; 115 k1 = n; 116 for (i=1; i<=newelm; i++) { // end: l40 117 k2 = k1-1; 118 k3 = k1-2; 119 res = epstab[k1+2]; 120 e0 = epstab[k3]; 121 e1 = epstab[k2]; 122 e2 = res; 123 e1abs = fabs(e1); 124 delta2 = e2-e1; 125 err2 = fabs(delta2); 126 tol2 = max(fabs(e2),e1abs)*epmach; 127 delta3 = e1-e0; 128 err3 = fabs(delta3); 129 tol3 = max(e1abs,fabs(e0))*epmach; 130 if (err2 > tol2 || err3 > tol3) goto l10; 131 // 132 // if e0, e1 and e2 are equal to within machine 133 // accuracy, convergence is assumed. 134 // result = e2 135 // abserr = abs(e1-e0)+abs(e2-e1) 136 // 137 result = res; 138 abserr = err2+err3; 139 // ***jump out of do-loop 140 goto l100; 141 l10: e3 = epstab[k1]; 142 epstab[k1] = e1; 143 delta1 = e1-e3; 144 err1 = fabs(delta1); 145 tol1 = max(e1abs,fabs(e3))*epmach; 146 // 147 // if two elements are very close to each other, omit 148 // a part of the table by adjusting the value of n 149 // 150 if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3) goto l20; 151 ss = 0.1e1/delta1+0.1e1/delta2-0.1e1/delta3; 152 epsinf = fabs(ss*e1); 153 // 154 // test to detect irregular behaviour in the table, and 155 // eventually omit a part of the table adjusting the value 156 // of n. 157 // 158 if(epsinf > 0.1e-3) goto l30; 159 l20: n = i+i-1; 160 // ***jump out of do-loop 161 goto l50; 162 // 163 // compute a new element and eventually adjust 164 // the value of result. 165 // 166 l30: res = e1+0.1e1/ss; 167 epstab[k1] = res; 168 k1 = k1-2; 169 error = err2+fabs(res-e2)+err3; 170 if (error > abserr) goto l40; 171 abserr = error; 172 result = res; 173 l40:;} 174 // 175 // shift the table. 176 // 177 l50: if (n == limexp) n = 2*(limexp/2)-1; 178 ib = 1; 179 if ((num/2)*2 == num) ib = 2; 180 ie = newelm+1; 181 for (i=1; i<=ie; i++) { // end: l60 182 ib2 = ib+2; 183 epstab[ib] = epstab[ib2]; 184 ib = ib2; 185 l60:;} 186 if (num == n) goto l80; 187 indx = num-n+1; 188 for (i=1; i<=n; i++) { // end: l70 189 epstab[i]= epstab[indx]; 190 indx = indx+1; 191 l70:;} 192 l80: if (nres >= 4) goto l90; 193 res3la[nres] = result; 194 abserr = oflow; 195 goto l100; 196 // 197 // compute error estimate 198 // 199 l90: abserr = fabs(result-res3la[3])+fabs(result-res3la[2]) 200 +fabs(result-res3la[1]); 201 res3la[1] = res3la[2]; 202 res3la[2] = res3la[3]; 203 res3la[3] = result; 204 l100: abserr = max(abserr, 0.5e1*epmach*fabs(result)); 205 return; 206 } 207 208 209 unittest 210 { 211 alias qelg!float fqelg; 212 alias qelg!double dqelg; 213 alias qelg!real rqelg; 214 }