From 4865cc317b25eac292cf92ecae2ff7b9c8b2c1e2 Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Thu, 14 Sep 2023 11:55:21 +0200
Subject: [PATCH] added: allow specifying offset and stride for vector in add

---
 src/LinAlg/matrix.h | 20 ++++++++++++--------
 1 file changed, 12 insertions(+), 8 deletions(-)

diff --git a/src/LinAlg/matrix.h b/src/LinAlg/matrix.h
index 2e7ba033..71989da3 100644
--- a/src/LinAlg/matrix.h
+++ b/src/LinAlg/matrix.h
@@ -128,7 +128,8 @@ namespace utl //! General utility classes and functions.
     //! \brief Subtract the given vector \b X from \a *this.
     vector<T>& operator-=(const vector<T>& X) { return this->add(X,T(-1)); }
     //! \brief Add the given vector \b X scaled by \a alfa to \a *this.
-    vector<T>& add(const std::vector<T>& X, const T& alfa = T(1));
+    vector<T>& add(const std::vector<T>& X, const T& alfa = T(1),
+                   int ofsx = 0, int stridex = 1);
 
     //! \brief Perform \f${\bf Y} = \alpha{\bf Y} + (1-\alpha){\bf X} \f$
     //! where \b Y = \a *this.
@@ -963,19 +964,21 @@ namespace utl //! General utility classes and functions.
 
   template<> inline
   vector<float>& vector<float>::add(const std::vector<float>& X,
-                                    const float& alfa)
+                                    const float& alfa,
+                                    int ofsx, int stridex)
   {
     size_t n = this->size() < X.size() ? this->size() : X.size();
-    cblas_saxpy(n,alfa,&X.front(),1,this->ptr(),1);
+    cblas_saxpy(n,alfa,X.data()+ofsx,stridex,this->data(),1);
     return *this;
   }
 
   template<> inline
   vector<double>& vector<double>::add(const std::vector<double>& X,
-                                      const double& alfa)
+                                      const double& alfa,
+                                      int ofsx, int stridex)
   {
     size_t n = this->size() < X.size() ? this->size() : X.size();
-    cblas_daxpy(n,alfa,&X.front(),1,this->ptr(),1);
+    cblas_daxpy(n,alfa,X.data()+ofsx,stridex,this->data(),1);
     return *this;
   }
 
@@ -1345,11 +1348,12 @@ namespace utl //! General utility classes and functions.
   }
 
   template<class T> inline
-  vector<T>& vector<T>::add(const std::vector<T>& X, const T& alfa)
+  vector<T>& vector<T>::add(const std::vector<T>& X, const T& alfa,
+                            int ofsx, int stridex)
   {
     T* p = this->ptr();
-    const T* q = &X.front();
-    for (size_t i = 0; i < this->size() && i < X.size(); i++, p++, q++)
+    const T* q = X.data() + ofsx;
+    for (size_t i = 0; i < this->size() && i < X.size(); i++, p++, q += stridex)
       *p += alfa*(*q);
     return *this;
   }