in modules/core/src/lapack.cpp [1050:1321]
bool cv::solve( InputArray _src, InputArray _src2arg, OutputArray _dst, int method )
{
bool result = true;
Mat src = _src.getMat(), _src2 = _src2arg.getMat();
int type = src.type();
bool is_normal = (method & DECOMP_NORMAL) != 0;
CV_Assert( type == _src2.type() && (type == CV_32F || type == CV_64F) );
method &= ~DECOMP_NORMAL;
CV_Assert( (method != DECOMP_LU && method != DECOMP_CHOLESKY) ||
is_normal || src.rows == src.cols );
// check case of a single equation and small matrix
if( (method == DECOMP_LU || method == DECOMP_CHOLESKY) && !is_normal &&
src.rows <= 3 && src.rows == src.cols && _src2.cols == 1 )
{
_dst.create( src.cols, _src2.cols, src.type() );
Mat dst = _dst.getMat();
#define bf(y) ((float*)(bdata + y*src2step))[0]
#define bd(y) ((double*)(bdata + y*src2step))[0]
const uchar* srcdata = src.ptr();
const uchar* bdata = _src2.ptr();
uchar* dstdata = dst.ptr();
size_t srcstep = src.step;
size_t src2step = _src2.step;
size_t dststep = dst.step;
if( src.rows == 2 )
{
if( type == CV_32FC1 )
{
double d = det2(Sf);
if( d != 0. )
{
double t;
d = 1./d;
t = (float)(((double)bf(0)*Sf(1,1) - (double)bf(1)*Sf(0,1))*d);
Df(1,0) = (float)(((double)bf(1)*Sf(0,0) - (double)bf(0)*Sf(1,0))*d);
Df(0,0) = (float)t;
}
else
result = false;
}
else
{
double d = det2(Sd);
if( d != 0. )
{
double t;
d = 1./d;
t = (bd(0)*Sd(1,1) - bd(1)*Sd(0,1))*d;
Dd(1,0) = (bd(1)*Sd(0,0) - bd(0)*Sd(1,0))*d;
Dd(0,0) = t;
}
else
result = false;
}
}
else if( src.rows == 3 )
{
if( type == CV_32FC1 )
{
double d = det3(Sf);
if( d != 0. )
{
float t[3];
d = 1./d;
t[0] = (float)(d*
(bf(0)*((double)Sf(1,1)*Sf(2,2) - (double)Sf(1,2)*Sf(2,1)) -
Sf(0,1)*((double)bf(1)*Sf(2,2) - (double)Sf(1,2)*bf(2)) +
Sf(0,2)*((double)bf(1)*Sf(2,1) - (double)Sf(1,1)*bf(2))));
t[1] = (float)(d*
(Sf(0,0)*(double)(bf(1)*Sf(2,2) - (double)Sf(1,2)*bf(2)) -
bf(0)*((double)Sf(1,0)*Sf(2,2) - (double)Sf(1,2)*Sf(2,0)) +
Sf(0,2)*((double)Sf(1,0)*bf(2) - (double)bf(1)*Sf(2,0))));
t[2] = (float)(d*
(Sf(0,0)*((double)Sf(1,1)*bf(2) - (double)bf(1)*Sf(2,1)) -
Sf(0,1)*((double)Sf(1,0)*bf(2) - (double)bf(1)*Sf(2,0)) +
bf(0)*((double)Sf(1,0)*Sf(2,1) - (double)Sf(1,1)*Sf(2,0))));
Df(0,0) = t[0];
Df(1,0) = t[1];
Df(2,0) = t[2];
}
else
result = false;
}
else
{
double d = det3(Sd);
if( d != 0. )
{
double t[9];
d = 1./d;
t[0] = ((Sd(1,1) * Sd(2,2) - Sd(1,2) * Sd(2,1))*bd(0) +
(Sd(0,2) * Sd(2,1) - Sd(0,1) * Sd(2,2))*bd(1) +
(Sd(0,1) * Sd(1,2) - Sd(0,2) * Sd(1,1))*bd(2))*d;
t[1] = ((Sd(1,2) * Sd(2,0) - Sd(1,0) * Sd(2,2))*bd(0) +
(Sd(0,0) * Sd(2,2) - Sd(0,2) * Sd(2,0))*bd(1) +
(Sd(0,2) * Sd(1,0) - Sd(0,0) * Sd(1,2))*bd(2))*d;
t[2] = ((Sd(1,0) * Sd(2,1) - Sd(1,1) * Sd(2,0))*bd(0) +
(Sd(0,1) * Sd(2,0) - Sd(0,0) * Sd(2,1))*bd(1) +
(Sd(0,0) * Sd(1,1) - Sd(0,1) * Sd(1,0))*bd(2))*d;
Dd(0,0) = t[0];
Dd(1,0) = t[1];
Dd(2,0) = t[2];
}
else
result = false;
}
}
else
{
assert( src.rows == 1 );
if( type == CV_32FC1 )
{
double d = Sf(0,0);
if( d != 0. )
Df(0,0) = (float)(bf(0)/d);
else
result = false;
}
else
{
double d = Sd(0,0);
if( d != 0. )
Dd(0,0) = (bd(0)/d);
else
result = false;
}
}
return result;
}
if( method == DECOMP_QR )
method = DECOMP_SVD;
int m = src.rows, m_ = m, n = src.cols, nb = _src2.cols;
size_t esz = CV_ELEM_SIZE(type), bufsize = 0;
size_t vstep = alignSize(n*esz, 16);
size_t astep = method == DECOMP_SVD && !is_normal ? alignSize(m*esz, 16) : vstep;
AutoBuffer<uchar> buffer;
Mat src2 = _src2;
_dst.create( src.cols, src2.cols, src.type() );
Mat dst = _dst.getMat();
if( m < n )
CV_Error(CV_StsBadArg, "The function can not solve under-determined linear systems" );
if( m == n )
is_normal = false;
else if( is_normal )
{
m_ = n;
if( method == DECOMP_SVD )
method = DECOMP_EIG;
}
size_t asize = astep*(method == DECOMP_SVD || is_normal ? n : m);
bufsize += asize + 32;
if( is_normal )
bufsize += n*nb*esz;
if( method == DECOMP_SVD || method == DECOMP_EIG )
bufsize += n*5*esz + n*vstep + nb*sizeof(double) + 32;
buffer.allocate(bufsize);
uchar* ptr = alignPtr((uchar*)buffer, 16);
Mat a(m_, n, type, ptr, astep);
if( is_normal )
mulTransposed(src, a, true);
else if( method != DECOMP_SVD )
src.copyTo(a);
else
{
a = Mat(n, m_, type, ptr, astep);
transpose(src, a);
}
ptr += asize;
if( !is_normal )
{
if( method == DECOMP_LU || method == DECOMP_CHOLESKY )
src2.copyTo(dst);
}
else
{
// a'*b
if( method == DECOMP_LU || method == DECOMP_CHOLESKY )
gemm( src, src2, 1, Mat(), 0, dst, GEMM_1_T );
else
{
Mat tmp(n, nb, type, ptr);
ptr += n*nb*esz;
gemm( src, src2, 1, Mat(), 0, tmp, GEMM_1_T );
src2 = tmp;
}
}
if( method == DECOMP_LU )
{
if( type == CV_32F )
result = hal::LU(a.ptr<float>(), a.step, n, dst.ptr<float>(), dst.step, nb) != 0;
else
result = hal::LU(a.ptr<double>(), a.step, n, dst.ptr<double>(), dst.step, nb) != 0;
}
else if( method == DECOMP_CHOLESKY )
{
if( type == CV_32F )
result = hal::Cholesky(a.ptr<float>(), a.step, n, dst.ptr<float>(), dst.step, nb);
else
result = hal::Cholesky(a.ptr<double>(), a.step, n, dst.ptr<double>(), dst.step, nb);
}
else
{
ptr = alignPtr(ptr, 16);
Mat v(n, n, type, ptr, vstep), w(n, 1, type, ptr + vstep*n), u;
ptr += n*(vstep + esz);
if( method == DECOMP_EIG )
{
if( type == CV_32F )
Jacobi(a.ptr<float>(), a.step, w.ptr<float>(), v.ptr<float>(), v.step, n, ptr);
else
Jacobi(a.ptr<double>(), a.step, w.ptr<double>(), v.ptr<double>(), v.step, n, ptr);
u = v;
}
else
{
if( type == CV_32F )
JacobiSVD(a.ptr<float>(), a.step, w.ptr<float>(), v.ptr<float>(), v.step, m_, n);
else
JacobiSVD(a.ptr<double>(), a.step, w.ptr<double>(), v.ptr<double>(), v.step, m_, n);
u = a;
}
if( type == CV_32F )
{
SVBkSb(m_, n, w.ptr<float>(), 0, u.ptr<float>(), u.step, true,
v.ptr<float>(), v.step, true, src2.ptr<float>(),
src2.step, nb, dst.ptr<float>(), dst.step, ptr);
}
else
{
SVBkSb(m_, n, w.ptr<double>(), 0, u.ptr<double>(), u.step, true,
v.ptr<double>(), v.step, true, src2.ptr<double>(),
src2.step, nb, dst.ptr<double>(), dst.step, ptr);
}
result = true;
}
if( !result )
dst = Scalar(0);
return result;
}