From 235a678458f721e477ede5e370c759cf7f32f927 Mon Sep 17 00:00:00 2001 From: Andrey Kamaev Date: Thu, 4 Apr 2013 13:55:36 +0400 Subject: [PATCH] SVD: always update W vector for better algorithm convergency --- modules/core/src/lapack.cpp | 36 +++++++++--------------------------- 1 file changed, 9 insertions(+), 27 deletions(-) diff --git a/modules/core/src/lapack.cpp b/modules/core/src/lapack.cpp index b20273f6eb..f5fe53ae9e 100644 --- a/modules/core/src/lapack.cpp +++ b/modules/core/src/lapack.cpp @@ -577,10 +577,10 @@ JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, continue; p *= 2; - double beta = a - b, gamma = hypot((double)p, beta), delta; + double beta = a - b, gamma = hypot((double)p, beta); if( beta < 0 ) { - delta = (gamma - beta)*0.5; + double delta = (gamma - beta)*0.5; s = (_Tp)std::sqrt(delta/gamma); c = (_Tp)(p/(gamma*s*2)); } @@ -588,36 +588,18 @@ JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, { c = (_Tp)std::sqrt((gamma + beta)/(gamma*2)); s = (_Tp)(p/(gamma*c*2)); - delta = p*p*0.5/(gamma + beta); } - W[i] += delta; - W[j] -= delta; - - if( iter % 2 != 0 && W[i] > 0 && W[j] > 0 ) - { - k = vblas.givens(Ai, Aj, m, c, s); - - 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; - } - } - else + a = b = 0; + for( k = 0; k < m; k++ ) { - a = b = 0; - for( k = 0; 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; + _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; - } - W[i] = a; W[j] = b; + a += (double)t0*t0; b += (double)t1*t1; } + W[i] = a; W[j] = b; changed = true;