Skip to content

Commit 63a8aef

Browse files
pvpv
committed
BUG: core/umath_linalg: ensure FP error status reflects LAPACK error status
1 parent bbdca51 commit 63a8aef

File tree

1 file changed

+69
-0
lines changed

1 file changed

+69
-0
lines changed

‎numpy/core/src/umath/umath_linalg.c.src‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -374,6 +374,12 @@ offset_ptr(void* ptr, ptrdiff_t offset)
374374
return (void*)((npy_uint8*)ptr + offset);
375375
}
376376

377+
static inline void
378+
clear_fp_errors()
379+
{
380+
PyUFunc_getfperr();
381+
}
382+
377383
/*
378384
*****************************************************************************
379385
** Some handy constants **
@@ -1780,6 +1786,7 @@ static inline void
17801786
size_t outer_dim = *dimensions++;
17811787
size_t op_count = (JOBZ=='N')?2:3;
17821788
EIGH_PARAMS_t eigh_params;
1789+
int error_occurred = 0;
17831790

17841791
for (iter=0; iter < op_count; ++iter) {
17851792
outer_steps[iter] = (ptrdiff_t) steps[iter];
@@ -1824,6 +1831,7 @@ static inline void
18241831
}
18251832
} else {
18261833
/* lapack fail, set result to nan */
1834+
error_occurred = 1;
18271835
nan_@BASETYPE@_matrix(args[1], &eigenvalues_out_ld);
18281836
if ('V' == eigh_params.JOBZ) {
18291837
nan_@TYPE@_matrix(args[2], &eigenvectors_out_ld);
@@ -1834,6 +1842,11 @@ static inline void
18341842

18351843
release_@lapack_func@(&eigh_params);
18361844
}
1845+
1846+
if (!error_occurred) {
1847+
/* lapack success, clear FP error status */
1848+
clear_fp_errors();
1849+
}
18371850
}
18381851
/**end repeat**/
18391852

@@ -1960,6 +1973,7 @@ static void
19601973
{
19611974
GESV_PARAMS_t params;
19621975
fortran_int n, nrhs;
1976+
int error_occurred = 0;
19631977
INIT_OUTER_LOOP_3
19641978

19651979
n = (fortran_int)dimensions[0];
@@ -1979,19 +1993,26 @@ static void
19791993
if (!not_ok) {
19801994
delinearize_@TYPE@_matrix(args[2], params.B, &r_out);
19811995
} else {
1996+
error_occurred = 1;
19821997
nan_@TYPE@_matrix(args[2], &r_out);
19831998
}
19841999
END_OUTER_LOOP
19852000

19862001
release_@lapack_func@(&params);
19872002
}
2003+
2004+
if (!error_occurred) {
2005+
/* lapack success, clear FP error status */
2006+
clear_fp_errors();
2007+
}
19882008
}
19892009

19902010
static void
19912011
@TYPE@_solve1(char **args, npy_intp *dimensions, npy_intp *steps,
19922012
void *NPY_UNUSED(func))
19932013
{
19942014
GESV_PARAMS_t params;
2015+
int error_occurred;
19952016
fortran_int n;
19962017
INIT_OUTER_LOOP_3
19972018

@@ -2010,12 +2031,18 @@ static void
20102031
if (!not_ok) {
20112032
delinearize_@TYPE@_matrix(args[2], params.B, &r_out);
20122033
} else {
2034+
error_occurred = 1;
20132035
nan_@TYPE@_matrix(args[2], &r_out);
20142036
}
20152037
END_OUTER_LOOP
20162038

20172039
release_@lapack_func@(&params);
20182040
}
2041+
2042+
if (!error_occurred) {
2043+
/* lapack success, clear FP error status */
2044+
clear_fp_errors();
2045+
}
20192046
}
20202047

20212048
static void
@@ -2024,6 +2051,7 @@ static void
20242051
{
20252052
GESV_PARAMS_t params;
20262053
fortran_int n;
2054+
int error_occurred = 0;
20272055
INIT_OUTER_LOOP_2
20282056

20292057
n = (fortran_int)dimensions[0];
@@ -2040,12 +2068,18 @@ static void
20402068
if (!not_ok) {
20412069
delinearize_@TYPE@_matrix(args[1], params.B, &r_out);
20422070
} else {
2071+
error_occurred = 1;
20432072
nan_@TYPE@_matrix(args[1], &r_out);
20442073
}
20452074
END_OUTER_LOOP
20462075

20472076
release_@lapack_func@(&params);
20482077
}
2078+
2079+
if (!error_occurred) {
2080+
/* lapack success, clear FP error status */
2081+
clear_fp_errors();
2082+
}
20492083
}
20502084
/**end repeat**/
20512085

