1 #include "common.h"
2 #ifdef FUNCTION_PROFILE
3 #include "functable.h"
4 #endif
5
6 #define GAM 4096.e0
7 #define GAMSQ 16777216.e0
8 #define RGAMSQ 5.9604645e-8
9
10 #ifndef CBLAS
11
NAME(FLOAT * dd1,FLOAT * dd2,FLOAT * dx1,FLOAT * DY1,FLOAT * dparam)12 void NAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT *DY1, FLOAT *dparam){
13
14 FLOAT dy1 = *DY1;
15
16 #else
17
18 void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
19
20 #endif
21
22 FLOAT du, dp1, dp2, dq2, dq1, dh11, dh21, dh12, dh22;
23 int igo, flag;
24 FLOAT dtemp;
25
26 #ifndef CBLAS
27 PRINT_DEBUG_NAME;
28 #else
29 PRINT_DEBUG_CNAME;
30 #endif
31
32 dh11 = ZERO;
33 dh12 = ZERO;
34 dh21 = ZERO;
35 dh22 = ZERO;
36
37 if (*dd1 < ZERO) goto L60;
38
39 dp2 = *dd2 * dy1;
40
41 if (dp2 == ZERO) {
42 flag = -2;
43 goto L260;
44 }
45
46 dp1 = *dd1 * *dx1;
47 dq2 = dp2 * dy1;
48 dq1 = dp1 * *dx1;
49
50 if (! (abs(dq1) > abs(dq2))) goto L40;
51
52 dh21 = -(dy1) / *dx1;
53 dh12 = dp2 / dp1;
54
55 du = ONE - dh12 * dh21;
56
57 if (du <= ZERO) goto L60;
58
59 flag = 0;
60 *dd1 /= du;
61 *dd2 /= du;
62 *dx1 *= du;
63
64 goto L100;
65
66 L40:
67 if (dq2 < ZERO) goto L60;
68
69 flag = 1;
70 dh11 = dp1 / dp2;
71 dh22 = *dx1 / dy1;
72 du = ONE + dh11 * dh22;
73 dtemp = *dd2 / du;
74 *dd2 = *dd1 / du;
75 *dd1 = dtemp;
76 *dx1 = dy1 * du;
77 goto L100;
78
79 L60:
80 flag = -1;
81 dh11 = ZERO;
82 dh12 = ZERO;
83 dh21 = ZERO;
84 dh22 = ZERO;
85
86 *dd1 = ZERO;
87 *dd2 = ZERO;
88 *dx1 = ZERO;
89 goto L220;
90
91
92 L70:
93 if (flag < 0) goto L90;
94
95 if (flag > 0) goto L80;
96
97 dh11 = ONE;
98 dh22 = ONE;
99 flag = -1;
100 goto L90;
101
102 L80:
103 dh21 = -ONE;
104 dh12 = ONE;
105 flag = -1;
106
107 L90:
108 switch (igo) {
109 case 0: goto L120;
110 case 1: goto L150;
111 case 2: goto L180;
112 case 3: goto L210;
113 }
114
115 L100:
116 if (!(*dd1 <= RGAMSQ)) goto L130;
117 if (*dd1 == ZERO) goto L160;
118 igo = 0;
119 goto L70;
120
121 L120:
122 *dd1 *= GAM * GAM;
123 *dx1 /= GAM;
124 dh11 /= GAM;
125 dh12 /= GAM;
126 goto L100;
127
128 L130:
129 if (! (*dd1 >= GAMSQ)) {
130 goto L160;
131 }
132 igo = 1;
133 goto L70;
134
135 L150:
136 *dd1 /= GAM * GAM;
137 *dx1 *= GAM;
138 dh11 *= GAM;
139 dh12 *= GAM;
140 goto L130;
141
142 L160:
143 if (! (abs(*dd2) <= RGAMSQ)) {
144 goto L190;
145 }
146 if (*dd2 == ZERO) {
147 goto L220;
148 }
149 igo = 2;
150 goto L70;
151
152 L180:
153 /* Computing 2nd power */
154 *dd2 *= GAM * GAM;
155 dh21 /= GAM;
156 dh22 /= GAM;
157 goto L160;
158
159 L190:
160 if (! (abs(*dd2) >= GAMSQ)) {
161 goto L220;
162 }
163 igo = 3;
164 goto L70;
165
166 L210:
167 /* Computing 2nd power */
168 *dd2 /= GAM * GAM;
169 dh21 *= GAM;
170 dh22 *= GAM;
171 goto L190;
172
173 L220:
174 if (flag < 0) {
175 goto L250;
176 } else if (flag == 0) {
177 goto L230;
178 } else {
179 goto L240;
180 }
181 L230:
182 dparam[2] = dh21;
183 dparam[3] = dh12;
184 goto L260;
185 L240:
186 dparam[2] = dh11;
187 dparam[4] = dh22;
188 goto L260;
189 L250:
190 dparam[1] = dh11;
191 dparam[2] = dh21;
192 dparam[3] = dh12;
193 dparam[4] = dh22;
194 L260:
195 dparam[0] = (FLOAT) flag;
196 return;
197 }
198
199
200