Skip to content

Commit 6940c59

Browse files
martin-frbgmartin-frbg
authored
Merge pull request #1454 from martin-frbg/issue1452
Keep the flag handling separate from the scaling loops in rotmg
2 parents 42514fa + 6500770 commit 6940c59

File tree

5 files changed

+253
-45
lines changed

5 files changed

+253
-45
lines changed

‎interface/rotmg.c‎

Lines changed: 24 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -136,7 +136,7 @@ void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
136136

137137
if(*dd1 != ZERO)
138138
{
139-
while( (*dd1 <= RGAMSQ) || (*dd1 >= GAMSQ) )
139+
if( (*dd1 <= RGAMSQ) || (*dd1 >= GAMSQ) )
140140
{
141141
if(dflag == ZERO)
142142
{
@@ -146,33 +146,34 @@ void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
146146
}
147147
else
148148
{
149-
if(dflag == ONE)
150-
{
151149
dh21 = -ONE;
152150
dh12 = ONE;
153151
dflag = -ONE;
154-
}
155152
}
156153
if( *dd1 <= RGAMSQ )
157154
{
158-
*dd1 = *dd1 * (GAM * GAM);
159-
*dx1 = *dx1 / GAM;
160-
dh11 = dh11 / GAM;
161-
dh12 = dh12 / GAM;
155+
while (ABS(*dd1) <= RGAMSQ) {
156+
*dd1 = *dd1 * (GAM * GAM);
157+
*dx1 = *dx1 / GAM;
158+
dh11 = dh11 / GAM;
159+
dh12 = dh12 / GAM;
160+
}
162161
}
163162
else
164163
{
165-
*dd1 = *dd1 / (GAM * GAM);
166-
*dx1 = *dx1 * GAM;
167-
dh11 = dh11 * GAM;
168-
dh12 = dh12 * GAM;
164+
while (ABS(*dd1) <= GAMSQ) {
165+
*dd1 = *dd1 / (GAM * GAM);
166+
*dx1 = *dx1 * GAM;
167+
dh11 = dh11 * GAM;
168+
dh12 = dh12 * GAM;
169+
}
169170
}
170171
}
171172
}
172173

173174
if(*dd2 != ZERO)
174175
{
175-
while( (ABS(*dd2) <= RGAMSQ) || (ABS(*dd2) >= GAMSQ) )
176+
if( (ABS(*dd2) <= RGAMSQ) || (ABS(*dd2) >= GAMSQ) )
176177
{
177178
if(dflag == ZERO)
178179
{
@@ -182,24 +183,25 @@ void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
182183
}
183184
else
184185
{
185-
if(dflag == ONE)
186-
{
187186
dh21 = -ONE;
188187
dh12 = ONE;
189188
dflag = -ONE;
190-
}
191189
}
192190
if( ABS(*dd2) <= RGAMSQ )
193191
{
194-
*dd2 = *dd2 * (GAM * GAM);
195-
dh21 = dh21 / GAM;
196-
dh22 = dh22 / GAM;
192+
while (ABS(*dd2) <= RGAMSQ) {
193+
*dd2 = *dd2 * (GAM * GAM);
194+
dh21 = dh21 / GAM;
195+
dh22 = dh22 / GAM;
196+
}
197197
}
198198
else
199199
{
200-
*dd2 = *dd2 / (GAM * GAM);
201-
dh21 = dh21 * GAM;
202-
dh22 = dh22 * GAM;
200+
while (ABS(*dd2) <= GAMSQ) {
201+
*dd2 = *dd2 / (GAM * GAM);
202+
dh21 = dh21 * GAM;
203+
dh22 = dh22 * GAM;
204+
}
203205
}
204206
}
205207
}

‎utest/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ else ()
77
set(OpenBLAS_utest_src
88
utest_main.c
99
test_amax.c
10+
test_rotmg.c
1011
)
1112
endif ()
1213

‎utest/Makefile‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,8 +8,8 @@ UTESTBIN=openblas_utest
88

99
include $(TOPDIR)/Makefile.system
1010

