This is an automated email from the git hooks/post-receive script. logari81 pushed a commit to branch master in repository getfem.
The following commit(s) were added to refs/heads/master by this push: new cbe1442c Remove hack for treating blas64 in gmm cbe1442c is described below commit cbe1442c4bb2b4fea804d47054ab55dc0bfe238c Author: Konstantinos Poulios <logar...@gmail.com> AuthorDate: Tue Mar 26 13:00:00 2024 +0100 Remove hack for treating blas64 in gmm - BLAS integer type must be fixed at configure time in gmm_arch_config.h through the GMM_USE_BLAS64_INTERFACE macro - lu_factor interface returns 0 if successful, or a pivot number if a pivot is negative or 0, it does not distinguish between the latter two --- src/gmm/gmm_dense_lu.h | 68 +++++++++++------------------------------- src/gmm/gmm_lapack_interface.h | 16 ++++------ 2 files changed, 24 insertions(+), 60 deletions(-) diff --git a/src/gmm/gmm_dense_lu.h b/src/gmm/gmm_dense_lu.h index 52102366..4823f3f3 100644 --- a/src/gmm/gmm_dense_lu.h +++ b/src/gmm/gmm_dense_lu.h @@ -73,47 +73,11 @@ namespace gmm { - /* ********************************************************************** */ - /* IPVT structure. */ - /* ********************************************************************** */ - // For compatibility with lapack version with 64 or 32 bit integer. - // Should be replaced by std::vector<size_type> if 32 bit integer version - // of lapack is not used anymore (and lapack_ipvt_int set to size_type) - - // Do not use iterators of this interface container - class lapack_ipvt : public std::vector<size_type> { - bool is_int64; - size_type &operator[](size_type i) - { return std::vector<size_type>::operator[](i); } - size_type operator[] (size_type i) const - { return std::vector<size_type>::operator[](i); } - void begin() const {} - void begin() {} - void end() const {} - void end() {} - - public: - void set_to_int32() { is_int64 = false; } - const size_type *pfirst() const - { return &(*(std::vector<size_type>::begin())); } - size_type *pfirst() { return &(*(std::vector<size_type>::begin())); } - - lapack_ipvt(size_type n) : std::vector<size_type>(n), is_int64(true) {} - - size_type get(size_type i) const { - const size_type *p = pfirst(); - return is_int64 ? p[i] : size_type(((const int *)(p))[i]); - } - void set(size_type i, size_type val) { - size_type *p = pfirst(); - if (is_int64) p[i] = val; else ((int *)(p))[i] = int(val); - } - }; -} - -#include "gmm_opt.h" - -namespace gmm { +#if defined(GMM_USES_BLAS) || defined(GMM_USES_LAPACK) + typedef std::vector<BLAS_INT> lapack_ipvt; +#else + typedef std::vector<size_type> lapack_ipvt; +#endif /** LU Factorization of a general (dense) matrix (real or complex). @@ -125,9 +89,10 @@ namespace gmm { The pivot indices in ipvt are indexed starting from 1 so that this is compatible with LAPACK (Fortran). */ - template <typename DenseMatrix> - size_type lu_factor(DenseMatrix& A, lapack_ipvt& ipvt) { + template <typename DenseMatrix, typename Pvector> + size_type lu_factor(DenseMatrix& A, Pvector& ipvt) { typedef typename linalg_traits<DenseMatrix>::value_type T; + typedef typename linalg_traits<Pvector>::value_type INT; typedef typename number_traits<T>::magnitude_type R; size_type info(0), i, j, jp, M(mat_nrows(A)), N(mat_ncols(A)); if (M == 0 || N == 0) @@ -136,14 +101,14 @@ namespace gmm { std::vector<T> c(M), r(N); GMM_ASSERT2(ipvt.size()+1 >= NN, "IPVT too small"); - for (i = 0; i+1 < NN; ++i) ipvt.set(i, i); + for (i = 0; i+1 < NN; ++i) ipvt[i] = INT(i); if (M || N) { for (j = 0; j+1 < NN; ++j) { R max = gmm::abs(A(j,j)); jp = j; for (i = j+1; i < M; ++i) /* find pivot. */ if (gmm::abs(A(i,j)) > max) { jp = i; max = gmm::abs(A(i,j)); } - ipvt.set(j, jp + 1); + ipvt[j] = INT(jp + 1); if (max == R(0)) { info = j + 1; break; } if (jp != j) for (i = 0; i < N; ++i) std::swap(A(jp, i), A(j, i)); @@ -153,7 +118,7 @@ namespace gmm { rank_one_update(sub_matrix(A, sub_interval(j+1, M-j-1), sub_interval(j+1, N-j-1)), c, conjugated(r)); } - ipvt.set(NN-1, NN); + ipvt[NN-1] = INT(NN); } return info; } @@ -167,7 +132,7 @@ namespace gmm { typedef typename linalg_traits<DenseMatrix>::value_type T; copy(b, x); for(size_type i = 0; i < pvector.size(); ++i) { - size_type perm = pvector.get(i)-1; // permutations stored in 1's offset + size_type perm = size_type(pvector[i]-1); // permutations stored in 1's offset if (i != perm) { T aux = x[i]; x[i] = x[perm]; x[perm] = aux; } } /* solve Ax = b -> LUx = b -> Ux = L^-1 b. */ @@ -198,7 +163,7 @@ namespace gmm { lower_tri_solve(transposed(LU), x, false); upper_tri_solve(transposed(LU), x, true); for (size_type i = pvector.size(); i > 0; --i) { - size_type perm = pvector.get(i-1)-1; // permutations stored in 1's offset + size_type perm = size_type(pvector[i-1]-1); // permutations stored in 1's offset if (i-1 != perm) { T aux = x[i-1]; x[i-1] = x[perm]; @@ -275,12 +240,13 @@ namespace gmm { typename linalg_traits<DenseMatrixLU>::value_type lu_det(const DenseMatrixLU& LU, const Pvector &pvector) { typedef typename linalg_traits<DenseMatrixLU>::value_type T; + typedef typename linalg_traits<Pvector>::value_type INT; T det(1); const size_type J=std::min(mat_nrows(LU), mat_ncols(LU)); for (size_type j = 0; j < J; ++j) det *= LU(j,j); - for(size_type i = 0; i < pvector.size(); ++i) - if (i != size_type(pvector.get(i)-1)) { det = -det; } + for(INT i = 0; i < INT(pvector.size()); ++i) + if (i != pvector[i]-1) { det = -det; } return det; } @@ -300,5 +266,7 @@ namespace gmm { } +#include "gmm_opt.h" + #endif diff --git a/src/gmm/gmm_lapack_interface.h b/src/gmm/gmm_lapack_interface.h index 9a71da8d..04f1802f 100644 --- a/src/gmm/gmm_lapack_interface.h +++ b/src/gmm/gmm_lapack_interface.h @@ -105,13 +105,9 @@ namespace gmm { size_type lu_factor(dense_matrix<base_type> &A, lapack_ipvt &ipvt){ \ GMMLAPACK_TRACE("getrf_interface"); \ BLAS_INT m = BLAS_INT(mat_nrows(A)), n = BLAS_INT(mat_ncols(A)), lda(m); \ - BLAS_INT info(-1); \ - if (m && n) lapack_name(&m, &n, &A(0,0), &lda, ipvt.pfirst(), &info); \ - if ((sizeof(BLAS_INT) == 4) || \ - ((info & 0xFFFFFFFF00000000L) && !(info & 0x00000000FFFFFFFFL))) \ - /* For compatibility with lapack version with 32 bit integer. */ \ - ipvt.set_to_int32(); \ - return size_type(int(info & 0x00000000FFFFFFFFL)); \ + BLAS_INT info=BLAS_INT(-1); \ + if (m && n) lapack_name(&m, &n, &A(0,0), &lda, &ipvt[0], &info); \ + return size_type(abs(info)); \ } getrf_interface(sgetrf_, BLAS_S) @@ -131,7 +127,7 @@ namespace gmm { BLAS_INT n = BLAS_INT(mat_nrows(A)), info(0), nrhs(1); \ gmm::copy(b, x); trans1; \ if (n) \ - lapack_name(&t,&n,&nrhs,&(A(0,0)),&n,ipvt.pfirst(),&x[0],&n,&info); \ + lapack_name(&t,&n,&nrhs,&(A(0,0)),&n,&ipvt[0],&x[0],&n,&info); \ } # define getrs_trans_n const char t = 'N' @@ -158,10 +154,10 @@ namespace gmm { base_type work1; \ if (n) { \ gmm::copy(LU, A); \ - lapack_name(&n, &A(0,0), &n, ipvt.pfirst(), &work1, &lwork, &info); \ + lapack_name(&n, &A(0,0), &n, &ipvt[0], &work1, &lwork, &info); \ lwork = int(gmm::real(work1)); \ std::vector<base_type> work(lwork); \ - lapack_name(&n, &A(0,0), &n, ipvt.pfirst(), &work[0], &lwork,&info); \ + lapack_name(&n, &A(0,0), &n, &ipvt[0], &work[0], &lwork, &info); \ } \ }