From 908c30ceb65c9b79add6c07c7e84adc5722f2334 Mon Sep 17 00:00:00 2001 From: Madan mohan Manokar Date: Fri, 21 Aug 2026 17:56:24 +0530 Subject: [PATCH] Merge pull request #29727 from amd:imp_jacobisvd_2 core: fix JacobiSVD SIMD accumulation to match scalar path - #29727 Replace FMA with mul-add in dotD/givensD to avoid Windows MSVC rounding drift. - Address the FMA drift introduced in https://github.com/opencv/opencv/pull/29720 ### Pull Request Readiness Checklist See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request - [x] I agree to contribute to the project under Apache 2 License. - [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV - [x] The PR is proposed to the proper branch - [ ] There is a reference to the original bug report and related work - [x] There is accuracy test, performance test and test data in opencv_extra repository, if applicable Patch to opencv_extra has the same branch name. - [x] The feature is well documented and sample code can be built with the project CMake --- modules/core/src/lapack.simd.hpp | 289 +++++++++++++++++++------------ modules/photo/test/test_hdr.cpp | 4 + 2 files changed, 182 insertions(+), 111 deletions(-) diff --git a/modules/core/src/lapack.simd.hpp b/modules/core/src/lapack.simd.hpp index 8c50480add..1dc5428e27 100644 --- a/modules/core/src/lapack.simd.hpp +++ b/modules/core/src/lapack.simd.hpp @@ -76,9 +76,9 @@ template static inline _Tp hypot(_Tp a, _Tp b) template struct VBLAS { int dot(const T*, const T*, int, T*) const { return 0; } - int givens(T*, T*, int, T, T) const { return 0; } - int dotD(const T*, const T*, int, double*) const { return 0; } - int givensD(T*, T*, int, T, T, double*, double*) const { return 0; } + void givens(T*, T*, int, T, T) const {} + void dotD(const T*, const T*, int, double*) const {} + void givensD(T*, T*, int, T, T, double*, double*) const {} }; #if CV_SIMD // TODO: enable for CV_SIMD_SCALABLE, GCC 13 related @@ -101,79 +101,149 @@ template<> inline int VBLAS::dot(const float* a, const float* b, int n, f } -template<> inline int VBLAS::givens(float* a, float* b, int n, float c, float s) const +template<> inline void VBLAS::givens(float* a, float* b, int n, float c, float s) const { - if( n < VTraits::vlanes()) - return 0; - int k = 0; + if( n <= 0 ) + return; + const int vl = VTraits::vlanes(); v_float32 c4 = vx_setall_f32(c), s4 = vx_setall_f32(s); v_float32 ns4 = vx_setall_f32(-s); - for( ; k <= n - VTraits::vlanes(); k += VTraits::vlanes() ) + int k = 0; + for( ; k <= n - vl; k += vl ) { v_float32 a0 = vx_load(a + k); v_float32 b0 = vx_load(b + k); - v_float32 t0 = v_fma(a0, c4, v_mul(b0, s4)); - v_float32 t1 = v_fma(a0, ns4, v_mul(b0, c4)); + v_float32 t0 = v_add(v_mul(a0, c4), v_mul(b0, s4)); + v_float32 t1 = v_add(v_mul(a0, ns4), v_mul(b0, c4)); v_store(a + k, t0); v_store(b + k, t1); } + for( ; k < n; k++ ) + { + float t0 = c*a[k] + s*b[k]; + float t1 = -s*a[k] + c*b[k]; + a[k] = t0; b[k] = t1; + } vx_cleanup(); - return k; } #if (CV_SIMD_64F || CV_SIMD_SCALABLE_64F) -template<> inline int VBLAS::dotD(const float* a, const float* b, int n, double* result) const +template<> inline void VBLAS::dotD(const float* a, const float* b, int n, double* result) const { - if( n < 2*VTraits::vlanes() ) - return 0; - int k = 0; + if( n <= 0 ) + return; + const int vl = VTraits::vlanes(); v_float64 s0 = vx_setzero_f64(), s1 = vx_setzero_f64(); - for( ; k <= n - VTraits::vlanes(); k += VTraits::vlanes() ) + int k = 0; + for( ; k <= n - vl; k += vl ) { v_float32 a0 = vx_load(a + k); v_float32 b0 = vx_load(b + k); - s0 = v_fma(v_cvt_f64(a0), v_cvt_f64(b0), s0); - s1 = v_fma(v_cvt_f64_high(a0), v_cvt_f64_high(b0), s1); + s0 = v_add(s0, v_mul(v_cvt_f64(a0), v_cvt_f64(b0))); + s1 = v_add(s1, v_mul(v_cvt_f64_high(a0), v_cvt_f64_high(b0))); } *result += v_reduce_sum(v_add(s0, s1)); + for( ; k < n; k++ ) + *result += (double)a[k]*(double)b[k]; vx_cleanup(); - return k; } -template<> inline int VBLAS::givensD(float* a, float* b, int n, float c, float s, +template<> inline void VBLAS::givensD(float* a, float* b, int n, float c, float s, double* na, double* nb) const { - if( n < VTraits::vlanes() ) - return 0; - int k = 0; + if( n <= 0 ) + return; + const int vl = VTraits::vlanes(); v_float32 c4 = vx_setall_f32(c), s4 = vx_setall_f32(s); v_float32 ns4 = vx_setall_f32(-s); v_float64 a0d = vx_setzero_f64(), a1d = vx_setzero_f64(); v_float64 b0d = vx_setzero_f64(), b1d = vx_setzero_f64(); - for( ; k <= n - VTraits::vlanes(); k += VTraits::vlanes() ) + int k = 0; + for( ; k <= n - vl; k += vl ) { v_float32 a0 = vx_load(a + k); v_float32 b0 = vx_load(b + k); - v_float32 t0 = v_fma(a0, c4, v_mul(b0, s4)); - v_float32 t1 = v_fma(a0, ns4, v_mul(b0, c4)); + v_float32 t0 = v_add(v_mul(a0, c4), v_mul(b0, s4)); + v_float32 t1 = v_add(v_mul(a0, ns4), v_mul(b0, c4)); v_store(a + k, t0); v_store(b + k, t1); - // accumulate squared norms in double precision v_float64 t0lo = v_cvt_f64(t0), t0hi = v_cvt_f64_high(t0); v_float64 t1lo = v_cvt_f64(t1), t1hi = v_cvt_f64_high(t1); - a0d = v_fma(t0lo, t0lo, a0d); - a1d = v_fma(t0hi, t0hi, a1d); - b0d = v_fma(t1lo, t1lo, b0d); - b1d = v_fma(t1hi, t1hi, b1d); + a0d = v_add(a0d, v_mul(t0lo, t0lo)); + a1d = v_add(a1d, v_mul(t0hi, t0hi)); + b0d = v_add(b0d, v_mul(t1lo, t1lo)); + b1d = v_add(b1d, v_mul(t1hi, t1hi)); } *na += v_reduce_sum(v_add(a0d, a1d)); *nb += v_reduce_sum(v_add(b0d, b1d)); + for( ; k < n; k++ ) + { + float t0 = c*a[k] + s*b[k]; + float t1 = -s*a[k] + c*b[k]; + a[k] = t0; b[k] = t1; + *na += (double)t0*t0; + *nb += (double)t1*t1; + } vx_cleanup(); - return k; } -#endif //CV_SIMD_64F for float + + +template<> inline void VBLAS::dotD(const double* a, const double* b, int n, double* result) const +{ + if( n <= 0 ) + return; + const int vl = VTraits::vlanes(); + v_float64 s0 = vx_setzero_f64(); + int k = 0; + for( ; k <= n - vl; k += vl ) + { + v_float64 a0 = vx_load(a + k); + v_float64 b0 = vx_load(b + k); + s0 = v_add(s0, v_mul(a0, b0)); + } + *result += v_reduce_sum(s0); + for( ; k < n; k++ ) + *result += a[k]*b[k]; + vx_cleanup(); +} + + +template<> inline void VBLAS::givensD(double* a, double* b, int n, double c, double s, + double* na, double* nb) const +{ + if( n <= 0 ) + return; + const int vl = VTraits::vlanes(); + v_float64 c2 = vx_setall_f64(c), s2 = vx_setall_f64(s); + v_float64 ns2 = vx_setall_f64(-s); + v_float64 nacc = vx_setzero_f64(), nbcc = vx_setzero_f64(); + int k = 0; + for( ; k <= n - vl; k += vl ) + { + v_float64 a0 = vx_load(a + k); + v_float64 b0 = vx_load(b + k); + v_float64 t0 = v_add(v_mul(a0, c2), v_mul(b0, s2)); + v_float64 t1 = v_add(v_mul(a0, ns2), v_mul(b0, c2)); + v_store(a + k, t0); + v_store(b + k, t1); + nacc = v_add(nacc, v_mul(t0, t0)); + nbcc = v_add(nbcc, v_mul(t1, t1)); + } + *na += v_reduce_sum(nacc); + *nb += v_reduce_sum(nbcc); + for( ; k < n; k++ ) + { + double t0 = c*a[k] + s*b[k]; + double t1 = -s*a[k] + c*b[k]; + a[k] = t0; b[k] = t1; + *na += t0*t0; + *nb += t1*t1; + } + vx_cleanup(); +} +#endif // CV_SIMD_64F #if (CV_SIMD_64F || CV_SIMD_SCALABLE_64F) @@ -196,72 +266,96 @@ template<> inline int VBLAS::dot(const double* a, const double* b, int n } -template<> inline int VBLAS::givens(double* a, double* b, int n, double c, double s) const +template<> inline void VBLAS::givens(double* a, double* b, int n, double c, double s) const { - int k = 0; + if( n <= 0 ) + return; + const int vl = VTraits::vlanes(); v_float64 c2 = vx_setall_f64(c), s2 = vx_setall_f64(s); v_float64 ns2 = vx_setall_f64(-s); - for( ; k <= n - VTraits::vlanes(); k += VTraits::vlanes() ) + int k = 0; + for( ; k <= n - vl; k += vl ) { v_float64 a0 = vx_load(a + k); v_float64 b0 = vx_load(b + k); - v_float64 t0 = v_fma(a0, c2, v_mul(b0, s2)); - v_float64 t1 = v_fma(a0, ns2, v_mul(b0, c2)); + v_float64 t0 = v_add(v_mul(a0, c2), v_mul(b0, s2)); + v_float64 t1 = v_add(v_mul(a0, ns2), v_mul(b0, c2)); v_store(a + k, t0); v_store(b + k, t1); } - vx_cleanup(); - return k; -} - - -template<> inline int VBLAS::dotD(const double* a, const double* b, int n, double* result) const -{ - if( n < 2*VTraits::vlanes() ) - return 0; - int k = 0; - v_float64 s0 = vx_setzero_f64(); - for( ; k <= n - VTraits::vlanes(); k += VTraits::vlanes() ) + for( ; k < n; k++ ) { - v_float64 a0 = vx_load(a + k); - v_float64 b0 = vx_load(b + k); - s0 = v_fma(a0, b0, s0); + double t0 = c*a[k] + s*b[k]; + double t1 = -s*a[k] + c*b[k]; + a[k] = t0; b[k] = t1; } - *result += v_reduce_sum(s0); vx_cleanup(); - return k; } +#endif // CV_SIMD_64F +#endif // CV_SIMD -template<> inline int VBLAS::givensD(double* a, double* b, int n, double c, double s, +#if !CV_SIMD +template<> inline void VBLAS::givens(float* a, float* b, int n, float c, float s) const +{ + for( int k = 0; k < n; k++ ) + { + float t0 = c*a[k] + s*b[k]; + float t1 = -s*a[k] + c*b[k]; + a[k] = t0; b[k] = t1; + } +} +#endif + +#if !(CV_SIMD_64F || CV_SIMD_SCALABLE_64F) +template<> inline void VBLAS::dotD(const float* a, const float* b, int n, double* result) const +{ + for( int k = 0; k < n; k++ ) + *result += (double)a[k]*(double)b[k]; +} + +template<> inline void VBLAS::givensD(float* a, float* b, int n, float c, float s, + double* na, double* nb) const +{ + for( int k = 0; k < n; k++ ) + { + float t0 = c*a[k] + s*b[k]; + float t1 = -s*a[k] + c*b[k]; + a[k] = t0; b[k] = t1; + *na += (double)t0*t0; + *nb += (double)t1*t1; + } +} + +template<> inline void VBLAS::dotD(const double* a, const double* b, int n, double* result) const +{ + for( int k = 0; k < n; k++ ) + *result += a[k]*b[k]; +} + +template<> inline void VBLAS::givens(double* a, double* b, int n, double c, double s) const +{ + for( int k = 0; k < n; k++ ) + { + double t0 = c*a[k] + s*b[k]; + double t1 = -s*a[k] + c*b[k]; + a[k] = t0; b[k] = t1; + } +} + +template<> inline void VBLAS::givensD(double* a, double* b, int n, double c, double s, double* na, double* nb) const { - if( n < VTraits::vlanes() ) - return 0; - int k = 0; - v_float64 c2 = vx_setall_f64(c), s2 = vx_setall_f64(s); - v_float64 ns2 = vx_setall_f64(-s); - v_float64 nacc = vx_setzero_f64(), nbcc = vx_setzero_f64(); - for( ; k <= n - VTraits::vlanes(); k += VTraits::vlanes() ) + for( int k = 0; k < n; k++ ) { - v_float64 a0 = vx_load(a + k); - v_float64 b0 = vx_load(b + k); - v_float64 t0 = v_fma(a0, c2, v_mul(b0, s2)); - v_float64 t1 = v_fma(a0, ns2, v_mul(b0, c2)); - v_store(a + k, t0); - v_store(b + k, t1); - nacc = v_fma(t0, t0, nacc); - nbcc = v_fma(t1, t1, nbcc); + double t0 = c*a[k] + s*b[k]; + double t1 = -s*a[k] + c*b[k]; + a[k] = t0; b[k] = t1; + *na += t0*t0; + *nb += t1*t1; } - *na += v_reduce_sum(nacc); - *nb += v_reduce_sum(nbcc); - vx_cleanup(); - return k; } - - -#endif //CV_SIMD_64F -#endif //CV_SIMD +#endif template void JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, @@ -280,12 +374,7 @@ JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, { _Tp* Ai = At + i*astep; sd = 0; - k = vblas.dotD(Ai, Ai, m, &sd); - for( ; k < m; k++ ) - { - _Tp t = Ai[k]; - sd += (double)t*t; - } + vblas.dotD(Ai, Ai, m, &sd); W[i] = sd; if( Vt ) @@ -306,9 +395,7 @@ JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, _Tp *Ai = At + i*astep, *Aj = At + j*astep; double a = W[i], p = 0, b = W[j]; - k = vblas.dotD(Ai, Aj, m, &p); - for( ; k < m; k++ ) - p += (double)Ai[k]*Aj[k]; + vblas.dotD(Ai, Aj, m, &p); if( std::abs(p) <= eps*std::sqrt((double)a*b) ) continue; @@ -328,15 +415,7 @@ JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, } a = b = 0; - k = vblas.givensD(Ai, Aj, m, c, s, &a, &b); - for( ; k < m; k++ ) - { - _Tp t0 = c*Ai[k] + s*Aj[k]; - _Tp t1 = -s*Ai[k] + c*Aj[k]; - Ai[k] = t0; Aj[k] = t1; - - a += (double)t0*t0; b += (double)t1*t1; - } + vblas.givensD(Ai, Aj, m, c, s, &a, &b); W[i] = a; W[j] = b; changed = true; @@ -344,14 +423,7 @@ JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, if( Vt ) { _Tp *Vi = Vt + i*vstep, *Vj = Vt + j*vstep; - k = vblas.givens(Vi, Vj, n, c, s); - - for( ; k < n; k++ ) - { - _Tp t0 = c*Vi[k] + s*Vj[k]; - _Tp t1 = -s*Vi[k] + c*Vj[k]; - Vi[k] = t0; Vj[k] = t1; - } + vblas.givens(Vi, Vj, n, c, s); } } if( !changed ) @@ -362,12 +434,7 @@ JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, { _Tp* Ai = At + i*astep; sd = 0; - k = vblas.dotD(Ai, Ai, m, &sd); - for( ; k < m; k++ ) - { - _Tp t = Ai[k]; - sd += (double)t*t; - } + vblas.dotD(Ai, Ai, m, &sd); W[i] = std::sqrt(sd); } diff --git a/modules/photo/test/test_hdr.cpp b/modules/photo/test/test_hdr.cpp index 8f11a4729a..9922a01fff 100644 --- a/modules/photo/test/test_hdr.cpp +++ b/modules/photo/test/test_hdr.cpp @@ -286,7 +286,11 @@ TEST(Photo_CalibrateDebevec, regression) diff = diff.mul(1.0f / response); double max; minMaxLoc(diff, NULL, &max); +#if defined(__arm__) || defined(__aarch64__) ASSERT_LT(max, 0.25); +#else + ASSERT_LT(max, 0.15); +#endif } TEST(Photo_CalibrateRobertson, regression)