source: sasview/sansmodels/src/libigor/2Y_cpoly.c @ 9714ff5

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 9714ff5 was a42ea15, checked in by Jae Cho <jhjcho@…>, 12 years ago

Yukawa structure model kernels (IGOR/NIST)
  • Property mode set to 100644
File size: 17.1 KB
RevLine 
[a42ea15]1/*
2 *******************************************************************************
3 *
4 *
5 *                       Copyright (c) 2002
6 *                       Henrik Vestermark
7 *                       Denmark
8 *
9 *                       All Rights Reserved
10 *
11 *   This source file is subject to the terms and conditions of the
12 *   Henrik Vestermark Software License Agreement which restricts the manner
13 *   in which it may be used.
14 *
15 *******************************************************************************
16 *
17 *
18 * Module name     :   cpoly.cpp
19 * Module ID Nbr   :   
20 * Description     :   cpoly.cpp -- Jenkins-Traub real polynomial root finder.
21 *                     Translation of TOMS493 from FORTRAN to C. This
22 *                     implementation of Jenkins-Traub partially adapts
23 *                     the original code to a C environment by restruction
24 *                     many of the 'goto' controls to better fit a block
25 *                     structured form. It also eliminates the global memory
26 *                     allocation in favor of local, dynamic memory management.
27 *
28 *                     The calling conventions are slightly modified to return
29 *                     the number of roots found as the function value.
30 *
31 *                     INPUT:
32 *                     opr - double precision vector of real coefficients in order of
33 *                          decreasing powers.
34 *                     opi - double precision vector of imaginary coefficients in order of
35 *                          decreasing powers.
36 *                     degree - integer degree of polynomial
37 *
38 *                     OUTPUT:
39 *                     zeror,zeroi - output double precision vectors of the
40 *                          real and imaginary parts of the zeros.
41 *                            to be consistent with rpoly.cpp the zeros is inthe index
42 *                            [0..max_degree-1]
43 *
44 *                     RETURN:
45 *                     returnval:   -1 if leading coefficient is zero, otherwise
46 *                          number of roots found.
47 * --------------------------------------------------------------------------
48 * Change Record   :   
49 *
50 * Version      Author/Date             Description of changes
51 * -------  -----------         ----------------------
52 * 01.01                HVE/021101              Initial release
53 *
54 * -- SRK changed name of *pi to *piw (and all occurances in this file) to avoid the conflict with pi in Apple's fp.h
55 *
56 * End of Change Record
57 * --------------------------------------------------------------------------
58 */
59
60
61/* define version string */
62
63#include <stdlib.h>
64#include <math.h>
65#include <float.h>
66
67static double sr, si, tr, ti, pvr, pvi, are, mre, eta, infin;
68static int nn;
69static double *pr, *piw, *hr, *hi, *qpr, *qpi, *qhr, *qhi, *shr, *shi; 
70
71static void noshft( const int l1 );
72static void fxshft( const int l2, double *zr, double *zi, int *conv );
73static void vrshft( const int l3, double *zr, double *zi, int *conv );
74static void calct( int *bol );
75static void nexth( const int bol );
76static void polyev( const int nn, const double sr, const double si, const double pr[], const double piw[], double qr[], double qi[], double *pvr, double *pvi );
77static double errev( const int nn, const double qr[], const double qi[], const double ms, const double mp, const double are, const double mre );
78static void cauchy( const int nn, double pt[], double q[], double *fn_val );
79static double scale( const int nn, const double pt[], const double eta, const double infin, const double smalno, const double base );
80static void cdivid( const double ar, const double ai, const double br, const double bi, double *cr, double *ci );
81static double cmod( const double r, const double i );
82static void mcon( double *eta, double *infiny, double *smalno, double *base );
83
84int cpoly( const double *opr, const double *opi, int degree, double *zeror, double *zeroi ) 
85{
86        int cnt1, cnt2, idnn2, i, conv;
87        double xx, yy, cosr, sinr, smalno, base, xxx, zr, zi, bnd;
88       
89        mcon( &eta, &infin, &smalno, &base );
90        are = eta;
91        mre = 2.0 * sqrt( 2.0 ) * eta;
92        xx = 0.70710678;
93        yy = -xx;
94        cosr = -0.060756474;
95        sinr = -0.99756405;
96        nn = degree; 
97       
98        // Algorithm fails if the leading coefficient is zero
99        if( opr[ 0 ] == 0 && opi[ 0 ] == 0 )
100                return -1;
101       
102        // Allocate arrays
103        pr = malloc( sizeof( double ) * ( degree + 1 ) );
104        piw = malloc( sizeof( double ) * ( degree + 1 ) );
105        hr =malloc( sizeof( double ) * ( degree + 1 ) );
106        hi = malloc( sizeof( double ) * ( degree + 1 ) );
107        qpr= malloc( sizeof( double ) * ( degree + 1 ) );
108        qpi= malloc( sizeof( double ) * ( degree + 1 ) );
109        qhr= malloc( sizeof( double ) * ( degree + 1 ) );
110        qhi= malloc( sizeof( double ) * ( degree + 1 ) );
111        shr= malloc( sizeof( double ) * ( degree + 1 ) );
112        shi= malloc( sizeof( double ) * ( degree + 1 ) );
113       
114        // Remove the zeros at the origin if any
115        while( opr[ nn ] == 0 && opi[ nn ] == 0 )
116        {
117                idnn2 = degree - nn;
118                zeror[ idnn2 ] = 0;
119                zeroi[ idnn2 ] = 0;
120                nn--;
121        }
122       
123        // Make a copy of the coefficients
124        for( i = 0; i <= nn; i++ )
125        {
126                pr[ i ] = opr[ i ];
127                piw[ i ] = opi[ i ];
128                shr[ i ] = cmod( pr[ i ], piw[ i ] );
129        }
130       
131        // Scale the polynomial
132        bnd = scale( nn, shr, eta, infin, smalno, base );
133        if( bnd != 1 )
134                for( i = 0; i <= nn; i++ )
135                {
136                        pr[ i ] *= bnd;
137                        piw[ i ] *= bnd;
138                }
139       
140search: 
141        if( nn <= 1 )
142        {
143                cdivid( -pr[ 1 ], -piw[ 1 ], pr[ 0 ], piw[ 0 ], &zeror[ degree-1 ], &zeroi[ degree-1 ] );
144                goto finish;
145        }
146       
147        // Calculate bnd, alower bound on the modulus of the zeros
148        for( i = 0; i<= nn; i++ )
149                shr[ i ] = cmod( pr[ i ], piw[ i ] );
150       
151        cauchy( nn, shr, shi, &bnd );
152       
153        // Outer loop to control 2 Major passes with different sequences of shifts
154        for( cnt1 = 1; cnt1 <= 2; cnt1++ )
155        {
156                // First stage  calculation , no shift
157                noshft( 5 );
158               
159                // Inner loop to select a shift
160                for( cnt2 = 1; cnt2 <= 9; cnt2++ )
161                {
162                        // Shift is chosen with modulus bnd and amplitude rotated by 94 degree from the previous shif
163                        xxx = cosr * xx - sinr * yy;
164                        yy = sinr * xx + cosr * yy;
165                        xx = xxx;
166                        sr = bnd * xx;
167                        si = bnd * yy;
168                       
169                        // Second stage calculation, fixed shift
170                        fxshft( 10 * cnt2, &zr, &zi, &conv );
171                        if( conv )
172            {
173                                // The second stage jumps directly to the third stage ieration
174                                // If successful the zero is stored and the polynomial deflated
175                                idnn2 = degree - nn;
176                                zeror[ idnn2 ] = zr;
177                                zeroi[ idnn2 ] = zi;
178                                nn--;
179                                for( i = 0; i <= nn; i++ )
180                                {
181                                        pr[ i ] = qpr[ i ];
182                                        piw[ i ] = qpi[ i ];
183                                }
184                                goto search;
185            }
186                        // If the iteration is unsuccessful another shift is chosen
187                }
188                // if 9 shifts fail, the outer loop is repeated with another sequence of shifts
189        }
190       
191        // The zerofinder has failed on two major passes
192        // return empty handed with the number of roots found (less than the original degree)
193        degree -= nn;
194       
195finish:
196        // Deallocate arrays
197        free( pr );
198        free( piw );
199        free( hr );
200        free( hi );
201        free( qpr );
202        free( qpi );
203        free( qhr );
204        free( qhi );
205        free( shr );
206        free( shi );
207       
208        return degree;       
209}
210
211
212// COMPUTES  THE DERIVATIVE  POLYNOMIAL AS THE INITIAL H
213// POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS.
214//
215static void noshft( const int l1 )
216{
217        int i, j, jj, n, nm1;
218        double xni, t1, t2;
219       
220        n = nn;
221        nm1 = n - 1;
222        for( i = 0; i < n; i++ )
223        {
224                xni = nn - i;
225                hr[ i ] = xni * pr[ i ] / n;
226                hi[ i ] = xni * piw[ i ] / n;
227        }
228        for( jj = 1; jj <= l1; jj++ )
229        {
230                if( cmod( hr[ n - 1 ], hi[ n - 1 ] ) > eta * 10 * cmod( pr[ n - 1 ], piw[ n - 1 ] ) )
231                {
232                        cdivid( -pr[ nn ], -piw[ nn ], hr[ n - 1 ], hi[ n - 1 ], &tr, &ti );
233                        for( i = 0; i < nm1; i++ )
234            {
235                                j = nn - i - 1;
236                                t1 = hr[ j - 1 ];
237                                t2 = hi[ j - 1 ];
238                                hr[ j ] = tr * t1 - ti * t2 + pr[ j ];
239                                hi[ j ] = tr * t2 + ti * t1 + piw[ j ];
240            }
241                        hr[ 0 ] = pr[ 0 ];
242                        hi[ 0 ] = piw[ 0 ];
243                }
244                else
245                {
246                        // If the constant term is essentially zero, shift H coefficients
247                        for( i = 0; i < nm1; i++ )
248            {
249                                j = nn - i - 1;
250                                hr[ j ] = hr[ j - 1 ];
251                                hi[ j ] = hi[ j - 1 ];
252            }
253                        hr[ 0 ] = 0;
254                        hi[ 0 ] = 0;
255                }
256        }
257}
258
259// COMPUTES L2 FIXED-SHIFT H POLYNOMIALS AND TESTS FOR CONVERGENCE.
260// INITIATES A VARIABLE-SHIFT ITERATION AND RETURNS WITH THE
261// APPROXIMATE ZERO IF SUCCESSFUL.
262// L2 - LIMIT OF FIXED SHIFT STEPS
263// ZR,ZI - APPROXIMATE ZERO IF CONV IS .TRUE.
264// CONV  - LOGICAL INDICATING CONVERGENCE OF STAGE 3 ITERATION
265//
266static void fxshft( const int l2, double *zr, double *zi, int *conv )
267{
268        int i, j, n;
269        int test, pasd, bol;
270        double otr, oti, svsr, svsi;
271       
272        n = nn;
273        polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
274        test = 1;
275        pasd = 0;
276       
277        // Calculate first T = -P(S)/H(S)
278        calct( &bol );
279       
280        // Main loop for second stage
281        for( j = 1; j <= l2; j++ )
282        {
283                otr = tr;
284                oti = ti;
285               
286                // Compute the next H Polynomial and new t
287                nexth( bol );
288                calct( &bol );
289                *zr = sr + tr;
290                *zi = si + ti;
291               
292                // Test for convergence unless stage 3 has failed once or this
293                // is the last H Polynomial
294                if( !( bol || !test || j == 12 ) )
295                        if( cmod( tr - otr, ti - oti ) < 0.5 * cmod( *zr, *zi ) )
296            {
297                                if( pasd )
298                                {
299                                        // The weak convergence test has been passwed twice, start the third stage
300                                        // Iteration, after saving the current H polynomial and shift
301                                        for( i = 0; i < n; i++ )
302                                        {
303                                                shr[ i ] = hr[ i ];
304                                                shi[ i ] = hi[ i ];
305                                        }
306                                        svsr = sr;
307                                        svsi = si;
308                                        vrshft( 10, zr, zi, conv );
309                                        if( *conv ) return;
310                                       
311                                        //The iteration failed to converge. Turn off testing and restore h,s,pv and T
312                                        test = 0;
313                                        for( i = 0; i < n; i++ )
314                                        {
315                                                hr[ i ] = shr[ i ];
316                                                hi[ i ] = shi[ i ];
317                                        }
318                                        sr = svsr;
319                                        si = svsi;
320                                        polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
321                                        calct( &bol );
322                                        continue;
323                                }
324                                pasd = 1;
325            }
326                        else
327                                pasd = 0;
328        }
329       
330        // Attempt an iteration with final H polynomial from second stage
331        vrshft( 10, zr, zi, conv );
332}
333
334// CARRIES OUT THE THIRD STAGE ITERATION.
335// L3 - LIMIT OF STEPS IN STAGE 3.
336// ZR,ZI   - ON ENTRY CONTAINS THE INITIAL ITERATE, IF THE
337//           ITERATION CONVERGES IT CONTAINS THE FINAL ITERATE ON EXIT.
338// CONV    -  .TRUE. IF ITERATION CONVERGES
339//
340static void vrshft( const int l3, double *zr, double *zi, int *conv )
341{
342        int b, bol;
343        int i, j;
344        double mp, ms, omp, relstp, r1, r2, tp;
345       
346        *conv = 0;
347        b = 0;
348        sr = *zr;
349        si = *zi;
350       
351        // Main loop for stage three
352        for( i = 1; i <= l3; i++ )
353        {
354                // Evaluate P at S and test for convergence
355                polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
356                mp = cmod( pvr, pvi );
357                ms = cmod( sr, si );
358                if( mp <= 20 * errev( nn, qpr, qpi, ms, mp, are, mre ) )
359                {
360                        // Polynomial value is smaller in value than a bound onthe error
361                        // in evaluationg P, terminate the ietartion
362                        *conv = 1;
363                        *zr = sr;
364                        *zi = si;
365                        return;
366                }
367                if( i != 1 )
368                {
369                        if( !( b || mp < omp || relstp >= 0.05 ) )
370            {
371                                // Iteration has stalled. Probably a cluster of zeros. Do 5 fixed
372                                // shift steps into the cluster to force one zero to dominate
373                                tp = relstp;
374                                b = 1;
375                                if( relstp < eta ) tp = eta;
376                                r1 = sqrt( tp );
377                                r2 = sr * ( 1 + r1 ) - si * r1;
378                                si = sr * r1 + si * ( 1 + r1 );
379                                sr = r2;
380                                polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
381                                for( j = 1; j <= 5; j++ )
382                                {
383                                        calct( &bol );
384                                        nexth( bol );
385                                }
386                                omp = infin;
387                                goto _20;
388            }
389                       
390                        // Exit if polynomial value increase significantly
391                        if( mp *0.1 > omp ) return;
392                }
393               
394                omp = mp;
395               
396                // Calculate next iterate
397        _20:  calct( &bol );
398                nexth( bol );
399                calct( &bol );
400                if( !bol )
401                {
402                        relstp = cmod( tr, ti ) / cmod( sr, si );
403                        sr += tr;
404                        si += ti;
405                }
406        }
407}
408
409// COMPUTES  T = -P(S)/H(S).
410// BOOL   - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO.
411static void calct( int *bol )
412{
413        int n;
414        double hvr, hvi;
415       
416        n = nn;
417       
418        // evaluate h(s)
419        polyev( n - 1, sr, si, hr, hi, qhr, qhi, &hvr, &hvi );
420        *bol = cmod( hvr, hvi ) <= are * 10 * cmod( hr[ n - 1 ], hi[ n - 1 ] ) ? 1 : 0;
421        if( !*bol )
422        {
423                cdivid( -pvr, -pvi, hvr, hvi, &tr, &ti );
424                return;
425        }
426       
427        tr = 0;
428        ti = 0;
429}
430
431// CALCULATES THE NEXT SHIFTED H POLYNOMIAL.
432// BOOL   -  LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO
433//
434static void nexth( const int bol )
435{
436        int j, n;
437        double t1, t2;
438       
439        n = nn;
440        if( !bol )
441        {
442                for( j = 1; j < n; j++ )
443                {
444                        t1 = qhr[ j - 1 ];
445                        t2 = qhi[ j - 1 ];
446                        hr[ j ] = tr * t1 - ti * t2 + qpr[ j ];
447                        hi[ j ] = tr * t2 + ti * t1 + qpi[ j ];
448                }
449                hr[ 0 ] = qpr[ 0 ];
450                hi[ 0 ] = qpi[ 0 ];
451                return;
452        }
453       
454        // If h[s] is zero replace H with qh
455        for( j = 1; j < n; j++ )
456        {
457                hr[ j ] = qhr[ j - 1 ];
458                hi[ j ] = qhi[ j - 1 ];
459        }
460        hr[ 0 ] = 0;
461        hi[ 0 ] = 0;
462}
463
464// EVALUATES A POLYNOMIAL  P  AT  S  BY THE HORNER RECURRENCE
465// PLACING THE PARTIAL SUMS IN Q AND THE COMPUTED VALUE IN PV.
466// 
467static void polyev( const int nn, const double sr, const double si, const double pr[], const double piw[], double qr[], double qi[], double *pvr, double *pvi ) 
468{
469        int i;
470        double t;
471       
472        qr[ 0 ] = pr[ 0 ];
473        qi[ 0 ] = piw[ 0 ];
474        *pvr = qr[ 0 ];
475        *pvi = qi[ 0 ];
476       
477        for( i = 1; i <= nn; i++ )
478        {
479                t = ( *pvr ) * sr - ( *pvi ) * si + pr[ i ];
480                *pvi = ( *pvr ) * si + ( *pvi ) * sr + piw[ i ];
481                *pvr = t;
482                qr[ i ] = *pvr;
483                qi[ i ] = *pvi;
484        }
485}
486
487// BOUNDS THE ERROR IN EVALUATING THE POLYNOMIAL BY THE HORNER RECURRENCE.
488// QR,QI - THE PARTIAL SUMS
489// MS    -MODULUS OF THE POINT
490// MP    -MODULUS OF POLYNOMIAL VALUE
491// ARE, MRE -ERROR BOUNDS ON COMPLEX ADDITION AND MULTIPLICATION
492//
493static double errev( const int nn, const double qr[], const double qi[], const double ms, const double mp, const double are, const double mre )
494{
495        int i;
496        double e;
497       
498        e = cmod( qr[ 0 ], qi[ 0 ] ) * mre / ( are + mre );
499        for( i = 0; i <= nn; i++ )
500                e = e * ms + cmod( qr[ i ], qi[ i ] );
501       
502        return e * ( are + mre ) - mp * mre;
503}
504
505// CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A
506// POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS.
507//
508static void cauchy( const int nn, double pt[], double q[], double *fn_val )
509{
510        int i, n;
511        double x, xm, f, dx, df;
512       
513        pt[ nn ] = -pt[ nn ];
514       
515        // Compute upper estimate bound
516        n = nn;
517        x = exp( log( -pt[ nn ] ) - log( pt[ 0 ] ) ) / n;
518        if( pt[ n - 1 ] != 0 )
519        {
520                // Newton step at the origin is better, use it
521                xm = -pt[ nn ] / pt[ n - 1 ];
522                if( xm < x ) x = xm;
523        }
524       
525        // Chop the interval (0,x) until f < 0
526        while(1)
527        {
528                xm = x * 0.1;
529                f = pt[ 0 ];
530                for( i = 1; i <= nn; i++ )
531                        f = f * xm + pt[ i ];
532                if( f <= 0 )
533                        break;
534                x = xm;
535        }
536        dx = x;
537       
538        // Do Newton iteration until x converges to two decimal places
539        while( fabs( dx / x ) > 0.005 )
540        {
541                q[ 0 ] = pt[ 0 ];
542                for( i = 1; i <= nn; i++ )
543                        q[ i ] = q[ i - 1 ] * x + pt[ i ];
544                f = q[ nn ];
545                df = q[ 0 ];
546                for( i = 1; i < n; i++ )
547                        df = df * x + q[ i ];
548                dx = f / df;
549                x -= dx;
550        }
551       
552        *fn_val = x;
553}
554
555// RETURNS A SCALE FACTOR TO MULTIPLY THE COEFFICIENTS OF THE POLYNOMIAL.
556// THE SCALING IS DONE TO AVOID OVERFLOW AND TO AVOID UNDETECTED UNDERFLOW
557// INTERFERING WITH THE CONVERGENCE CRITERION.  THE FACTOR IS A POWER OF THE
558// BASE.
559// PT - MODULUS OF COEFFICIENTS OF P
560// ETA, INFIN, SMALNO, BASE - CONSTANTS DESCRIBING THE FLOATING POINT ARITHMETIC.
561//
562static double scale( const int nn, const double pt[], const double eta, const double infin, const double smalno, const double base )
563{
564        int i, l;
565        double hi, lo, max, min, x, sc;
566        double fn_val;
567       
568        // Find largest and smallest moduli of coefficients
569        hi = sqrt( infin );
570        lo = smalno / eta;
571        max = 0;
572        min = infin;
573       
574        for( i = 0; i <= nn; i++ )
575        {
576                x = pt[ i ];
577                if( x > max ) max = x;
578                if( x != 0 && x < min ) min = x;
579        }
580       
581        // Scale only if there are very large or very small components
582        fn_val = 1;
583        if( min >= lo && max <= hi ) return fn_val;
584        x = lo / min;
585        if( x <= 1 )
586                sc = 1 / ( sqrt( max )* sqrt( min ) );
587        else
588        {
589                sc = x;
590                if( infin / sc > max ) sc = 1;
591        }
592        l = (int)( log( sc ) / log(base ) + 0.5 );
593        fn_val = pow( base , l );
594        return fn_val;
595}
596
597// COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW.
598//
599static void cdivid( const double ar, const double ai, const double br, const double bi, double *cr, double *ci )
600{
601        double r, d, t, infin;
602       
603        if( br == 0 && bi == 0 )
604        {
605                // Division by zero, c = infinity
606                mcon( &t, &infin, &t, &t );
607                *cr = infin;
608                *ci = infin;
609                return;
610        }
611       
612        if( fabs( br ) < fabs( bi ) )
613        {
614                r = br/ bi;
615                d = bi + r * br;
616                *cr = ( ar * r + ai ) / d;
617                *ci = ( ai * r - ar ) / d;
618                return;
619        }
620       
621        r = bi / br;
622        d = br + r * bi;
623        *cr = ( ar + ai * r ) / d;
624        *ci = ( ai - ar * r ) / d;
625}
626
627// MODULUS OF A COMPLEX NUMBER AVOIDING OVERFLOW.
628//
629static double cmod( const double r, const double i )
630{
631        double ar, ai;
632       
633        ar = fabs( r );
634        ai = fabs( i );
635        if( ar < ai )
636                return ai * sqrt( 1.0 + pow( ( ar / ai ), 2.0 ) );
637       
638        if( ar > ai )
639                return ar * sqrt( 1.0 + pow( ( ai / ar ), 2.0 ) );
640       
641        return ar * sqrt( 2.0 );
642}
643
644// MCON PROVIDES MACHINE CONSTANTS USED IN VARIOUS PARTS OF THE PROGRAM.
645// THE USER MAY EITHER SET THEM DIRECTLY OR USE THE STATEMENTS BELOW TO
646// COMPUTE THEM. THE MEANING OF THE FOUR CONSTANTS ARE -
647// ETA       THE MAXIMUM RELATIVE REPRESENTATION ERROR WHICH CAN BE DESCRIBED
648//           AS THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
649//           1.0_dp + ETA &gt; 1.0.
650// INFINY    THE LARGEST FLOATING-POINT NUMBER
651// SMALNO    THE SMALLEST POSITIVE FLOATING-POINT NUMBER
652// BASE      THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED
653//
654static void mcon( double *eta, double *infiny, double *smalno, double *base )
655{
656        *base = 10;//DBL_RADIX;
657        *eta = DBL_EPSILON;
658        *infiny = DBL_MAX;
659        *smalno = DBL_MIN;
660}
Note: See TracBrowser for help on using the repository browser.