@@ -2116,6 +2150,7 @@ static void
21162150
@TYPE@_cholesky(char uplo, char **args, npy_intp *dimensions, npy_intp *steps)
21172151
{
21182152
POTR_PARAMS_t params;
2153+
int error_occurred = 0;
21192154
fortran_int n;
21202155
INIT_OUTER_LOOP_2
21212156

@@ -2133,11 +2168,17 @@ static void
21332168
triu_@TYPE@_matrix(params.A, params.N);
21342169
delinearize_@TYPE@_matrix(args[1], params.A, &r_out);
21352170
} else {
2171+
error_occurred = 1;
21362172
nan_@TYPE@_matrix(args[1], &r_out);
21372173
}
21382174
END_OUTER_LOOP
21392175
release_@lapack_func@(&params);
21402176
}
2177+
2178+
if (!error_occurred) {
2179+
/* lapack success, clear FP error status */
2180+
clear_fp_errors();
2181+
}
21412182
}
21422183

21432184
static void
@@ -2551,6 +2592,7 @@ static inline void
25512592
size_t iter;
25522593
size_t outer_dim = *dimensions++;
25532594
size_t op_count = 2;
2595+
int error_occurred = 0;
25542596
GEEV_PARAMS_t geev_params;
25552597

25562598
STACK_TRACE;
@@ -2613,6 +2655,7 @@ static inline void
26132655
&vr_out);
26142656
} else {
26152657
/* geev failed */
2658+
error_occurred = 1;
26162659
nan_@COMPLEXTYPE@_matrix(*arg_iter++, &w_out);
26172660
if ('V' == geev_params.JOBVL)
26182661
nan_@COMPLEXTYPE@_matrix(*arg_iter++, &vl_out);
@@ -2624,6 +2667,11 @@ static inline void
26242667

26252668
release_@lapack_func@(&geev_params);
26262669
}
2670+
2671+
if (!error_occurred) {
2672+
/* lapack success, clear FP error status */
2673+
clear_fp_errors();
2674+
}
26272675
}
26282676

26292677
static void
@@ -2990,6 +3038,7 @@ static inline void
29903038
npy_intp* steps)
29913039
{
29923040
ptrdiff_t outer_steps[3];
3041+
int error_occurred = 0;
29933042
size_t iter;
29943043
size_t outer_dim = *dimensions++;
29953044
size_t op_count = (JOBZ=='N')?2:4;
@@ -3046,6 +3095,7 @@ static inline void
30463095
delinearize_@TYPE@_matrix(args[3], params.VT, &v_out);
30473096
}
30483097
} else {
3098+
error_occurred = 1;
30493099
if ('N' == params.JOBZ) {
30503100
nan_@REALTYPE@_matrix(args[1], &s_out);
30513101
} else {
@@ -3059,6 +3109,11 @@ static inline void
30593109

30603110
release_@lapack_func@(&params);
30613111
}
3112+
3113+
if (!error_occurred) {
3114+
/* lapack success, clear FP error status */
3115+
clear_fp_errors();
3116+
}
30623117
}
30633118
/**end repeat*/
30643119

@@ -3464,6 +3519,7 @@ static void
34643519
)
34653520
{
34663521
POTRS_PARAMS_t params;
3522+
int error_occurred = 0;
34673523
fortran_int n, nrhs;
34683524
INIT_OUTER_LOOP_3
34693525

@@ -3489,12 +3545,18 @@ static void
34893545
if (!not_ok) {
34903546
delinearize_@TYPE@_matrix(args[2], params.B, &r_out);
34913547
} else {
3548+
error_occurred = 1;
34923549
nan_@TYPE@_matrix(args[2], &r_out);
34933550
}
34943551
END_OUTER_LOOP
34953552

34963553
release_@lapack_func@(&params);
34973554
}
3555+
3556+
if (!error_occurred) {
3557+
/* lapack success, clear FP error status */
3558+
clear_fp_errors();
3559+
}
34983560
}
34993561

35003562
static void
@@ -3617,6 +3679,7 @@ static inline void
36173679
npy_intp *steps)
36183680
{
36193681
POTRI_PARAMS_t params;
3682+
int error_occurred = 0;
36203683
fortran_int n;
36213684
INIT_OUTER_LOOP_2
36223685

@@ -3635,12 +3698,18 @@ static inline void
36353698
make_symmetric_@TYPE@_matrix(params.UPLO, params.N, params.A);
36363699
delinearize_@TYPE@_matrix(args[1], params.A, &r_out);
36373700
} else {
3701+
error_occurred = 1;
36383702
nan_@TYPE@_matrix(args[1], &r_out);
36393703
}
36403704
END_OUTER_LOOP
36413705

36423706
release_@potri@(&params);
36433707
}
3708+
3709+
if (!error_occurred) {
3710+
/* lapack success, clear FP error status */
3711+
clear_fp_errors();
3712+
}
36443713
}
36453714

36463715
static void

0 commit comments

Comments
 (0)