[Getfem-commits] [getfem-commits] branch master updated: Remove hack for treating blas64 in gmm

[Konstantinos Poulios]

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 <logari81@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);     \
     }                                                                      \
   }