[Getfem-commits] [getfem-commits] branch master updated: Allow building and running GetFEM without external BLAS

[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 70981fd9 Allow building and running GetFEM without external BLAS
70981fd9 is described below

commit 70981fd95eb3445dfba81015e56235ebce58e613
Author: Konstantinos Poulios <logari81@gmail.com>
AuthorDate: Tue Dec 17 09:44:13 2024 +0100

    Allow building and running GetFEM without external BLAS
---
 configure.ac                     | 18 +++++++---------
 src/bgeot_sparse_tensors.cc      | 45 +++++++++++++++++++++++++++++++++-------
 src/getfem_assembling_tensors.cc | 14 ++++++++++---
 src/getfem_mat_elem.cc           | 13 ++++++++++--
 src/gmm/gmm_lapack_interface.h   | 29 +++-----------------------
 5 files changed, 71 insertions(+), 48 deletions(-)

diff --git a/configure.ac b/configure.ac
index 887242f3..81afda21 100644
--- a/configure.ac
+++ b/configure.ac
@@ -672,15 +672,13 @@ LIBS="$acx_blas_save_LIBS"
 
 # Finally, execute ACTION-IF-FOUND/ACTION-IF-NOT-FOUND:
 if test x"$acx_blas_ok" = xyes; then
-  echo "OK, You have working BLAS libs ! Using $BLAS_LIBS" ; HAVE_VENDOR_BLAS=1
+  echo "OK, You have working BLAS libs ! Using $BLAS_LIBS"
+  AC_DEFINE(GMM_USES_BLAS,,[defined if GMM is linked to a blas library])
+  LIBS="$LIBS $BLAS_LIBS"
 else
-  echo " *** YOU DONT HAVE BLAS! *** Using a cheap replacement" ; 
HAVE_VENDOR_BLAS=0
+  echo " *** YOU DONT HAVE BLAS! *** Using a cheap replacement"
 fi
 
-dnl ACX_BLAS([ echo "OK, You have working BLAS libs !"; HAVE_VENDOR_BLAS=1 
],[echo "YOU DONT HAVE BLAS! Using a cheap replacement" ; HAVE_VENDOR_BLAS=0])
-LIBS="$LIBS $BLAS_LIBS"
-AC_DEFINE(GMM_USES_BLAS,,[defined if GMM is linked to a blas library])
-
 useblas64support=NO
 AC_ARG_ENABLE(blas64-support,
         [AS_HELP_STRING([--enable-blas64-support],[enable the 64 bits integer 
blas and lapack support])],
@@ -1327,15 +1325,15 @@ fi;
 if test x"$acx_lapack_ok" = xyes; then
   echo "- Lapack library found: $LAPACK_LIBS"
 else
-  echo "- Lapack library not found: generic (less effective) algorithms will 
be used"
+  echo "- Lapack library not found: generic (less efficient) algorithms will 
be used"
 fi
 
-if test "x$HAVE_VENDOR_BLAS" = "x0"; then
+if test x"$acx_blas_ok" = xyes; then
+  echo "- BLAS library found. Link options: $BLAS_LIBS"
+else
   echo "- *** No usable blas library was found ***"
   echo "  A generic BLAS implementation will be used, however you should "
   echo "  consider installing a faster BLAS, such as ATLAS"
-else
-  echo "- BLAS library found. Link options: $BLAS_LIBS"
 fi;
 echo "  You can give the location of your prefered blas library with either"
 echo "  the --with-blas=<lib> option, or the BLAS_LIBS environment variable"
diff --git a/src/bgeot_sparse_tensors.cc b/src/bgeot_sparse_tensors.cc
index 7d08de4a..4ccbf9a8 100644
--- a/src/bgeot_sparse_tensors.cc
+++ b/src/bgeot_sparse_tensors.cc
@@ -927,17 +927,33 @@ namespace bgeot {
   }
 
   static bool do_reduction2v(bgeot::multi_tensor_iterator &mti) {
-    BLAS_INT n = mti.vectorized_size();
+    size_type n = mti.vectorized_size();
     const std::vector<stride_type> &s = mti.vectorized_strides();
     if (n && s[0] && s[1] && s[2] == 0) {
-      BLAS_INT incx = s[1], incy = s[0];
+#if defined(GMM_USES_BLAS)
+      BLAS_INT nn = n, incx = s[1], incy = s[0];
+#else
+      dim_type incx = dim_type(s[1]), incy = dim_type(s[0]);
+#endif
       /*mti.print();
         scalar_type *b[3];
         for (int i=0; i < 3; ++i)       b[i] = &mti.p(i);*/
       do {
         /*cout << "vectorized_ reduction2a : n=" << n << ", s = " << s << " 
mti.p=" << &mti.p(0)-b[0] << ","
           << &mti.p(1)-b[1] << "," << &mti.p(2)-b[2] << "\n";*/
-        gmm::daxpy_(&n, &mti.p(2), &mti.p(1), &incx, &mti.p(0), &incy);
+#if defined(GMM_USES_BLAS)
+        gmm::daxpy_(&nn, &mti.p(2), &mti.p(1), &incx, &mti.p(0), &incy);
+#else
+        double a = mti.p(2);
+        scalar_type* itx = &mti.p(1);
+        scalar_type* ity = &mti.p(0);
+        *ity += a * (*itx);
+        for (size_type i=1; i < n; ++i) {
+          itx += incx;
+          ity += incy;
+          *ity += a * (*itx);
+        }
+#endif
       } while (mti.vnext());
       return true;
     } else
@@ -966,13 +982,28 @@ namespace bgeot {
   }
 
   static bool do_reduction3v(bgeot::multi_tensor_iterator &mti) {
-    BLAS_INT n = mti.vectorized_size();
+    size_type 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) {
-      BLAS_INT incx = s[1], incy = s[0];
+#if defined(GMM_USES_BLAS)
+      BLAS_INT nn = n, incx = s[1], incy = s[0];
+#else
+      dim_type incx = dim_type(s[1]), incy = dim_type(s[0]);
+#endif
       do {
-        double v = mti.p(2)*mti.p(3);
-        gmm::daxpy_(&n, &v, &mti.p(1), &incx, &mti.p(0), &incy);
+        double a = mti.p(2)*mti.p(3);
+#if defined(GMM_USES_BLAS)
+        gmm::daxpy_(&nn, &a, &mti.p(1), &incx, &mti.p(0), &incy);
+#else
+        scalar_type* itx = &mti.p(1);
+        scalar_type* ity = &mti.p(0);
+        *ity += a * (*itx);
+        for (size_type i=1; i < n; ++i) {
+          itx += incx;
+          ity += incy;
+          *ity += a * (*itx);
+        }
+#endif
       } while (mti.vnext());
       return true;
     } else
diff --git a/src/getfem_assembling_tensors.cc b/src/getfem_assembling_tensors.cc
index 91963dae..bd5bb389 100644
--- a/src/getfem_assembling_tensors.cc
+++ b/src/getfem_assembling_tensors.cc
@@ -331,7 +331,7 @@ namespace getfem {
   /* called (if possible, i.e. if not an exact integration) for each
   integration point during mat_elem->compute() */
   struct computed_tensor_integration_callback
-    : public mat_elem_integration_callback {
+  : public mat_elem_integration_callback {
       bgeot::tensor_reduction red;
       bool was_called;
       std::vector<TDIter> tensor_bases; /* each tref of 'red' has a   */
@@ -347,10 +347,18 @@ namespace getfem {
           tensor_bases[k] = const_cast<TDIter>(&(*eltm[k]->begin()));
         }
         red.do_reduction();
+#if defined(GMM_USES_BLAS)
         BLAS_INT one = BLAS_INT(1), n = BLAS_INT(red.out_data.size());
         assert(n);
-       gmm::daxpy_(&n, &c, const_cast<double*>(&(red.out_data[0])),
-                   &one, (double*)&(t[0]), &one);
+        gmm::daxpy_(&n, &c, const_cast<double*>(&(red.out_data[0])),
+                    &one, (double*)&(t[0]), &one);
+#else
+        size_type n = red.out_data.size();
+        assert(n);
+        for (size_type k=0; k < n; ++k)
+          t[k] += c * red.out_data[k];
+        // gmm::add(gmm::scaled(red.out_data, c), t.as_vector())
+#endif
       }
       void resize_t(bgeot::base_tensor &t) {
         bgeot::multi_index r;
diff --git a/src/getfem_mat_elem.cc b/src/getfem_mat_elem.cc
index 66557426..0cd9ff83 100644
--- a/src/getfem_mat_elem.cc
+++ b/src/getfem_mat_elem.cc
@@ -400,20 +400,29 @@ namespace getfem {
       if (nm == 0)
         t[0] += J;
       else {
+#if defined(GMM_USES_BLAS)
         BLAS_INT n0 = BLAS_INT(es_end[0] - es_beg[0]);
+        BLAS_INT one = BLAS_INT(1);
+#else
+        size_type n0 = size_type(es_end[0] - es_beg[0]);
+#endif
         base_tensor::const_iterator pts0 = pts[0];
 
         /* very heavy reduction .. takes much time */
         k = nm-1; Vtab[k] = J;
-        BLAS_INT one = BLAS_INT(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");
+#if defined(GMM_USES_BLAS)
           gmm::daxpy_(&n0, &V, const_cast<double*>(&(pts0[0])), &one,
                       (double*)&(*pt), &one);
-          pt+=n0;
+          pt += n0;
+#else
+          for (k=0; k < n0; ++k)
+            *pt++ += V*pts0[k];
+#endif
           for (k=1; k != nm && ++pts[k] == es_end[k]; ++k)
             pts[k] = es_beg[k];
         } while (k != nm);
diff --git a/src/gmm/gmm_lapack_interface.h b/src/gmm/gmm_lapack_interface.h
index 135721fb..ff427b1e 100644
--- a/src/gmm/gmm_lapack_interface.h
+++ b/src/gmm/gmm_lapack_interface.h
@@ -467,37 +467,14 @@ void schur(const MAT &, MAT &, MAT &)
                      "to be built with Lapack");
 }
 
-inline void svd(dense_matrix<BLAS_S> &, dense_matrix<BLAS_S> &,
-         dense_matrix<BLAS_S> &, std::vector<BLAS_S> &)
+template <typename BLAS_TYPE>
+inline void svd(dense_matrix<BLAS_TYPE> &, dense_matrix<BLAS_TYPE> &,
+         dense_matrix<BLAS_TYPE> &, std::vector<BLAS_TYPE> &)
 {
   GMM_ASSERT1(false, "Use of function svd(X,U,Vtransposed,sigma) requires 
GetFEM "
                      "to be built with Lapack");
 }
 
-inline void svd(dense_matrix<BLAS_D> &, dense_matrix<BLAS_D> &,
-         dense_matrix<BLAS_D> &, std::vector<BLAS_D> &)
-{
-  GMM_ASSERT1(false, "Use of function svd(X,U,Vtransposed,sigma) requires 
GetFEM "
-              "to be built with Lapack");
-}
-
-inline void svd(dense_matrix<BLAS_C> &, dense_matrix<BLAS_C> &,
-         dense_matrix<BLAS_C> &, std::vector<BLAS_S> &)
-{
-  GMM_ASSERT1(false, "Use of function svd(X,U,Vtransposed,sigma) requires 
GetFEM "
-              "to be built with Lapack");
-}
-
-inline void svd(dense_matrix<BLAS_Z> &, dense_matrix<BLAS_Z> &,
-         dense_matrix<BLAS_Z> &, std::vector<BLAS_D> &)
-{
-  GMM_ASSERT1(false, "Use of function svd(X,U,Vtransposed,sigma) requires 
GetFEM "
-              "to be built with Lapack");
-}
-
-
-  
-
 }// namespace gmm
 
 #endif // GMM_USES_LAPACK