[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Thu, 21 Sep 2017 06:34:20 -0400 (EDT) |
branch: devel-yves-fix-matlab-pb
commit c4f955057c37fe0c9771c7b37b4a05bf60e4d9aa
Author: Yves Renard <address@hidden>
Date: Thu Sep 21 10:28:01 2017 +0200
Fix for compatibility with 64 bit integer BLAS and LAPACK versions
---
src/bgeot_geometric_trans.cc | 20 +-
src/bgeot_sparse_tensors.cc | 8 +-
src/getfem/bgeot_geometric_trans.h | 2 +-
src/getfem_assembling_tensors.cc | 9 +-
src/getfem_contact_and_friction_large_sliding.cc | 2 +-
src/getfem_generic_assembly.cc | 8 +-
src/getfem_mat_elem.cc | 11 +-
src/getfem_mesher.cc | 2 +-
src/gmm/gmm_blas_interface.h | 72 +++----
src/gmm/gmm_dense_lu.h | 6 +-
src/gmm/gmm_lapack_interface.h | 244 +++++++++++------------
src/gmm/gmm_opt.h | 2 +-
tests/test_gmm_matrix_functions.cc | 26 +--
13 files changed, 200 insertions(+), 212 deletions(-)
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index 7db4bd4..e260003 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -43,8 +43,8 @@ namespace bgeot {
return v;
}
- std::vector<int>& __ipvt_aux(){
- DEFINE_STATIC_THREAD_LOCAL(std::vector<int>, vi);
+ std::vector<long>& __ipvt_aux(){
+ DEFINE_STATIC_THREAD_LOCAL(std::vector<long>, vi);
return vi;
}
@@ -111,7 +111,7 @@ namespace bgeot {
// Optimized lu_factor for small square matrices
- size_type lu_factor(scalar_type *A, std::vector<int> &ipvt,
+ size_type lu_factor(scalar_type *A, std::vector<long> &ipvt,
size_type N) {
size_type info(0), i, j, jp, N_1 = N-1;
@@ -123,7 +123,7 @@ namespace bgeot {
scalar_type ap = gmm::abs(*(++it));
if (ap > max) { jp = i; max = ap; }
}
- ipvt[j] = int(jp + 1);
+ ipvt[j] = long(jp + 1);
if (max == scalar_type(0)) { info = j + 1; break; }
if (jp != j) {
@@ -142,7 +142,7 @@ namespace bgeot {
}
}
- ipvt[N-1] = int(N);
+ ipvt[N-1] = long(N);
}
return info;
}
@@ -170,20 +170,20 @@ namespace bgeot {
}
}
- static void lu_solve(const scalar_type *LU, const std::vector<int> &ipvt,
+ static void lu_solve(const scalar_type *LU, const std::vector<long> &ipvt,
scalar_type *x, scalar_type *b, int N) {
std::copy(b, b+N, x);
for(int i = 0; i < N; ++i)
- { int perm = ipvt[i]-1; if(i != perm) std::swap(x[i], x[perm]); }
+ { long perm = ipvt[i]-1; if(i != perm) std::swap(x[i], x[perm]); }
bgeot::lower_tri_solve(LU, x, N, true);
bgeot::upper_tri_solve(LU, x, N, false);
}
- scalar_type lu_det(const scalar_type *LU, const std::vector<int> &ipvt,
+ scalar_type lu_det(const scalar_type *LU, const std::vector<long> &ipvt,
size_type N) {
scalar_type det(1);
for (size_type j = 0; j < N; ++j) det *= *(LU+j*(N+1));
- for(int i = 0; i < int(N); ++i) if (i != ipvt[i]-1) { det = -det; }
+ for(long i = 0; i < long(N); ++i) if (i != ipvt[i]-1) { det = -det; }
return det;
}
@@ -209,7 +209,7 @@ namespace bgeot {
}
}
- void lu_inverse(const scalar_type *LU, const std::vector<int> &ipvt,
+ void lu_inverse(const scalar_type *LU, const std::vector<long> &ipvt,
scalar_type *A, size_type N) {
__aux2().resize(N); gmm::clear(__aux2());
__aux3().resize(N);
diff --git a/src/bgeot_sparse_tensors.cc b/src/bgeot_sparse_tensors.cc
index 0ac040b..fdf6ca7 100644
--- a/src/bgeot_sparse_tensors.cc
+++ b/src/bgeot_sparse_tensors.cc
@@ -907,10 +907,10 @@ namespace bgeot {
}
static bool do_reduction2v(bgeot::multi_tensor_iterator &mti) {
- int n = mti.vectorized_size();
+ long n = mti.vectorized_size();
const std::vector<stride_type> &s = mti.vectorized_strides();
if (n && s[0] && s[1] && s[2] == 0) {
- int incx = s[1], incy = s[0];
+ long incx = s[1], incy = s[0];
/*mti.print();
scalar_type *b[3];
for (int i=0; i < 3; ++i) b[i] = &mti.p(i);*/
@@ -945,10 +945,10 @@ namespace bgeot {
}
static bool do_reduction3v(bgeot::multi_tensor_iterator &mti) {
- int n = mti.vectorized_size();
+ long n = mti.vectorized_size();
const std::vector<stride_type> &s = mti.vectorized_strides();
if (n && s[0] && s[1] && s[2] == 0 && s[3] == 0) {
- int incx = s[1], incy = s[0];
+ long incx = s[1], incy = s[0];
do {
double v = mti.p(2)*mti.p(3);
gmm::daxpy_(&n, &v, &mti.p(1), &incx, &mti.p(0), &incy);
diff --git a/src/getfem/bgeot_geometric_trans.h
b/src/getfem/bgeot_geometric_trans.h
index c024a49..a8c2998 100644
--- a/src/getfem/bgeot_geometric_trans.h
+++ b/src/getfem/bgeot_geometric_trans.h
@@ -415,7 +415,7 @@ namespace bgeot {
mutable scalar_type J_, J__; /** Jacobian */
mutable base_matrix PC, B_factors;
mutable base_vector aux1, aux2;
- mutable std::vector<int> ipvt;
+ mutable std::vector<long> ipvt;
mutable bool have_J_, have_B_, have_B3_, have_B32_, have_K_,
have_cv_center_;
void compute_J() const;
public:
diff --git a/src/getfem_assembling_tensors.cc b/src/getfem_assembling_tensors.cc
index a58704d..af5cd84 100644
--- a/src/getfem_assembling_tensors.cc
+++ b/src/getfem_assembling_tensors.cc
@@ -22,9 +22,6 @@
#include "getfem/getfem_assembling_tensors.h"
#include "getfem/getfem_mat_elem.h"
-extern "C" void daxpy_(const int *n, const double *alpha, const double *x,
- const int *incx, double *y, const int *incy);
-
namespace getfem {
size_type vdim_specif_list::nb_mf() const {
return std::count_if(begin(),end(),
@@ -355,9 +352,9 @@ namespace getfem {
tensor_bases[k] = const_cast<TDIter>(&(*eltm[k]->begin()));
}
red.do_reduction();
- int one = 1, n = int(red.out_data.size()); assert(n);
- daxpy_(&n, &c, const_cast<double*>(&(red.out_data[0])),
- &one, (double*)&(t[0]), &one);
+ long one = 1, n = int(red.out_data.size()); assert(n);
+ gmm::daxpy_(&n, &c, const_cast<double*>(&(red.out_data[0])),
+ &one, (double*)&(t[0]), &one);
}
void resize_t(bgeot::base_tensor &t) {
bgeot::multi_index r;
diff --git a/src/getfem_contact_and_friction_large_sliding.cc
b/src/getfem_contact_and_friction_large_sliding.cc
index b1bb1e9..c52bfd1 100644
--- a/src/getfem_contact_and_friction_large_sliding.cc
+++ b/src/getfem_contact_and_friction_large_sliding.cc
@@ -1628,7 +1628,7 @@ namespace getfem {
gmm::add(gmm::identity_matrix(), grad);
gmm::mult(grad, ctx_y0.K(), gradtot);
- std::vector<int> ipvt(N);
+ std::vector<long> ipvt(N);
size_type info = gmm::lu_factor(gradtot, ipvt);
GMM_ASSERT1(!info, "Singular system, pivot = " << info); // il
faudrait faire qlq chose d'autre ... perturber par exemple
gmm::lu_solve(gradtot, ipvt, h, val);
diff --git a/src/getfem_generic_assembly.cc b/src/getfem_generic_assembly.cc
index 5843a7c..cec2231 100644
--- a/src/getfem_generic_assembly.cc
+++ b/src/getfem_generic_assembly.cc
@@ -4734,8 +4734,8 @@ namespace getfem {
virtual int exec() {
GA_DEBUG_INFO("Instruction: reduction operation of size " << nn);
#if GA_USES_BLAS
- int m = int(tc1.size()/nn), k = int(nn), n = int(tc2.size()/nn);
- int lda = m, ldb = n, ldc = m;
+ long m = int(tc1.size()/nn), k = int(nn), n = int(tc2.size()/nn);
+ long lda = m, ldb = n, ldc = m;
char T = 'T', N = 'N';
scalar_type alpha(1), beta(0);
gmm::dgemm_(&N, &T, &m, &n, &k, &alpha, &(tc1[0]), &lda, &(tc2[0]), &ldb,
@@ -13860,7 +13860,11 @@ namespace getfem {
loc_M.base_resize(nd, nd); gmm::resize(loc_U, nd);
gmm::sub_index J(mf.ind_basic_dof_of_element(v.cv()));
gmm::copy(gmm::sub_matrix(M, J, J), loc_M);
+ cout << "will loch lu_solve" << endl;
+ cout << "loc_M = " << loc_M << endl;
+ cout << "sub_F = " << gmm::sub_vector(F, J) << endl;
gmm::lu_solve(loc_M, loc_U, gmm::sub_vector(F, J));
+ cout << "result : " << vref(loc_U) << endl;
gmm::copy(loc_U, gmm::sub_vector(result, J));
}
}
diff --git a/src/getfem_mat_elem.cc b/src/getfem_mat_elem.cc
index a5498f6..ae7f3f7 100644
--- a/src/getfem_mat_elem.cc
+++ b/src/getfem_mat_elem.cc
@@ -26,9 +26,6 @@
#include "getfem/getfem_mat_elem.h"
#include "getfem/getfem_omp.h"
-extern "C" void daxpy_(const int *n, const double *alpha, const double *x,
- const int *incx, double *y, const int *incy);
-
namespace getfem {
/* ********************************************************************* */
/* Elementary matrices computation. */
@@ -384,19 +381,19 @@ namespace getfem {
if (nm == 0) {
t[0] += J;
} else {
- int n0 = int(es_end[0] - es_beg[0]);
+ long n0 = int(es_end[0] - es_beg[0]);
base_tensor::const_iterator pts0 = pts[0];
/* very heavy reduction .. takes much time */
k = nm-1; Vtab[k] = J;
- int one = 1;
+ long one = 1;
scalar_type V;
do {
for (V = Vtab[k]; k; --k)
Vtab[k-1] = V = *pts[k] * V;
GMM_ASSERT1(pt+n0 <= t.end(), "Internal error");
- daxpy_(&n0, &V, const_cast<double*>(&(pts0[0])), &one,
- (double*)&(*pt), &one);
+ gmm::daxpy_(&n0, &V, const_cast<double*>(&(pts0[0])), &one,
+ (double*)&(*pt), &one);
pt+=n0;
for (k=1; k != nm && ++pts[k] == es_end[k]; ++k)
pts[k] = es_beg[k];
diff --git a/src/getfem_mesher.cc b/src/getfem_mesher.cc
index 0127b2c..c9064e9 100644
--- a/src/getfem_mesher.cc
+++ b/src/getfem_mesher.cc
@@ -109,7 +109,7 @@ namespace getfem {
// cout << "nbco = " << nbco << endl;
std::vector<const mesher_signed_distance*> ls(nbco);
std::vector<scalar_type> d(nbco), v(nbco);
- std::vector<int> ipvt(nbco);
+ std::vector<long> ipvt(nbco);
gmm::col_matrix<base_node> G(N, nbco);
gmm::dense_matrix<scalar_type> H(nbco, nbco);
base_small_vector dd(N);
diff --git a/src/gmm/gmm_blas_interface.h b/src/gmm/gmm_blas_interface.h
index c41ae95..8144293 100644
--- a/src/gmm/gmm_blas_interface.h
+++ b/src/gmm/gmm_blas_interface.h
@@ -151,13 +151,13 @@ namespace gmm {
/* BLAS functions used. */
/* ********************************************************************* */
extern "C" {
- void daxpy_(const int *n, const double *alpha, const double *x,
- const int *incx, double *y, const int *incy);
- void dgemm_(const char *tA, const char *tB, const int *m,
- const int *n, const int *k, const double *alpha,
- const double *A, const int *ldA, const double *B,
- const int *ldB, const double *beta, double *C,
- const int *ldC);
+ void daxpy_(const long *n, const double *alpha, const double *x,
+ const long *incx, double *y, const long *incy);
+ void dgemm_(const char *tA, const char *tB, const long *m,
+ const long *n, const long *k, const double *alpha,
+ const double *A, const long *ldA, const double *B,
+ const long *ldB, const double *beta, double *C,
+ const long *ldC);
void sgemm_(...); void cgemm_(...); void zgemm_(...);
void sgemv_(...); void dgemv_(...); void cgemv_(...); void zgemv_(...);
void strsv_(...); void dtrsv_(...); void ctrsv_(...); void ztrsv_(...);
@@ -180,7 +180,7 @@ namespace gmm {
inline number_traits<base_type >::magnitude_type \
vect_norm2(param1(base_type)) { \
GMMLAPACK_TRACE("nrm2_interface"); \
- int inc(1), n(int(vect_size(x))); trans1(base_type); \
+ long inc(1), n(long(vect_size(x))); trans1(base_type); \
return blas_name(&n, &x[0], &inc); \
}
@@ -200,7 +200,7 @@ namespace gmm {
blas_name, base_type) \
inline base_type vect_sp(param1(base_type), param2(base_type)) { \
GMMLAPACK_TRACE("dot_interface"); \
- trans1(base_type); trans2(base_type); int inc(1), n(int(vect_size(y)));\
+ trans1(base_type); trans2(base_type); long inc(1), n(long(vect_size(y)));\
return mult1 mult2 blas_name(&n, &x[0], &inc, &y[0], &inc); \
}
@@ -267,7 +267,7 @@ namespace gmm {
blas_name, base_type) \
inline base_type vect_hp(param1(base_type), param2(base_type)) { \
GMMLAPACK_TRACE("dotc_interface"); \
- trans1(base_type); trans2(base_type); int inc(1), n(int(vect_size(y)));\
+ trans1(base_type); trans2(base_type); long inc(1), n(long(vect_size(y)));\
return mult1 mult2 blas_name(&n, &x[0], &inc, &y[0], &inc); \
}
@@ -332,7 +332,7 @@ namespace gmm {
# define axpy_interface(param1, trans1, blas_name, base_type) \
inline void add(param1(base_type), std::vector<base_type > &y) { \
GMMLAPACK_TRACE("axpy_interface"); \
- int inc(1), n(int(vect_size(y))); trans1(base_type); \
+ long inc(1), n(long(vect_size(y))); trans1(base_type); \
if (n == 0) return;
\
blas_name(&n, &a, &x[0], &inc, &y[0], &inc); \
}
@@ -367,7 +367,7 @@ namespace gmm {
std::vector<base_type > &z, orien) { \
GMMLAPACK_TRACE("gemv_interface"); \
trans1(base_type); trans2(base_type); base_type beta(1); \
- int m(int(mat_nrows(A))), lda(m), n(int(mat_ncols(A))), inc(1); \
+ long m(long(mat_nrows(A))), lda(m), n(long(mat_ncols(A))), inc(1); \
if (m && n) blas_name(&t, &m, &n, &alpha, &A(0,0), &lda, &x[0], &inc, \
&beta, &z[0], &inc); \
else gmm::clear(z); \
@@ -489,7 +489,7 @@ namespace gmm {
std::vector<base_type > &z, orien) { \
GMMLAPACK_TRACE("gemv_interface2"); \
trans1(base_type); trans2(base_type); base_type beta(0); \
- int m(int(mat_nrows(A))), lda(m), n(int(mat_ncols(A))), inc(1); \
+ long m(long(mat_nrows(A))), lda(m), n(long(mat_ncols(A))), inc(1); \
if (m && n) \
blas_name(&t, &m, &n, &alpha, &A(0,0), &lda, &x[0], &inc, &beta, \
&z[0], &inc); \
@@ -586,8 +586,8 @@ namespace gmm {
const std::vector<base_type > &V, \
const std::vector<base_type > &W) { \
GMMLAPACK_TRACE("ger_interface"); \
- int m(int(mat_nrows(A))), lda = m, n(int(mat_ncols(A))); \
- int incx = 1, incy = 1; \
+ long m(long(mat_nrows(A))), lda = m, n(long(mat_ncols(A)));
\
+ long incx = 1, incy = 1; \
base_type alpha(1); \
if (m && n)
\
blas_name(&m, &n, &alpha, &V[0], &incx, &W[0], &incy, &A(0,0), &lda);\
@@ -604,8 +604,8 @@ namespace gmm {
const std::vector<base_type > &W) { \
GMMLAPACK_TRACE("ger_interface"); \
gemv_trans2_s(base_type); \
- int m(int(mat_nrows(A))), lda = m, n(int(mat_ncols(A))); \
- int incx = 1, incy = 1; \
+ long m(long(mat_nrows(A))), lda = m, n(long(mat_ncols(A)));
\
+ long incx = 1, incy = 1; \
if (m && n)
\
blas_name(&m, &n, &alpha, &x[0], &incx, &W[0], &incy, &A(0,0), &lda);\
}
@@ -621,8 +621,8 @@ namespace gmm {
gemv_p2_s(base_type)) { \
GMMLAPACK_TRACE("ger_interface"); \
gemv_trans2_s(base_type); \
- int m(int(mat_nrows(A))), lda = m, n(int(mat_ncols(A))); \
- int incx = 1, incy = 1; \
+ long m(long(mat_nrows(A))), lda = m, n(long(mat_ncols(A)));
\
+ long incx = 1, incy = 1; \
base_type al2 = gmm::conj(alpha); \
if (m && n)
\
blas_name(&m, &n, &al2, &V[0], &incx, &x[0], &incy, &A(0,0), &lda); \
@@ -643,9 +643,9 @@ namespace gmm {
dense_matrix<base_type > &C, c_mult) { \
GMMLAPACK_TRACE("gemm_interface_nn"); \
const char t = 'N'; \
- int m(int(mat_nrows(A))), lda = m, k(int(mat_ncols(A))); \
- int n(int(mat_ncols(B))); \
- int ldb = k, ldc = m; \
+ long m(long(mat_nrows(A))), lda = m, k(long(mat_ncols(A)));
\
+ long n(long(mat_ncols(B)));
\
+ long ldb = k, ldc = m; \
base_type alpha(1), beta(0); \
if (m && k && n) \
blas_name(&t, &t, &m, &n, &k, &alpha, \
@@ -671,8 +671,8 @@ namespace gmm {
dense_matrix<base_type > &A \
= const_cast<dense_matrix<base_type > &>(*(linalg_origin(A_))); \
const char t = 'T', u = 'N'; \
- int m(int(mat_ncols(A))), k(int(mat_nrows(A))), n(int(mat_ncols(B))); \
- int lda = k, ldb = k, ldc = m; \
+ long m(long(mat_ncols(A))), k(long(mat_nrows(A))), n(long(mat_ncols(B))); \
+ long lda = k, ldb = k, ldc = m; \
base_type alpha(1), beta(0); \
if (m && k && n) \
blas_name(&t, &u, &m, &n, &k, &alpha, \
@@ -701,9 +701,9 @@ namespace gmm {
dense_matrix<base_type > &B \
= const_cast<dense_matrix<base_type > &>(*(linalg_origin(B_))); \
const char t = 'N', u = 'T'; \
- int m(int(mat_nrows(A))), lda = m, k(int(mat_ncols(A))); \
- int n(int(mat_nrows(B))); \
- int ldb = n, ldc = m; \
+ long m(long(mat_nrows(A))), lda = m, k(long(mat_ncols(A))); \
+ long n(long(mat_nrows(B)));
\
+ long ldb = n, ldc = m; \
base_type alpha(1), beta(0); \
if (m && k && n) \
blas_name(&t, &u, &m, &n, &k, &alpha, \
@@ -735,8 +735,8 @@ namespace gmm {
dense_matrix<base_type > &B \
= const_cast<dense_matrix<base_type > &>(*(linalg_origin(B_))); \
const char t = 'T', u = 'T'; \
- int m(int(mat_ncols(A))), k(int(mat_nrows(A))), n(int(mat_nrows(B))); \
- int lda = k, ldb = n, ldc = m; \
+ long m(long(mat_ncols(A))), k(long(mat_nrows(A))), n(long(mat_nrows(B))); \
+ long lda = k, ldb = n, ldc = m; \
base_type alpha(1), beta(0); \
if (m && k && n) \
blas_name(&t, &u, &m, &n, &k, &alpha, \
@@ -775,8 +775,8 @@ namespace gmm {
dense_matrix<base_type > &A \
= const_cast<dense_matrix<base_type > &>(*(linalg_origin(A_))); \
const char t = 'C', u = 'N'; \
- int m(int(mat_ncols(A))), k(int(mat_nrows(A))), n(int(mat_ncols(B))); \
- int lda = k, ldb = k, ldc = m; \
+ long m(long(mat_ncols(A))), k(long(mat_nrows(A))), n(long(mat_ncols(B))); \
+ long lda = k, ldb = k, ldc = m; \
base_type alpha(1), beta(0); \
if (m && k && n) \
blas_name(&t, &u, &m, &n, &k, &alpha, \
@@ -801,8 +801,8 @@ namespace gmm {
dense_matrix<base_type > &B \
= const_cast<dense_matrix<base_type > &>(*(linalg_origin(B_))); \
const char t = 'N', u = 'C'; \
- int m(int(mat_nrows(A))), lda = m, k(int(mat_ncols(A))); \
- int n(int(mat_nrows(B))), ldb = n, ldc = m;
\
+ long m(long(mat_nrows(A))), lda = m, k(long(mat_ncols(A)));
\
+ long n(long(mat_nrows(B))), ldb = n, ldc = m; \
base_type alpha(1), beta(0); \
if (m && k && n) \
blas_name(&t, &u, &m, &n, &k, &alpha, \
@@ -830,8 +830,8 @@ namespace gmm {
dense_matrix<base_type > &B \
= const_cast<dense_matrix<base_type > &>(*(linalg_origin(B_))); \
const char t = 'C', u = 'C'; \
- int m(int(mat_ncols(A))), k(int(mat_nrows(A))), lda = k; \
- int n(int(mat_nrows(B))), ldb = n, ldc = m;
\
+ long m(long(mat_ncols(A))), k(long(mat_nrows(A))), lda = k;
\
+ long n(long(mat_nrows(B))), ldb = n, ldc = m; \
base_type alpha(1), beta(0); \
if (m && k && n) \
blas_name(&t, &u, &m, &n, &k, &alpha, \
@@ -853,7 +853,7 @@ namespace gmm {
size_type k, bool is_unit) { \
GMMLAPACK_TRACE("trsv_interface"); \
loru; trans1(base_type); char d = is_unit ? 'U' : 'N'; \
- int lda(int(mat_nrows(A))), inc(1), n = int(k); \
+ long lda(long(mat_nrows(A))), inc(1), n = long(k); \
if (lda) blas_name(&l, &t, &d, &n, &A(0,0), &lda, &x[0], &inc); \
}
diff --git a/src/gmm/gmm_dense_lu.h b/src/gmm/gmm_dense_lu.h
index 5107abe..cb53771 100644
--- a/src/gmm/gmm_dense_lu.h
+++ b/src/gmm/gmm_dense_lu.h
@@ -138,7 +138,7 @@ namespace gmm {
void lu_solve(const DenseMatrix &A, VectorX &x, const VectorB &b) {
typedef typename linalg_traits<DenseMatrix>::value_type T;
dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
- std::vector<int> ipvt(mat_nrows(A));
+ std::vector<long> ipvt(mat_nrows(A));
gmm::copy(A, B);
size_type info = lu_factor(B, ipvt);
GMM_ASSERT1(!info, "Singular system, pivot = " << info);
@@ -212,7 +212,7 @@ namespace gmm {
typedef typename linalg_traits<DenseMatrix>::value_type T;
DenseMatrix& A = const_cast<DenseMatrix&>(A_);
dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
- std::vector<int> ipvt(mat_nrows(A));
+ std::vector<long> ipvt(mat_nrows(A));
gmm::copy(A, B);
size_type info = lu_factor(B, ipvt);
if (doassert) GMM_ASSERT1(!info, "Non invertible matrix, pivot = "<<info);
@@ -238,7 +238,7 @@ namespace gmm {
lu_det(const DenseMatrix& A) {
typedef typename linalg_traits<DenseMatrix>::value_type T;
dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
- std::vector<int> ipvt(mat_nrows(A));
+ std::vector<long> ipvt(mat_nrows(A));
gmm::copy(A, B);
lu_factor(B, ipvt);
return lu_det(B, ipvt);
diff --git a/src/gmm/gmm_lapack_interface.h b/src/gmm/gmm_lapack_interface.h
index 7888aea..3995569 100644
--- a/src/gmm/gmm_lapack_interface.h
+++ b/src/gmm/gmm_lapack_interface.h
@@ -47,42 +47,42 @@
namespace gmm {
- /* ********************************************************************* */
- /* Operations interfaced for T = float, double, std::complex<float> */
- /* or std::complex<double> : */
- /* */
- /* lu_factor(dense_matrix<T>, std::vector<int>) */
- /* lu_solve(dense_matrix<T>, std::vector<T>, std::vector<T>) */
- /* lu_solve(dense_matrix<T>, std::vector<int>, std::vector<T>, */
- /* std::vector<T>) */
- /* lu_solve_transposed(dense_matrix<T>, std::vector<int>, std::vector<T>,*/
- /* std::vector<T>) */
- /* lu_inverse(dense_matrix<T>) */
- /* lu_inverse(dense_matrix<T>, std::vector<int>, dense_matrix<T>) */
- /* */
- /* qr_factor(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>) */
- /* */
- /* implicit_qr_algorithm(dense_matrix<T>, std::vector<T>) */
- /* implicit_qr_algorithm(dense_matrix<T>, std::vector<T>, */
- /* dense_matrix<T>) */
- /* implicit_qr_algorithm(dense_matrix<T>, std::vector<std::complex<T> >) */
- /* implicit_qr_algorithm(dense_matrix<T>, std::vector<std::complex<T> >, */
- /* dense_matrix<T>) */
- /* */
- /* geev_interface_right */
- /* geev_interface_left */
- /* */
- /* schur(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>) */
- /* */
- /* svd(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>, std::vector<T>)*/
- /* svd(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>, */
- /* std::vector<std::complex<T> >) */
- /* */
- /* ********************************************************************* */
-
- /* ********************************************************************* */
- /* LAPACK functions used. */
- /* ********************************************************************* */
+ /* ********************************************************************** */
+ /* Operations interfaced for T = float, double, std::complex<float> */
+ /* or std::complex<double> : */
+ /* */
+ /* lu_factor(dense_matrix<T>, std::vector<long>) */
+ /* lu_solve(dense_matrix<T>, std::vector<T>, std::vector<T>) */
+ /* lu_solve(dense_matrix<T>, std::vector<long>, std::vector<T>, */
+ /* std::vector<T>) */
+ /* lu_solve_transposed(dense_matrix<T>, std::vector<long>, std::vector<T>,*/
+ /* std::vector<T>) */
+ /* lu_inverse(dense_matrix<T>) */
+ /* lu_inverse(dense_matrix<T>, std::vector<long>, dense_matrix<T>) */
+ /* */
+ /* qr_factor(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>) */
+ /* */
+ /* implicit_qr_algorithm(dense_matrix<T>, std::vector<T>) */
+ /* implicit_qr_algorithm(dense_matrix<T>, std::vector<T>, */
+ /* dense_matrix<T>) */
+ /* implicit_qr_algorithm(dense_matrix<T>, std::vector<std::complex<T> >) */
+ /* implicit_qr_algorithm(dense_matrix<T>, std::vector<std::complex<T> >, */
+ /* dense_matrix<T>) */
+ /* */
+ /* geev_interface_right */
+ /* geev_interface_left */
+ /* */
+ /* schur(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>) */
+ /* */
+ /* svd(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>, std::vector<T>) */
+ /* svd(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>, */
+ /* std::vector<std::complex<T> >) */
+ /* */
+ /* ********************************************************************** */
+
+ /* ********************************************************************** */
+ /* LAPACK functions used. */
+ /* ********************************************************************** */
extern "C" {
void sgetrf_(...); void dgetrf_(...); void cgetrf_(...); void zgetrf_(...);
@@ -97,16 +97,33 @@ namespace gmm {
void sgesvd_(...); void dgesvd_(...); void cgesvd_(...); void zgesvd_(...);
}
- /* ********************************************************************* */
- /* LU decomposition. */
- /* ********************************************************************* */
-
-# define getrf_interface(lapack_name, base_type) inline \
- size_type lu_factor(dense_matrix<base_type > &A, std::vector<int> &ipvt){\
- GMMLAPACK_TRACE("getrf_interface"); \
- int m = int(mat_nrows(A)), n = int(mat_ncols(A)), lda(m), info(0); \
- if (m && n) lapack_name(&m, &n, &A(0,0), &lda, &ipvt[0], &info); \
- return size_type(info); \
+ // For compatibility with lapack version with 32 bit integer.
+# define int_to_long_ipvt(ipvt) \
+ { \
+ for (size_type ii = ipvt.size(); ii != 0; --ii) \
+ ipvt[ii-1] = long(((int *)(&ipvt[0]))[ii-1]); \
+ }
+
+// # define long_to_int_ipvt(ipvt)
+// {
+// long *ipvt_ = const_cast<long *>(&ipvt[0]);
+// for (size_type ii = 0; ii < ipvt.size(); ++ii)
+// ((int *)(&ipvt_[0]))[ii] = int(ipvt_[ii]);
+// }
+
+ /* ********************************************************************** */
+ /* LU decomposition. */
+ /* ********************************************************************** */
+
+# define getrf_interface(lapack_name, base_type) inline \
+ size_type lu_factor(dense_matrix<base_type > &A, std::vector<long> &ipvt){ \
+ GMMLAPACK_TRACE("getrf_interface"); \
+ long m = long(mat_nrows(A)), n = long(mat_ncols(A)), lda(m), info(-1L); \
+ if (m && n) lapack_name(&m, &n, &A(0,0), &lda, &ipvt[0], &info); \
+ if ((info & 0xFFFFFFFF00000000L) && !(info & 0x00000000FFFFFFFFL)) \
+ /* For compatibility with lapack version with 32 bit integer. */ \
+ int_to_long_ipvt(ipvt); \
+ return size_type(int(info & 0x00000000FFFFFFFFL)); \
}
getrf_interface(sgetrf_, BLAS_S)
@@ -114,47 +131,24 @@ namespace gmm {
getrf_interface(cgetrf_, BLAS_C)
getrf_interface(zgetrf_, BLAS_Z)
- /* ********************************************************************* */
- /* LU solve. */
- /* ********************************************************************* */
-
-# define getrs_interface(f_name, trans1, lapack_name, base_type) inline \
- void f_name(const dense_matrix<base_type > &A, \
- const std::vector<int> &ipvt, std::vector<base_type > &x, \
- const std::vector<base_type > &b) { \
- GMMLAPACK_TRACE("getrs_interface"); \
- int n = int(mat_nrows(A)), info, nrhs(1); \
- gmm::copy(b, x); trans1; \
- if (n) \
- lapack_name(&t, &n, &nrhs, &(A(0,0)),&n,&ipvt[0], &x[0], &n, &info); \
- }
+ /* ********************************************************************** */
+ /* LU inverse. */
+ /* ********************************************************************** */
-# define getrs_trans_n const char t = 'N'
-# define getrs_trans_t const char t = 'T'
-
- getrs_interface(lu_solve, getrs_trans_n, sgetrs_, BLAS_S)
- getrs_interface(lu_solve, getrs_trans_n, dgetrs_, BLAS_D)
- getrs_interface(lu_solve, getrs_trans_n, cgetrs_, BLAS_C)
- getrs_interface(lu_solve, getrs_trans_n, zgetrs_, BLAS_Z)
- getrs_interface(lu_solve_transposed, getrs_trans_t, sgetrs_, BLAS_S)
- getrs_interface(lu_solve_transposed, getrs_trans_t, dgetrs_, BLAS_D)
- getrs_interface(lu_solve_transposed, getrs_trans_t, cgetrs_, BLAS_C)
- getrs_interface(lu_solve_transposed, getrs_trans_t, zgetrs_, BLAS_Z)
-
- /* ********************************************************************* */
- /* LU inverse. */
- /* ********************************************************************* */
-
# define getri_interface(lapack_name, base_type) inline \
void lu_inverse(const dense_matrix<base_type > &LU, \
- std::vector<int> &ipvt, const dense_matrix<base_type > &A_) { \
+ std::vector<long> &ipvt, const dense_matrix<base_type > &A_) { \
GMMLAPACK_TRACE("getri_interface"); \
dense_matrix<base_type> &A \
= const_cast<dense_matrix<base_type > &>(A_); \
- int n = int(mat_nrows(A)), info, lwork(10000); base_type work[10000]; \
+ long n=long(mat_nrows(A)), info(-1L), lwork(10000);
\
+ base_type work[10000]; \
if (n) { \
std::copy(LU.begin(), LU.end(), A.begin()); \
lapack_name(&n, &A(0,0), &n, &ipvt[0], &work[0], &lwork, &info); \
+ if ((info & 0xFFFFFFFF00000000L) && !(info & 0x00000000FFFFFFFFL)) \
+ /* For compatibility with lapack version with 32 bit integer. */ \
+ int_to_long_ipvt(ipvt);
\
} \
}
@@ -164,19 +158,19 @@ namespace gmm {
getri_interface(zgetri_, BLAS_Z)
- /* ********************************************************************* */
- /* QR factorization. */
- /* ********************************************************************* */
+ /* ********************************************************************** */
+ /* QR factorization. */
+ /* ********************************************************************** */
# define geqrf_interface(lapack_name1, base_type) inline \
void qr_factor(dense_matrix<base_type > &A){ \
GMMLAPACK_TRACE("geqrf_interface"); \
- int m = int(mat_nrows(A)), n = int(mat_ncols(A)), info, lwork(-1); \
+ long m = long(mat_nrows(A)), n=long(mat_ncols(A)), info(0), lwork(-1); \
base_type work1; \
if (m && n) { \
std::vector<base_type > tau(n); \
lapack_name1(&m, &n, &A(0,0), &m, &tau[0], &work1 , &lwork, &info); \
- lwork = int(gmm::real(work1)); \
+ lwork = long(gmm::real(work1)); \
std::vector<base_type > work(lwork); \
lapack_name1(&m, &n, &A(0,0), &m, &tau[0], &work[0], &lwork, &info); \
GMM_ASSERT1(!info, "QR factorization failed"); \
@@ -194,19 +188,19 @@ namespace gmm {
void qr_factor(const dense_matrix<base_type > &A, \
dense_matrix<base_type > &Q, dense_matrix<base_type > &R) { \
GMMLAPACK_TRACE("geqrf_interface2"); \
- int m = int(mat_nrows(A)), n = int(mat_ncols(A)), info, lwork(-1); \
+ long m = long(mat_nrows(A)), n=long(mat_ncols(A)), info(0), lwork(-1); \
base_type work1; \
if (m && n) { \
std::copy(A.begin(), A.end(), Q.begin());
\
std::vector<base_type > tau(n); \
lapack_name1(&m, &n, &Q(0,0), &m, &tau[0], &work1 , &lwork, &info); \
- lwork = int(gmm::real(work1)); \
+ lwork = long(gmm::real(work1)); \
std::vector<base_type > work(lwork); \
lapack_name1(&m, &n, &Q(0,0), &m, &tau[0], &work[0], &lwork, &info); \
GMM_ASSERT1(!info, "QR factorization failed"); \
base_type *p = &R(0,0), *q = &Q(0,0); \
- for (int j = 0; j < n; ++j, q += m-n) \
- for (int i = 0; i < n; ++i, ++p, ++q) \
+ for (long j = 0; j < n; ++j, q += m-n) \
+ for (long i = 0; i < n; ++i, ++p, ++q) \
*p = (j < i) ? base_type(0) : *q; \
lapack_name2(&m, &n, &n, &Q(0,0), &m,&tau[0],&work[0],&lwork,&info); \
} \
@@ -218,9 +212,9 @@ namespace gmm {
geqrf_interface2(cgeqrf_, cungqr_, BLAS_C)
geqrf_interface2(zgeqrf_, zungqr_, BLAS_Z)
- /* ********************************************************************* */
- /* QR algorithm for eigenvalues search. */
- /* ********************************************************************* */
+ /* ********************************************************************** */
+ /* QR algorithm for eigenvalues search. */
+ /* ********************************************************************** */
# define gees_interface(lapack_name, base_type) \
template <typename VECT> inline void implicit_qr_algorithm( \
@@ -229,14 +223,14 @@ namespace gmm {
double tol=gmm::default_tol(base_type()), bool compvect = true) { \
GMMLAPACK_TRACE("gees_interface"); \
typedef bool (*L_fp)(...); L_fp p = 0; \
- int n = int(mat_nrows(A)), info, lwork(-1), sdim; base_type work1; \
+ long n=long(mat_nrows(A)), info(0), lwork(-1), sdim; base_type work1; \
if (!n) return; \
dense_matrix<base_type > H(n,n); gmm::copy(A, H); \
char jobvs = (compvect ? 'V' : 'N'), sort = 'N'; \
std::vector<double> rwork(n), eigv1(n), eigv2(n); \
lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigv1[0], \
&eigv2[0], &Q(0,0), &n, &work1, &lwork, &rwork[0], &info); \
- lwork = int(gmm::real(work1)); \
+ lwork = long(gmm::real(work1)); \
std::vector<base_type > work(lwork); \
lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigv1[0], \
&eigv2[0], &Q(0,0), &n, &work[0], &lwork, &rwork[0],&info);\
@@ -251,14 +245,14 @@ namespace gmm {
double tol=gmm::default_tol(base_type()), bool compvect = true) { \
GMMLAPACK_TRACE("gees_interface2"); \
typedef bool (*L_fp)(...); L_fp p = 0; \
- int n = int(mat_nrows(A)), info, lwork(-1), sdim; base_type work1; \
+ long n=long(mat_nrows(A)), info(0), lwork(-1), sdim; base_type work1; \
if (!n) return; \
dense_matrix<base_type > H(n,n); gmm::copy(A, H); \
char jobvs = (compvect ? 'V' : 'N'), sort = 'N'; \
std::vector<double> rwork(n), eigvv(n*2); \
lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigvv[0], \
&Q(0,0), &n, &work1, &lwork, &rwork[0], &rwork[0], &info); \
- lwork = int(gmm::real(work1)); \
+ lwork = long(gmm::real(work1)); \
std::vector<base_type > work(lwork); \
lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigvv[0], \
&Q(0,0), &n, &work[0], &lwork, &rwork[0], &rwork[0],&info);\
@@ -276,18 +270,18 @@ namespace gmm {
# define jobv_left char jobvl = 'V', jobvr = 'N';
# define geev_interface(lapack_name, base_type, side) \
- template <typename VECT> inline void geev_interface_ ## side( \
+ template <typename VECT> inline void geev_interface_ ## side( \
const dense_matrix<base_type > &A, const VECT &eigval_, \
dense_matrix<base_type > &Q) { \
GMMLAPACK_TRACE("geev_interface"); \
- int n = int(mat_nrows(A)), info, lwork(-1); base_type work1; \
+ long n = long(mat_nrows(A)), info(0), lwork(-1); base_type work1; \
if (!n) return; \
dense_matrix<base_type > H(n,n); gmm::copy(A, H); \
jobv_ ## side \
std::vector<base_type > eigvr(n), eigvi(n); \
lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigvr[0], &eigvi[0], \
&Q(0,0), &n, &Q(0,0), &n, &work1, &lwork, &info); \
- lwork = int(gmm::real(work1)); \
+ lwork = long(gmm::real(work1)); \
std::vector<base_type > work(lwork); \
lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigvr[0], &eigvi[0], \
&Q(0,0), &n, &Q(0,0), &n, &work[0], &lwork, &info); \
@@ -301,7 +295,7 @@ namespace gmm {
const dense_matrix<base_type > &A, const VECT &eigval_, \
dense_matrix<base_type > &Q) { \
GMMLAPACK_TRACE("geev_interface"); \
- int n = int(mat_nrows(A)), info, lwork(-1); base_type work1; \
+ long n = long(mat_nrows(A)), info(0), lwork(-1); base_type work1; \
if (!n) return; \
dense_matrix<base_type > H(n,n); gmm::copy(A, H); \
jobv_ ## side \
@@ -309,7 +303,7 @@ namespace gmm {
std::vector<base_type> eigv(n); \
lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigv[0], &Q(0,0), &n, \
&Q(0,0), &n, &work1, &lwork, &rwork[0], &info); \
- lwork = int(gmm::real(work1)); \
+ lwork = long(gmm::real(work1)); \
std::vector<base_type > work(lwork); \
lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigv[0], &Q(0,0), &n, \
&Q(0,0), &n, &work[0], &lwork, &rwork[0], &info); \
@@ -328,32 +322,32 @@ namespace gmm {
geev_interface2(zgeev_, BLAS_Z, left)
- /* ********************************************************************* */
- /* SCHUR algorithm: */
- /* A = Q*S*(Q^T), with Q orthogonal and S upper quasi-triangula */
- /* ********************************************************************* */
+ /* ********************************************************************** */
+ /* SCHUR algorithm: */
+ /* A = Q*S*(Q^T), with Q orthogonal and S upper quasi-triangula */
+ /* ********************************************************************** */
# define geesx_interface(lapack_name, base_type) inline \
void schur(dense_matrix<base_type> &A, \
dense_matrix<base_type> &S, \
dense_matrix<base_type> &Q) { \
GMMLAPACK_TRACE("geesx_interface"); \
- int m = int(mat_nrows(A)), n = int(mat_ncols(A)); \
+ long m = long(mat_nrows(A)), n = long(mat_ncols(A)); \
GMM_ASSERT1(m == n, "Schur decomposition requires square matrix"); \
char jobvs = 'V', sort = 'N', sense = 'N'; \
bool select = false; \
- int lwork = 8*n, sdim = 0, liwork = 1; \
+ long lwork = 8*n, sdim = 0, liwork = 1; \
std::vector<base_type> work(lwork), wr(n), wi(n); \
- std::vector<int> iwork(liwork); \
- std::vector<int> bwork(1); \
+ std::vector<long> iwork(liwork); \
+ std::vector<long> bwork(1); \
resize(S, n, n); copy(A, S); \
resize(Q, n, n); \
base_type rconde(0), rcondv(0); \
- int info = -1; \
+ long info(0); \
lapack_name(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \
&sdim, &wr[0], &wi[0], &Q(0,0), &n, &rconde, &rcondv, \
&work[0], &lwork, &iwork[0], &liwork, &bwork[0], &info);\
- GMM_ASSERT1(!info, "SCHUR algorithm failed"); \
+ GMM_ASSERT1(!info, "SCHUR algorithm failed"); \
}
# define geesx_interface2(lapack_name, base_type) inline \
@@ -361,18 +355,18 @@ namespace gmm {
dense_matrix<base_type> &S, \
dense_matrix<base_type> &Q) { \
GMMLAPACK_TRACE("geesx_interface"); \
- int m = int(mat_nrows(A)), n = int(mat_ncols(A)); \
+ long m = long(mat_nrows(A)), n = long(mat_ncols(A)); \
GMM_ASSERT1(m == n, "Schur decomposition requires square matrix"); \
char jobvs = 'V', sort = 'N', sense = 'N'; \
bool select = false; \
- int lwork = 8*n, sdim = 0; \
+ long lwork = 8*n, sdim = 0; \
std::vector<base_type::value_type> rwork(lwork); \
std::vector<base_type> work(lwork), w(n); \
- std::vector<int> bwork(1); \
+ std::vector<long> bwork(1); \
resize(S, n, n); copy(A, S); \
resize(Q, n, n); \
base_type rconde(0), rcondv(0); \
- int info = -1; \
+ long info(0); \
lapack_name(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \
&sdim, &w[0], &Q(0,0), &n, &rconde, &rcondv, \
&work[0], &lwork, &rwork[0], &bwork[0], &info); \
@@ -391,10 +385,10 @@ namespace gmm {
}
- /* ********************************************************************* */
- /* Interface to SVD. Does not correspond to a Gmm++ functionnality. */
- /* Author : Sebastian Nowozin <address@hidden> */
- /* ********************************************************************* */
+ /* ********************************************************************** */
+ /* Interface to SVD. Does not correspond to a Gmm++ functionnality. */
+ /* Author : Sebastian Nowozin <address@hidden> */
+ /* ********************************************************************** */
# define gesvd_interface(lapack_name, base_type) inline \
void svd(dense_matrix<base_type> &X, \
@@ -402,15 +396,15 @@ namespace gmm {
dense_matrix<base_type> &Vtransposed, \
std::vector<base_type> &sigma) { \
GMMLAPACK_TRACE("gesvd_interface"); \
- int m = int(mat_nrows(X)), n = int(mat_ncols(X)); \
- int mn_min = m < n ? m : n; \
+ long m = long(mat_nrows(X)), n = long(mat_ncols(X)); \
+ long mn_min = m < n ? m : n; \
sigma.resize(mn_min); \
std::vector<base_type> work(15 * mn_min); \
- int lwork = int(work.size()); \
+ long lwork = long(work.size()); \
resize(U, m, m); \
resize(Vtransposed, n, n); \
char job = 'A'; \
- int info = -1; \
+ long info(0); \
lapack_name(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \
&m, &Vtransposed(0,0), &n, &work[0], &lwork, &info); \
}
@@ -421,16 +415,16 @@ namespace gmm {
dense_matrix<base_type> &Vtransposed, \
std::vector<base_type2> &sigma) { \
GMMLAPACK_TRACE("gesvd_interface"); \
- int m = int(mat_nrows(X)), n = int(mat_ncols(X)); \
- int mn_min = m < n ? m : n; \
+ long m = long(mat_nrows(X)), n = long(mat_ncols(X)); \
+ long mn_min = m < n ? m : n; \
sigma.resize(mn_min); \
std::vector<base_type> work(15 * mn_min); \
std::vector<base_type2> rwork(5 * mn_min); \
- int lwork = int(work.size()); \
+ long lwork = long(work.size()); \
resize(U, m, m); \
resize(Vtransposed, n, n); \
char job = 'A'; \
- int info = -1; \
+ long info(0); \
lapack_name(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \
&m, &Vtransposed(0,0), &n, &work[0], &lwork, \
&rwork[0], &info); \
diff --git a/src/gmm/gmm_opt.h b/src/gmm/gmm_opt.h
index 4dbe05f..31e830e 100644
--- a/src/gmm/gmm_opt.h
+++ b/src/gmm/gmm_opt.h
@@ -110,7 +110,7 @@ namespace gmm {
}
default : {
dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
- std::vector<int> ipvt(mat_nrows(A));
+ std::vector<long> ipvt(mat_nrows(A));
gmm::copy(A, B);
size_type info = lu_factor(B, ipvt);
GMM_ASSERT1(!info, "non invertible matrix");
diff --git a/tests/test_gmm_matrix_functions.cc
b/tests/test_gmm_matrix_functions.cc
index 9589926..b13c27f 100644
--- a/tests/test_gmm_matrix_functions.cc
+++ b/tests/test_gmm_matrix_functions.cc
@@ -58,8 +58,8 @@ template <class T> void test_logm(T) {
gmm::copy(gmm::identity_matrix(), X);
X(0,1) = -T(1);
Z(0,1) = -T(1);
- gmm::logm(X, LOGMX);
std::cout << "X = " << X << std::endl;
+ gmm::logm(X, LOGMX);
std::cout << "logm(X) = " << LOGMX << std::endl;
gmm::add(gmm::scaled(Z, -T(1)), LOGMX);
std::cout << "err(logm(X)) = " << LOGMX << std::endl;
@@ -94,19 +94,15 @@ int main(void)
{
srand(1459);
- try {
-
- test_sqrtm(float());
- test_sqrtm(double());
- test_sqrtm(std::complex<float>());
- test_sqrtm(std::complex<double>());
-
- test_logm(float());
- test_logm(double());
- test_logm(std::complex<float>());
- test_logm(std::complex<double>());
- }
- GMM_STANDARD_CATCH_ERROR;
-
+ test_sqrtm(float());
+ test_sqrtm(double());
+ test_sqrtm(std::complex<float>());
+ test_sqrtm(std::complex<double>());
+
+ test_logm(float());
+ test_logm(double());
+ test_logm(std::complex<float>());
+ test_logm(std::complex<double>());
+
return 0;
}