1 #include <cmath>
2 #include "../gmx_lapack.h"
3 #include "lapack_limits.h"
4
5 #include "gromacs/utility/real.h"
6
7 void
F77_FUNC(dlartg,DLARTG)8 F77_FUNC(dlartg,DLARTG)(double *f,
9 double *g,
10 double *cs,
11 double *sn,
12 double *r)
13 {
14 double minval,safemin, safemin2, safemx2, eps;
15 double f1,g1,f1a,g1a,scale;
16 int i,n,count;
17
18 eps = GMX_DOUBLE_EPS;
19 minval = GMX_DOUBLE_MIN;
20 safemin = minval*(1.0+eps);
21 n = static_cast<int>(0.5*std::log( safemin/eps ) / std::log(2.0));
22 safemin2 = std::pow(2.0,static_cast<double>(n));
23
24 safemx2 = 1.0 / safemin2;
25
26 if(std::abs(*g)<GMX_DOUBLE_MIN) {
27 *cs = 1.0;
28 *sn = 0.0;
29 *r = *f;
30 } else if (std::abs(*f)<GMX_DOUBLE_MIN) {
31 *cs = 0.0;
32 *sn = 1.0;
33 *r = *g;
34 } else {
35 f1 = *f;
36 g1 = *g;
37 f1a = std::abs(f1);
38 g1a = std::abs(g1);
39 scale = (f1a > g1a) ? f1a : g1a;
40 if(scale >= safemx2) {
41 count = 0;
42 while(scale >= safemx2) {
43 count++;
44 f1 *= safemin2;
45 g1 *= safemin2;
46 f1a = std::abs(f1);
47 g1a = std::abs(g1);
48 scale = (f1a > g1a) ? f1a : g1a;
49 }
50 *r = std::sqrt(f1*f1 + g1*g1);
51 *cs = f1 / *r;
52 *sn = g1 / *r;
53 for(i=0;i<count;i++)
54 *r *= safemx2;
55 } else if (scale<=safemin2) {
56 count = 0;
57 while(scale <= safemin2) {
58 count++;
59 f1 *= safemx2;
60 g1 *= safemx2;
61 f1a = std::abs(f1);
62 g1a = std::abs(g1);
63 scale = (f1a > g1a) ? f1a : g1a;
64 }
65 *r = std::sqrt(f1*f1 + g1*g1);
66 *cs = f1 / *r;
67 *sn = g1 / *r;
68 for(i=0;i<count;i++)
69 *r *= safemin2;
70 } else {
71 *r = std::sqrt(f1*f1 + g1*g1);
72 *cs = f1 / *r;
73 *sn = g1 / *r;
74 }
75 if(std::abs(*f)>std::abs(*g) && *cs<0.0) {
76 *cs *= -1.0;
77 *sn *= -1.0;
78 *r *= -1.0;
79 }
80 }
81 return;
82 }
83
84