11-
OBJS=utest_main.o test_amax.o
12-
#test_rot.o test_swap.o test_axpy.o test_dotu.o test_rotmg.o test_dsdot.o test_fork.o
11+
OBJS=utest_main.o test_amax.o test_rotmg.o
12+
#test_rot.o test_swap.o test_axpy.o test_dotu.o test_dsdot.o test_fork.o
1313

1414
ifneq ($(NO_LAPACK), 1)
1515
#OBJS += test_potrs.o

‎utest/test_rotmg.c‎

Lines changed: 92 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -31,9 +31,9 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
3131
3232
**********************************************************************************/
3333

34-
#include "common_utest.h"
34+
#include "openblas_utest.h"
3535

36-
void test_drotmg()
36+
CTEST (drotmg,rotmg)
3737
{
3838
double te_d1, tr_d1;
3939
double te_d2, tr_d2;
@@ -42,31 +42,92 @@ void test_drotmg()
4242
double te_param[5];
4343
double tr_param[5];
4444
int i=0;
45-
te_d1= tr_d1=0.21149573940783739;
46-
te_d2= tr_d2=0.046892057172954082;
47-
te_x1= tr_x1=-0.42272687517106533;
48-
te_y1= tr_y1=0.42211309121921659;
45+
// original test case for libGoto bug fixed by feb2014 rewrite
46+
te_d1= 0.21149573940783739;
47+
te_d2= 0.046892057172954082;
48+
te_x1= -0.42272687517106533;
49+
te_y1= 0.42211309121921659;
50+
4951

5052
for(i=0; i<5; i++){
5153
te_param[i]=tr_param[i]=0.0;
5254
}
5355

56+
//reference values as calulated by netlib blas
57+
58+
tr_d1= 0.1732048;
59+
tr_d2= 0.03840234;
60+
tr_x1= -0.516180;
61+
tr_y1= 0.422113;
62+
tr_d1= 0.17320483687975;
63+
tr_d2= 0.03840233915037;
64+
tr_x1= -0.51618034832329;
65+
tr_y1= 0.42211309121922;
66+
67+
tr_param[0]= 0.0;
68+
tr_param[1]= 0.0;
69+
tr_param[2]= 0.99854803659786;
70+
tr_param[3]= -0.22139439665872;
71+
tr_param[4]= 0.0;
72+
73+
BLASFUNC(drotmg)(&te_d1, &te_d2, &te_x1, &te_y1, te_param);
74+
ASSERT_DBL_NEAR_TOL(te_d1, tr_d1, DOUBLE_EPS);
75+
ASSERT_DBL_NEAR_TOL(te_d2, tr_d2, DOUBLE_EPS);
76+
ASSERT_DBL_NEAR_TOL(te_x1, tr_x1, DOUBLE_EPS);
77+
ASSERT_DBL_NEAR_TOL(te_y1, tr_y1, DOUBLE_EPS);
78+
79+
for(i=0; i<5; i++){
80+
ASSERT_DBL_NEAR_TOL(te_param[i], tr_param[i], DOUBLE_EPS);
81+
}
82+
}
83+
84+
CTEST (drotmg,rotmg_issue1452)
85+
{
86+
double te_d1, tr_d1;
87+
double te_d2, tr_d2;
88+
double te_x1, tr_x1;
89+
double te_y1, tr_y1;
90+
double te_param[5];
91+
double tr_param[5];
92+
int i=0;
93+
94+
// from issue #1452, buggy version returned 0.000244 for param[3]
95+
te_d1 = 5.9e-8;
96+
te_d2 = 5.960464e-8;
97+
te_x1 = 1.0;
98+
te_y1 = 150.0;
99+
100+
for(i=0; i<5; i++){
101+
te_param[i]=tr_param[i]=0.0;
102+
}
103+
104+
//reference values as calulated by netlib blas
105+
tr_d1= 0.99995592822897;
106+
tr_d2= 0.98981219860583;
107+
tr_x1= 0.03662270484346;
108+
tr_y1= 150.000000000000;
109+
110+
tr_param[0]= -1.0;
111+
tr_param[1]= 0.00000161109346;
112+
tr_param[2]= -0.00024414062500;
113+
tr_param[3]= 1.0;
114+
tr_param[4]= 0.00000162760417;
115+
54116
//OpenBLAS
55117
BLASFUNC(drotmg)(&te_d1, &te_d2, &te_x1, &te_y1, te_param);
56-
//reference
57-
BLASFUNC_REF(drotmg)(&tr_d1, &tr_d2, &tr_x1, &tr_y1, tr_param);
58118

59-
CU_ASSERT_DOUBLE_EQUAL(te_d1, tr_d1, CHECK_EPS);
60-
CU_ASSERT_DOUBLE_EQUAL(te_d2, tr_d2, CHECK_EPS);
61-
CU_ASSERT_DOUBLE_EQUAL(te_x1, tr_x1, CHECK_EPS);
62-
CU_ASSERT_DOUBLE_EQUAL(te_y1, tr_y1, CHECK_EPS);
119+
ASSERT_DBL_NEAR_TOL(te_d1, tr_d1, DOUBLE_EPS);
120+
ASSERT_DBL_NEAR_TOL(te_d2, tr_d2, DOUBLE_EPS);
121+
ASSERT_DBL_NEAR_TOL(te_x1, tr_x1, DOUBLE_EPS);
122+
ASSERT_DBL_NEAR_TOL(te_y1, tr_y1, DOUBLE_EPS);
63123

64124
for(i=0; i<5; i++){
65-
CU_ASSERT_DOUBLE_EQUAL(te_param[i], tr_param[i], CHECK_EPS);
125+
ASSERT_DBL_NEAR_TOL(te_param[i], tr_param[i], DOUBLE_EPS);
66126
}
127+
67128
}
68129

