Commit 1a2a524a authored by Davis King's avatar Davis King

Changed LAPACK bindings so that really small matrices don't use LAPACK

for eigenvalue_decomposition or the triangular solver.
parent 2bbab691
......@@ -18,6 +18,8 @@
#include "lapack/syevr.h"
#endif
#define DLIB_LAPACK_EIGENVALUE_DECOMP_SIZE_THRESH 4
namespace dlib
{
......@@ -178,35 +180,41 @@ namespace dlib
V = A;
#ifdef DLIB_USE_LAPACK
e = 0;
if (A.nr() > DLIB_LAPACK_EIGENVALUE_DECOMP_SIZE_THRESH)
{
e = 0;
// We could compute the result using syev()
//lapack::syev('V', 'L', V, d);
// We could compute the result using syev()
//lapack::syev('V', 'L', V, d);
// Instead, we use syevr because its faster and maybe more stable.
matrix_type tempA(A);
matrix<lapack::integer,0,0,mem_manager_type,layout_type> isupz;
// Instead, we use syevr because its faster and maybe more stable.
matrix_type tempA(A);
matrix<lapack::integer,0,0,mem_manager_type,layout_type> isupz;
lapack::integer temp;
lapack::syevr('V','A','L',tempA,0,0,0,0,-1,temp,d,V,isupz);
#else
lapack::integer temp;
lapack::syevr('V','A','L',tempA,0,0,0,0,-1,temp,d,V,isupz);
}
#endif
// Tridiagonalize.
tred2();
// Diagonalize.
tql2();
#endif
}
else
{
#ifdef DLIB_USE_LAPACK
matrix<type,0,0,mem_manager_type, column_major_layout> temp, vl, vr;
temp = A;
lapack::geev('N', 'V', temp, d, e, vl, vr);
V = vr;
#else
if (A.nr() > DLIB_LAPACK_EIGENVALUE_DECOMP_SIZE_THRESH)
{
matrix<type,0,0,mem_manager_type, column_major_layout> temp, vl, vr;
temp = A;
lapack::geev('N', 'V', temp, d, e, vl, vr);
V = vr;
return;
}
#endif
H = A;
ort.set_size(n);
......@@ -216,7 +224,6 @@ namespace dlib
// Reduce Hessenberg to real Schur form.
hqr2();
#endif
}
}
......@@ -252,25 +259,27 @@ namespace dlib
V = A;
#ifdef DLIB_USE_LAPACK
e = 0;
// We could compute the result using syev()
//lapack::syev('V', 'L', V, d);
if (A.nr() > DLIB_LAPACK_EIGENVALUE_DECOMP_SIZE_THRESH)
{
e = 0;
// Instead, we use syevr because its faster and maybe more stable.
matrix_type tempA(A);
matrix<lapack::integer,0,0,mem_manager_type,layout_type> isupz;
// We could compute the result using syev()
//lapack::syev('V', 'L', V, d);
lapack::integer temp;
lapack::syevr('V','A','L',tempA,0,0,0,0,-1,temp,d,V,isupz);
// Instead, we use syevr because its faster and maybe more stable.
matrix_type tempA(A);
matrix<lapack::integer,0,0,mem_manager_type,layout_type> isupz;
#else
lapack::integer temp;
lapack::syevr('V','A','L',tempA,0,0,0,0,-1,temp,d,V,isupz);
return;
}
#endif
// Tridiagonalize.
tred2();
// Diagonalize.
tql2();
#endif
}
......
......@@ -519,10 +519,13 @@ namespace dlib
float *B, const int ldb)
{
#ifdef DLIB_USE_BLAS
cblas_strsm(Order, Side, Uplo, TransA, Diag, M, N, alpha, A, lda, B, ldb);
#else
local_trsm(Order, Side, Uplo, TransA, Diag, M, N, alpha, A, lda, B, ldb);
if (M > 4)
{
cblas_strsm(Order, Side, Uplo, TransA, Diag, M, N, alpha, A, lda, B, ldb);
return;
}
#endif
local_trsm(Order, Side, Uplo, TransA, Diag, M, N, alpha, A, lda, B, ldb);
}
inline void cblas_trsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
......@@ -532,10 +535,13 @@ namespace dlib
double *B, const int ldb)
{
#ifdef DLIB_USE_BLAS
cblas_dtrsm(Order, Side, Uplo, TransA, Diag, M, N, alpha, A, lda, B, ldb);
#else
local_trsm(Order, Side, Uplo, TransA, Diag, M, N, alpha, A, lda, B, ldb);
if (M > 4)
{
cblas_dtrsm(Order, Side, Uplo, TransA, Diag, M, N, alpha, A, lda, B, ldb);
return;
}
#endif
local_trsm(Order, Side, Uplo, TransA, Diag, M, N, alpha, A, lda, B, ldb);
}
inline void cblas_trsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment