mirror of
https://github.com/opencv/opencv.git
synced 2026-07-29 15:23:05 +04:00
Merge pull request #18571 from vpisarev:add_lapack
Added clapack * bring a small subset of Lapack, automatically converted to C, into OpenCV * added missing lsame_ prototype * * small fix in make_clapack script * trying to fix remaining CI problems * fixed character arrays' initializers * get rid of F2C_STR_MAX * * added back single-precision versions for QR, LU and Cholesky decompositions. It adds very little extra overhead. * added stub version of sdesdd. * uncommented calls to all the single-precision Lapack functions from opencv/core/src/hal_internal.cpp. * fixed warning from Visual Studio + cleaned f2c runtime a bit * * regenerated Lapack w/o forward declarations of intrinsic functions (such as sqrt(), r_cnjg() etc.) * at once, trailing whitespaces are removed from the generated sources, just in case * since there is no declarations of intrinsic functions anymore, we could turn some of them into inline functions * trying to eliminate the crash on ARM * fixed API and semantics of s_copy * * CLapack has been tested successfully. It's now time to restore the standard LAPACK detection procedure * removed some more trailing whitespaces * * retained only the essential stuff in CLapack * added checks to lapack calls to gracefully return "not implemented" instead of returning invalid results with "ok" status * disabled warning when building lapack * cmake: update LAPACK detection Co-authored-by: Alexander Alekhin <alexander.a.alekhin@gmail.com>
This commit is contained in:
@@ -100,7 +100,7 @@ template <typename fptype> static inline int
|
||||
lapack_LU(fptype* a, size_t a_step, int m, fptype* b, size_t b_step, int n, int* info)
|
||||
{
|
||||
int lda = (int)(a_step / sizeof(fptype)), sign = 0;
|
||||
int* piv = new int[m];
|
||||
std::vector<int> piv(m+1);
|
||||
|
||||
transpose_square_inplace(a, lda, m);
|
||||
|
||||
@@ -109,34 +109,35 @@ lapack_LU(fptype* a, size_t a_step, int m, fptype* b, size_t b_step, int n, int*
|
||||
if(n == 1 && b_step == sizeof(fptype))
|
||||
{
|
||||
if(typeid(fptype) == typeid(float))
|
||||
sgesv_(&m, &n, (float*)a, &lda, piv, (float*)b, &m, info);
|
||||
sgesv_(&m, &n, (float*)a, &lda, &piv[0], (float*)b, &m, info);
|
||||
else if(typeid(fptype) == typeid(double))
|
||||
dgesv_(&m, &n, (double*)a, &lda, piv, (double*)b, &m, info);
|
||||
dgesv_(&m, &n, (double*)a, &lda, &piv[0], (double*)b, &m, info);
|
||||
}
|
||||
else
|
||||
{
|
||||
int ldb = (int)(b_step / sizeof(fptype));
|
||||
fptype* tmpB = new fptype[m*n];
|
||||
std::vector<fptype> tmpB(m*n+1);
|
||||
|
||||
transpose(b, ldb, tmpB, m, m, n);
|
||||
transpose(b, ldb, &tmpB[0], m, m, n);
|
||||
|
||||
if(typeid(fptype) == typeid(float))
|
||||
sgesv_(&m, &n, (float*)a, &lda, piv, (float*)tmpB, &m, info);
|
||||
sgesv_(&m, &n, (float*)a, &lda, &piv[0], (float*)&tmpB[0], &m, info);
|
||||
else if(typeid(fptype) == typeid(double))
|
||||
dgesv_(&m, &n, (double*)a, &lda, piv, (double*)tmpB, &m, info);
|
||||
dgesv_(&m, &n, (double*)a, &lda, &piv[0], (double*)&tmpB[0], &m, info);
|
||||
|
||||
transpose(tmpB, m, b, ldb, n, m);
|
||||
delete[] tmpB;
|
||||
transpose(&tmpB[0], m, b, ldb, n, m);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if(typeid(fptype) == typeid(float))
|
||||
sgetrf_(&m, &m, (float*)a, &lda, piv, info);
|
||||
sgetrf_(&m, &m, (float*)a, &lda, &piv[0], info);
|
||||
else if(typeid(fptype) == typeid(double))
|
||||
dgetrf_(&m, &m, (double*)a, &lda, piv, info);
|
||||
dgetrf_(&m, &m, (double*)a, &lda, &piv[0], info);
|
||||
}
|
||||
|
||||
int retcode = *info >= 0 ? CV_HAL_ERROR_OK : CV_HAL_ERROR_NOT_IMPLEMENTED;
|
||||
|
||||
if(*info == 0)
|
||||
{
|
||||
for(int i = 0; i < m; i++)
|
||||
@@ -146,8 +147,7 @@ lapack_LU(fptype* a, size_t a_step, int m, fptype* b, size_t b_step, int n, int*
|
||||
else
|
||||
*info = 0; //in opencv LU function zero means error
|
||||
|
||||
delete[] piv;
|
||||
return CV_HAL_ERROR_OK;
|
||||
return retcode;
|
||||
}
|
||||
|
||||
template <typename fptype> static inline int
|
||||
@@ -192,7 +192,7 @@ lapack_Cholesky(fptype* a, size_t a_step, int m, fptype* b, size_t b_step, int n
|
||||
if(lapackStatus == 0) *info = true;
|
||||
else *info = false; //in opencv Cholesky function false means error
|
||||
|
||||
return CV_HAL_ERROR_OK;
|
||||
return lapackStatus >= 0 ? CV_HAL_ERROR_OK : CV_HAL_ERROR_NOT_IMPLEMENTED;
|
||||
}
|
||||
|
||||
template <typename fptype> static inline int
|
||||
@@ -202,7 +202,8 @@ lapack_SVD(fptype* a, size_t a_step, fptype *w, fptype* u, size_t u_step, fptype
|
||||
int ldv = (int)(v_step / sizeof(fptype));
|
||||
int ldu = (int)(u_step / sizeof(fptype));
|
||||
int lwork = -1;
|
||||
int* iworkBuf = new int[8*std::min(m, n)];
|
||||
std::vector<int> iworkBuf(8*std::min(m, n));
|
||||
std::vector<fptype> ubuf;
|
||||
fptype work1 = 0;
|
||||
|
||||
//A already transposed and m>=n
|
||||
@@ -221,22 +222,30 @@ lapack_SVD(fptype* a, size_t a_step, fptype *w, fptype* u, size_t u_step, fptype
|
||||
|
||||
if((flags & CV_HAL_SVD_MODIFY_A) && (flags & CV_HAL_SVD_FULL_UV)) //U stored in a
|
||||
{
|
||||
u = new fptype[m*m];
|
||||
ubuf.resize(m*m);
|
||||
u = &ubuf[0];
|
||||
ldu = m;
|
||||
}
|
||||
|
||||
if(typeid(fptype) == typeid(float))
|
||||
sgesdd_(mode, &m, &n, (float*)a, &lda, (float*)w, (float*)u, &ldu, (float*)vt, &ldv, (float*)&work1, &lwork, iworkBuf, info);
|
||||
sgesdd_(mode, &m, &n, (float*)a, &lda, (float*)w, (float*)u, &ldu,
|
||||
(float*)vt, &ldv, (float*)&work1, &lwork, &iworkBuf[0], info);
|
||||
else if(typeid(fptype) == typeid(double))
|
||||
dgesdd_(mode, &m, &n, (double*)a, &lda, (double*)w, (double*)u, &ldu, (double*)vt, &ldv, (double*)&work1, &lwork, iworkBuf, info);
|
||||
dgesdd_(mode, &m, &n, (double*)a, &lda, (double*)w, (double*)u, &ldu,
|
||||
(double*)vt, &ldv, (double*)&work1, &lwork, &iworkBuf[0], info);
|
||||
|
||||
if(*info < 0)
|
||||
return CV_HAL_ERROR_NOT_IMPLEMENTED;
|
||||
|
||||
lwork = (int)round(work1); //optimal buffer size
|
||||
fptype* buffer = new fptype[lwork + 1];
|
||||
std::vector<fptype> buffer(lwork + 1);
|
||||
|
||||
if(typeid(fptype) == typeid(float))
|
||||
sgesdd_(mode, &m, &n, (float*)a, &lda, (float*)w, (float*)u, &ldu, (float*)vt, &ldv, (float*)buffer, &lwork, iworkBuf, info);
|
||||
sgesdd_(mode, &m, &n, (float*)a, &lda, (float*)w, (float*)u, &ldu,
|
||||
(float*)vt, &ldv, (float*)&buffer[0], &lwork, &iworkBuf[0], info);
|
||||
else if(typeid(fptype) == typeid(double))
|
||||
dgesdd_(mode, &m, &n, (double*)a, &lda, (double*)w, (double*)u, &ldu, (double*)vt, &ldv, (double*)buffer, &lwork, iworkBuf, info);
|
||||
dgesdd_(mode, &m, &n, (double*)a, &lda, (double*)w, (double*)u, &ldu,
|
||||
(double*)vt, &ldv, (double*)&buffer[0], &lwork, &iworkBuf[0], info);
|
||||
|
||||
if(!(flags & CV_HAL_SVD_NO_UV))
|
||||
transpose_square_inplace(vt, ldv, n);
|
||||
@@ -246,11 +255,10 @@ lapack_SVD(fptype* a, size_t a_step, fptype *w, fptype* u, size_t u_step, fptype
|
||||
for(int i = 0; i < m; i++)
|
||||
for(int j = 0; j < m; j++)
|
||||
a[i*lda + j] = u[i*m + j];
|
||||
delete[] u;
|
||||
}
|
||||
|
||||
delete[] iworkBuf;
|
||||
delete[] buffer;
|
||||
if(*info < 0)
|
||||
return CV_HAL_ERROR_NOT_IMPLEMENTED;
|
||||
return CV_HAL_ERROR_OK;
|
||||
}
|
||||
|
||||
@@ -291,6 +299,9 @@ lapack_QR(fptype* a, size_t a_step, int m, int n, int k, fptype* b, size_t b_ste
|
||||
else if (typeid(fptype) == typeid(double))
|
||||
dgels_(mode, &m, &n, &k, (double*)tmpA, &ldtmpA, (double*)b, &m, (double*)&work1, &lwork, info);
|
||||
|
||||
if (*info < 0)
|
||||
return CV_HAL_ERROR_NOT_IMPLEMENTED;
|
||||
|
||||
lwork = cvRound(work1); //optimal buffer size
|
||||
std::vector<fptype> workBufMemHolder(lwork + 1);
|
||||
fptype* buffer = &workBufMemHolder.front();
|
||||
@@ -312,6 +323,9 @@ lapack_QR(fptype* a, size_t a_step, int m, int n, int k, fptype* b, size_t b_ste
|
||||
else if (typeid(fptype) == typeid(double))
|
||||
dgels_(mode, &m, &n, &k, (double*)tmpA, &ldtmpA, (double*)tmpB, &m, (double*)&work1, &lwork, info);
|
||||
|
||||
if (*info < 0)
|
||||
return CV_HAL_ERROR_NOT_IMPLEMENTED;
|
||||
|
||||
lwork = cvRound(work1); //optimal buffer size
|
||||
std::vector<fptype> workBufMemHolder(lwork + 1);
|
||||
fptype* buffer = &workBufMemHolder.front();
|
||||
@@ -331,6 +345,9 @@ lapack_QR(fptype* a, size_t a_step, int m, int n, int k, fptype* b, size_t b_ste
|
||||
else if (typeid(fptype) == typeid(double))
|
||||
dgeqrf_(&m, &n, (double*)tmpA, &ldtmpA, (double*)dst, (double*)&work1, &lwork, info);
|
||||
|
||||
if (*info < 0)
|
||||
return CV_HAL_ERROR_NOT_IMPLEMENTED;
|
||||
|
||||
lwork = cvRound(work1); //optimal buffer size
|
||||
std::vector<fptype> workBufMemHolder(lwork + 1);
|
||||
fptype* buffer = &workBufMemHolder.front();
|
||||
@@ -538,19 +555,17 @@ int lapack_Cholesky64f(double* a, size_t a_step, int m, double* b, size_t b_step
|
||||
|
||||
int lapack_SVD32f(float* a, size_t a_step, float *w, float* u, size_t u_step, float* vt, size_t v_step, int m, int n, int flags)
|
||||
{
|
||||
|
||||
if(m < HAL_SVD_SMALL_MATRIX_THRESH)
|
||||
if(m < HAL_SVD_SMALL_MATRIX_THRESH || n <= 0)
|
||||
return CV_HAL_ERROR_NOT_IMPLEMENTED;
|
||||
int info;
|
||||
int info = 0;
|
||||
return lapack_SVD(a, a_step, w, u, u_step, vt, v_step, m, n, flags, &info);
|
||||
}
|
||||
|
||||
int lapack_SVD64f(double* a, size_t a_step, double *w, double* u, size_t u_step, double* vt, size_t v_step, int m, int n, int flags)
|
||||
{
|
||||
|
||||
if(m < HAL_SVD_SMALL_MATRIX_THRESH)
|
||||
if(m < HAL_SVD_SMALL_MATRIX_THRESH || n <= 0)
|
||||
return CV_HAL_ERROR_NOT_IMPLEMENTED;
|
||||
int info;
|
||||
int info = 0;
|
||||
return lapack_SVD(a, a_step, w, u, u_step, vt, v_step, m, n, flags, &info);
|
||||
}
|
||||
|
||||
|
||||
Reference in New Issue
Block a user