69-
void test_drotmg_D1eqD2_X1eqX2()
130+
CTEST(drotmg, rotmg_D1eqD2_X1eqX2)
70131
{
71132
double te_d1, tr_d1;
72133
double te_d2, tr_d2;
@@ -83,18 +144,28 @@ void test_drotmg_D1eqD2_X1eqX2()
83144
for(i=0; i<5; i++){
84145
te_param[i]=tr_param[i]=0.0;
85146
}
147+
148+
//reference values as calulated by netlib blas
149+
tr_d1= 1.0;
150+
tr_d2= 1.0;
151+
tr_x1= 16.0;
152+
tr_y1= 8.0;
153+
154+
tr_param[0]=1.0;
155+
tr_param[1]=1.0;
156+
tr_param[2]=0.0;
157+
tr_param[3]=0.0;
158+
tr_param[4]=1.0;
86159

87160
//OpenBLAS
88161
BLASFUNC(drotmg)(&te_d1, &te_d2, &te_x1, &te_y1, te_param);
89-
//reference
90-
BLASFUNC_REF(drotmg)(&tr_d1, &tr_d2, &tr_x1, &tr_y1, tr_param);
91162

92-
CU_ASSERT_DOUBLE_EQUAL(te_d1, tr_d1, CHECK_EPS);
93-
CU_ASSERT_DOUBLE_EQUAL(te_d2, tr_d2, CHECK_EPS);
94-
CU_ASSERT_DOUBLE_EQUAL(te_x1, tr_x1, CHECK_EPS);
95-
CU_ASSERT_DOUBLE_EQUAL(te_y1, tr_y1, CHECK_EPS);
163+
ASSERT_DBL_NEAR_TOL(te_d1, tr_d1, DOUBLE_EPS);
164+
ASSERT_DBL_NEAR_TOL(te_d2, tr_d2, DOUBLE_EPS);
165+
ASSERT_DBL_NEAR_TOL(te_x1, tr_x1, DOUBLE_EPS);
166+
ASSERT_DBL_NEAR_TOL(te_y1, tr_y1, DOUBLE_EPS);
96167

97168
for(i=0; i<5; i++){
98-
CU_ASSERT_DOUBLE_EQUAL(te_param[i], tr_param[i], CHECK_EPS);
169+
ASSERT_DBL_NEAR_TOL(te_param[i], tr_param[i], DOUBLE_EPS);
99170
}
100171
}

0 commit comments

Comments
 (0)