#endif
template<typename _Tp> void
-JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, int m, int n, int n1)
+JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep, int m, int n, int n1, double minval)
{
VBLAS<_Tp> vblas;
AutoBuffer<double> Wbuf(n);
delta = p*p*0.5/(gamma + beta);
}
- if( iter % 2 )
+ W[i] += delta;
+ W[j] -= delta;
+
+ if( iter % 2 != 0 && W[i] > 0 && W[j] > 0 )
{
- W[i] += delta;
- W[j] -= delta;
-
k = vblas.givens(Ai, Aj, m, c, s);
for( ; k < m; k++ )
if( !Vt )
return;
+
RNG rng(0x12345678);
for( i = 0; i < n1; i++ )
{
sd = i < n ? W[i] : 0;
- while( sd == 0 )
+ while( sd <= minval )
{
// if we got a zero singular value, then in order to get the corresponding left singular vector
// we generate a random vector, project it to the previously computed left singular vectors,
static void JacobiSVD(float* At, size_t astep, float* W, float* Vt, size_t vstep, int m, int n, int n1=-1)
{
- JacobiSVDImpl_(At, astep, W, Vt, vstep, m, n, !Vt ? 0 : n1 < 0 ? n : n1);
+ JacobiSVDImpl_(At, astep, W, Vt, vstep, m, n, !Vt ? 0 : n1 < 0 ? n : n1, FLT_MIN);
}
static void JacobiSVD(double* At, size_t astep, double* W, double* Vt, size_t vstep, int m, int n, int n1=-1)
{
- JacobiSVDImpl_(At, astep, W, Vt, vstep, m, n, !Vt ? 0 : n1 < 0 ? n : n1);
+ JacobiSVDImpl_(At, astep, W, Vt, vstep, m, n, !Vt ? 0 : n1 < 0 ? n : n1, DBL_MIN);
}
/* y[0:m,0:n] += diag(a[0:1,0:m]) * x[0:m,0:n] */
SVD svd(A, SVD::FULL_UV);
if( norm(A*svd.vt.row(3).t(), CV_C) > FLT_EPSILON )
throw test_excep();
+
+ Mat Dp(3,3,CV_32FC1);
+ Mat Dc(3,3,CV_32FC1);
+ Mat Q(3,3,CV_32FC1);
+ Mat U,Vt,R,T,W;
+
+ Dp.at<float>(0,0)=0.86483884; Dp.at<float>(0,1)= -0.3077251; Dp.at<float>(0,2)=-0.55711365;
+ Dp.at<float>(1,0)=0.49294353; Dp.at<float>(1,1)=-0.24209651; Dp.at<float>(1,2)=-0.25084701;
+ Dp.at<float>(2,0)=0; Dp.at<float>(2,1)=0; Dp.at<float>(2,2)=0;
+
+ Dc.at<float>(0,0)=0.75632739; Dc.at<float>(0,1)= -0.38859656; Dc.at<float>(0,2)=-0.36773083;
+ Dc.at<float>(1,0)=0.9699229; Dc.at<float>(1,1)=-0.49858192; Dc.at<float>(1,2)=-0.47134098;
+ Dc.at<float>(2,0)=0.10566688; Dc.at<float>(2,1)=-0.060333252; Dc.at<float>(2,2)=-0.045333147;
+
+ Q=Dp*Dc.t();
+ SVD decomp;
+ decomp=SVD(Q);
+ U=decomp.u;
+ Vt=decomp.vt;
+ W=decomp.w;
+ Mat I = Mat::eye(3, 3, CV_32F);
+
+ if( norm(U*U.t(), I, CV_C) > FLT_EPSILON ||
+ norm(Vt*Vt.t(), I, CV_C) > FLT_EPSILON ||
+ W.at<float>(2) < 0 || W.at<float>(1) < W.at<float>(2) ||
+ W.at<float>(0) < W.at<float>(1) ||
+ norm(U*Mat::diag(W)*Vt, Q, CV_C) > FLT_EPSILON )
+ throw test_excep();
}
catch(const test_excep&)
{