Skip to content

Commit 80920a3

Browse files
martin-frbgmartin-frbg
authored
Merge pull request #5944 from lujiaweics/fix/zrotg-missing-profile-end
Fix build error in zrotg.c with FUNCTION_PROFILE=1 from Makefile.rule
2 parents 2621df7 + b6ad3e2 commit 80920a3

1 file changed

Lines changed: 137 additions & 132 deletions

File tree

‎interface/zrotg.c‎

Lines changed: 137 additions & 132 deletions
Original file line numberDiff line numberDiff line change
@@ -44,140 +44,145 @@ void CNAME(void *VDA, void *VDB, FLOAT *C, void *VS) {
4444

4545
FUNCTION_PROFILE_START();
4646

47-
if (db_r == ZERO && db_i == ZERO) {
48-
*C = ONE;
49-
*(S + 0) = ZERO;
50-
*(S + 1) = ZERO;
51-
return;
52-
}
47+
do {
48+
if (db_r == ZERO && db_i == ZERO) {
49+
*C = ONE;
50+
*(S + 0) = ZERO;
51+
*(S + 1) = ZERO;
52+
break;
53+
}
5354

54-
long double safmax = 1./safmin;
55-
#if defined DOUBLE
56-
long double rtmax = safmax /DBL_EPSILON;
57-
#else
58-
long double rtmax = safmax /FLT_EPSILON;
59-
#endif
60-
*(S1 + 0) = *(DB + 0);
61-
*(S1 + 1) = *(DB + 1) *-1;
62-
if (da_r == ZERO && da_i == ZERO) {
63-
*C = ZERO;
64-
if (db_r == ZERO) {
65-
(*DA) = fabsl(db_i);
66-
*S = *S1 /(*DA);
67-
*(S+1) = *(S1+1) /(*DA);
68-
return;
69-
} else if ( db_i == ZERO) {
70-
*DA = fabsl(db_r);
71-
*S = *S1 /(*DA);
72-
*(S+1) = *(S1+1) /(*DA);
73-
return;
74-
} else {
75-
long double g1 = MAX( fabsl(db_r), fabsl(db_i));
76-
rtmax =sqrt(safmax/2.);
77-
if (g1 > rtmin && g1 < rtmax) { // unscaled
78-
d = sqrt(adb);
79-
*S = *S1 /d;
80-
*(S+1) = *(S1+1) /d;
81-
*DA = d ;
82-
*(DA+1) = ZERO;
83-
return;
84-
} else { // scaled algorithm
85-
long double u = MIN ( safmax, MAX ( safmin, g1));
86-
FLOAT gs_r = db_r/u;
87-
FLOAT gs_i = db_i/u;
88-
d = sqrt ( gs_r*gs_r + gs_i*gs_i);
89-
*S = gs_r / d;
90-
*(S + 1) = (gs_i * -1) / d;
91-
*DA = d * u;
92-
*(DA+1) = ZERO;
93-
return;
94-
}
95-
}
96-
} else {
97-
FLOAT f1 = MAX ( fabsl(da_r), fabsl(da_i));
98-
FLOAT g1 = MAX ( fabsl(db_r), fabsl(db_i));
99-
rtmax = sqrt(safmax / 4.);
100-
if ( f1 > rtmin && f1 < rtmax && g1 > rtmin && g1 < rtmax) { //unscaled
101-
long double h = ada + adb;
102-
double adahsq = sqrt(ada * h);
103-
if (ada >= h *safmin) {
104-
*C = sqrt(ada/h);
105-
*R = *DA / *C;
106-
*(R+1) = *(DA+1) / *C;
107-
rtmax *= 2.;
108-
if ( ada > rtmin && h < rtmax) { // no risk of intermediate overflow
109-
*S = *S1 * (*DA / adahsq) - *(S1+1)* (*(DA+1)/adahsq);
110-
*(S+1) = *S1 * (*(DA+1) / adahsq) + *(S1+1) * (*DA/adahsq);
111-
} else {
112-
*S = *S1 * (*R/h) - *(S1+1) * (*(R+1)/h);
113-
*(S+1) = *S1 * (*(R+1)/h) + *(S1+1) * (*(R)/h);
55+
long double safmax = 1./safmin;
56+
#if defined DOUBLE
57+
long double rtmax = safmax /DBL_EPSILON;
58+
#else
59+
long double rtmax = safmax /FLT_EPSILON;
60+
#endif
61+
*(S1 + 0) = *(DB + 0);
62+
*(S1 + 1) = *(DB + 1) *-1;
63+
if (da_r == ZERO && da_i == ZERO) {
64+
*C = ZERO;
65+
if (db_r == ZERO) {
66+
(*DA) = fabsl(db_i);
67+
*S = *S1 /(*DA);
68+
*(S+1) = *(S1+1) /(*DA);
69+
break;
70+
} else if ( db_i == ZERO) {
71+
*DA = fabsl(db_r);
72+
*S = *S1 /(*DA);
73+
*(S+1) = *(S1+1) /(*DA);
74+
break;
75+
} else {
76+
long double g1 = MAX( fabsl(db_r), fabsl(db_i));
77+
rtmax =sqrt(safmax/2.);
78+
if (g1 > rtmin && g1 < rtmax) { // unscaled
79+
d = sqrt(adb);
80+
*S = *S1 /d;
81+
*(S+1) = *(S1+1) /d;
82+
*DA = d ;
83+
*(DA+1) = ZERO;
84+
break;
85+
} else { // scaled algorithm
86+
long double u = MIN ( safmax, MAX ( safmin, g1));
87+
FLOAT gs_r = db_r/u;
88+
FLOAT gs_i = db_i/u;
89+
d = sqrt ( gs_r*gs_r + gs_i*gs_i);
90+
*S = gs_r / d;
91+
*(S + 1) = (gs_i * -1) / d;
92+
*DA = d * u;
93+
*(DA+1) = ZERO;
94+
break;
95+
}
96+
}
97+
} else {
98+
FLOAT f1 = MAX ( fabsl(da_r), fabsl(da_i));
99+
FLOAT g1 = MAX ( fabsl(db_r), fabsl(db_i));
100+
rtmax = sqrt(safmax / 4.);
101+
if ( f1 > rtmin && f1 < rtmax && g1 > rtmin && g1 < rtmax) { //unscaled
102+
long double h = ada + adb;
103+
double adahsq = sqrt(ada * h);
104+
if (ada >= h *safmin) {
105+
*C = sqrt(ada/h);
106+
*R = *DA / *C;
107+
*(R+1) = *(DA+1) / *C;
108+
rtmax *= 2.;
109+
if ( ada > rtmin && h < rtmax) { // no risk of intermediate overflow
110+
*S = *S1 * (*DA / adahsq) - *(S1+1)* (*(DA+1)/adahsq);
111+
*(S+1) = *S1 * (*(DA+1) / adahsq) + *(S1+1) * (*DA/adahsq);
112+
} else {
113+
*S = *S1 * (*R/h) - *(S1+1) * (*(R+1)/h);
114+
*(S+1) = *S1 * (*(R+1)/h) + *(S1+1) * (*(R)/h);
115+
}
116+
} else {
117+
*C = ada / adahsq;
118+
if (*C >= safmin) {
119+
*R = *DA / *C;
120+
*(R+1) = *(DA+1) / *C;
121+
} else {
122+
*R = *DA * (h / adahsq);
123+
*(R+1) = *(DA+1) * (h / adahsq);
124+
}
125+
*S = *S1 * ada / adahsq;
126+
*(S+1) = *(S1+1) * ada / adahsq;
127+
}
128+
*DA=*R;
129+
*(DA+1)=*(R+1);
130+
break;
131+
} else { // scaled
132+
FLOAT fs_r, fs_i, gs_r, gs_i;
133+
long double v,w,f2,g2,h;
134+
long double u = MIN ( safmax, MAX ( safmin, MAX(f1,g1)));
135+
gs_r = db_r/u;
136+
gs_i = db_i/u;
137+
g2 = sqrt ( gs_r*gs_r + gs_i*gs_i);
138+
if (f1 /u < rtmin) {
139+
v = MIN (safmax, MAX (safmin, f1));
140+
w = v / u;
141+
fs_r = *DA/ v;
142+
fs_i = *(DA+1) / v;
143+
f2 = sqrt ( fs_r*fs_r + fs_i*fs_i);
144+
h = f2 * w * w + g2;
145+
} else { // use same scaling for both
146+
w = 1.;
147+
fs_r = *DA/ u;
148+
fs_i = *(DA+1) / u;
149+
f2 = sqrt ( fs_r*fs_r + fs_i*fs_i);
150+
h = f2 + g2;
151+
}
152+
if ( f2 >= h * safmin) {
153+
*C = sqrt ( f2 / h );
154+
*DA = fs_r / *C;
155+
*(DA+1) = fs_i / *C;
156+
rtmax *= 2;
157+
if ( f2 > rtmin && h < rtmax) {
158+
*S = gs_r * (fs_r /sqrt(f2*h)) - gs_i * (fs_i / sqrt(f2*h));
159+
*(S+1) = gs_r * (fs_i /sqrt(f2*h)) + gs_i * -1. * (fs_r / sqrt(f2*h));
160+
} else {
161+
*S = gs_r * (*DA/h) - gs_i * (*(DA+1) / h);
162+
*(S+1) = gs_r * (*(DA+1) /h) + gs_i * -1. * (*DA / h);
163+
}
164+
} else { // intermediates might overflow
165+
d = sqrt ( f2 * h);
166+
*C = f2 /d;
167+
if (*C >= safmin) {
168+
*DA = fs_r / *C;
169+
*(DA+1) = fs_i / *C;
170+
} else {
171+
*DA = fs_r * (h / d);
172+
*(DA+1) = fs_i / (h / d);
173+
}
174+
*S = gs_r * (fs_r /d) - gs_i * (fs_i / d);
175+
*(S+1) = gs_r * (fs_i /d) + gs_i * -1. * (fs_r / d);
176+
}
177+
*C *= w;
178+
*DA *= u;
179+
*(DA+1) *= u;
180+
break;
114181
}
115-
} else {
116-
*C = ada / adahsq;
117-
if (*C >= safmin) {
118-
*R = *DA / *C;
119-
*(R+1) = *(DA+1) / *C;
120-
} else {
121-
*R = *DA * (h / adahsq);
122-
*(R+1) = *(DA+1) * (h / adahsq);
123-
}
124-
*S = *S1 * ada / adahsq;
125-
*(S+1) = *(S1+1) * ada / adahsq;
126-
}
127-
*DA=*R;
128-
*(DA+1)=*(R+1);
129-
return;
130-
} else { // scaled
131-
FLOAT fs_r, fs_i, gs_r, gs_i;
132-
long double v,w,f2,g2,h;
133-
long double u = MIN ( safmax, MAX ( safmin, MAX(f1,g1)));
134-
gs_r = db_r/u;
135-
gs_i = db_i/u;
136-
g2 = sqrt ( gs_r*gs_r + gs_i*gs_i);
137-
if (f1 /u < rtmin) {
138-
v = MIN (safmax, MAX (safmin, f1));
139-
w = v / u;
140-
fs_r = *DA/ v;
141-
fs_i = *(DA+1) / v;
142-
f2 = sqrt ( fs_r*fs_r + fs_i*fs_i);
143-
h = f2 * w * w + g2;
144-
} else { // use same scaling for both
145-
w = 1.;
146-
fs_r = *DA/ u;
147-
fs_i = *(DA+1) / u;
148-
f2 = sqrt ( fs_r*fs_r + fs_i*fs_i);
149-
h = f2 + g2;
150-
}
151-
if ( f2 >= h * safmin) {
152-
*C = sqrt ( f2 / h );
153-
*DA = fs_r / *C;
154-
*(DA+1) = fs_i / *C;
155-
rtmax *= 2;
156-
if ( f2 > rtmin && h < rtmax) {
157-
*S = gs_r * (fs_r /sqrt(f2*h)) - gs_i * (fs_i / sqrt(f2*h));
158-
*(S+1) = gs_r * (fs_i /sqrt(f2*h)) + gs_i * -1. * (fs_r / sqrt(f2*h));
159-
} else {
160-
*S = gs_r * (*DA/h) - gs_i * (*(DA+1) / h);
161-
*(S+1) = gs_r * (*(DA+1) /h) + gs_i * -1. * (*DA / h);
162-
}
163-
} else { // intermediates might overflow
164-
d = sqrt ( f2 * h);
165-
*C = f2 /d;
166-
if (*C >= safmin) {
167-
*DA = fs_r / *C;
168-
*(DA+1) = fs_i / *C;
169-
} else {
170-
*DA = fs_r * (h / d);
171-
*(DA+1) = fs_i / (h / d);
172-
}
173-
*S = gs_r * (fs_r /d) - gs_i * (fs_i / d);
174-
*(S+1) = gs_r * (fs_i /d) + gs_i * -1. * (fs_r / d);
175-
}
176-
*C *= w;
177-
*DA *= u;
178-
*(DA+1) *= u;
179-
return;
180182
}
181-
}
183+
} while (0);
184+
185+
FUNCTION_PROFILE_END(4, 4, 4);
186+
IDEBUG_END;
182187
}
183188

0 commit comments

Comments
 (0)