GitHub Viewer
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "precomp.hpp"
#include
#ifdef HAVE_EIGEN
#include
#include
#include "opencv2/core/eigen.hpp"
#endif
#if defined _M_IX86 && defined _MSC_VER && _MSC_VER < 1700
#pragma float_control(precise, on)
#endif
namespace cv
{
int LU(float* A, size_t astep, int m, float* b, size_t bstep, int n)
{
CV_INSTRUMENT_REGION()
return hal::LU32f(A, astep, m, b, bstep, n);
}
int LU(double* A, size_t astep, int m, double* b, size_t bstep, int n)
{
CV_INSTRUMENT_REGION()
return hal::LU64f(A, astep, m, b, bstep, n);
}
bool Cholesky(float* A, size_t astep, int m, float* b, size_t bstep, int n)
{
CV_INSTRUMENT_REGION()
return hal::Cholesky32f(A, astep, m, b, bstep, n);
}
bool Cholesky(double* A, size_t astep, int m, double* b, size_t bstep, int n)
{
CV_INSTRUMENT_REGION()
return hal::Cholesky64f(A, astep, m, b, bstep, n);
}
template static inline _Tp hypot(_Tp a, _Tp b)
{
a = std::abs(a);
b = std::abs(b);
if( a > b )
{
b /= a;
return a*std::sqrt(1 + b*b);
}
if( b > 0 )
{
a /= b;
return b*std::sqrt(1 + a*a);
}
return 0;
}
template bool
JacobiImpl_( _Tp* A, size_t astep, _Tp* W, _Tp* V, size_t vstep, int n, uchar* buf )
{
const _Tp eps = std::numeric_limits::epsilon();
int i, j, k, m;
astep /= sizeof(A[0]);
if( V )
{
vstep /= sizeof(V[0]);
for( i = 0; i < n; i++ )
{
for( j = 0; j < n; j++ )
V[i*vstep + j] = (_Tp)0;
V[i*vstep + i] = (_Tp)1;
}
}
int iters, maxIters = n*n*30;
int* indR = (int*)alignPtr(buf, sizeof(int));
int* indC = indR + n;
_Tp mv = (_Tp)0;
for( k = 0; k < n; k++ )
{
W[k] = A[(astep + 1)*k];
if( k < n - 1 )
{
for( m = k+1, mv = std::abs(A[astep*k + m]), i = k+2; i < n; i++ )
{
_Tp val = std::abs(A[astep*k+i]);
if( mv < val )
mv = val, m = i;
}
indR[k] = m;
}
if( k > 0 )
{
for( m = 0, mv = std::abs(A[k]), i = 1; i < k; i++ )
{
_Tp val = std::abs(A[astep*i+k]);
if( mv < val )
mv = val, m = i;
}
indC[k] = m;
}
}
if( n > 1 ) for( iters = 0; iters < maxIters; iters++ )
{
// find index (k,l) of pivot p
for( k = 0, mv = std::abs(A[indR[0]]), i = 1; i < n-1; i++ )
{
_Tp val = std::abs(A[astep*i + indR[i]]);
if( mv < val )
mv = val, k = i;
}
int l = indR[k];
for( i = 1; i < n; i++ )
{
_Tp val = std::abs(A[astep*indC[i] + i]);
if( mv < val )
mv = val, k = indC[i], l = i;
}
_Tp p = A[astep*k + l];
if( std::abs(p) 0 )
{
for( m = 0, mv = std::abs(A[idx]), i = 1; i < idx; i++ )
{
_Tp val = std::abs(A[astep*i+idx]);
if( mv < val )
mv = val, m = i;
}
indC[idx] = m;
}
}
}
// sort eigenvalues & eigenvectors
for( k = 0; k < n-1; k++ )
{
m = k;
for( i = k+1; i < n; i++ )
{
if( W[m] < W[i] )
m = i;
}
if( k != m )
{
std::swap(W[m], W[k]);
if( V )
for( i = 0; i < n; i++ )
std::swap(V[vstep*m + i], V[vstep*k + i]);
}
}
return true;
}
static bool Jacobi( float* S, size_t sstep, float* e, float* E, size_t estep, int n, uchar* buf )
{
return JacobiImpl_(S, sstep, e, E, estep, n, buf);
}
static bool Jacobi( double* S, size_t sstep, double* e, double* E, size_t estep, int n, uchar* buf )
{
return JacobiImpl_(S, sstep, e, E, estep, n, buf);
}
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 givensx(T*, T*, int, T, T, T*, T*) const { return 0; }
};
#if CV_SIMD128
template inline int VBLAS::dot(const float* a, const float* b, int n, float* result) const
{
if( n < 8 )
return 0;
int k = 0;
v_float32x4 s0 = v_setzero_f32();
for( ; k