From 83bf100bb2fa40a817c5ffef53e3a62e631caf7f Mon Sep 17 00:00:00 2001
From: conrad <conrad@localhost.localdomain>
Date: Sun, 2 Apr 2023 23:40:51 +1000
Subject: [PATCH] speedup via inplace processing

---
 include/armadillo_bits/SpSubview_meat.hpp | 63 ++++++++++++++++++++++-
 1 file changed, 62 insertions(+), 1 deletion(-)

diff --git a/include/armadillo_bits/SpSubview_meat.hpp b/include/armadillo_bits/SpSubview_meat.hpp
index 4e350ff9..481359be 100644
--- a/include/armadillo_bits/SpSubview_meat.hpp
+++ b/include/armadillo_bits/SpSubview_meat.hpp
@@ -420,7 +420,68 @@ SpSubview<eT>::operator/=(const Base<eT, T1>& x)
   {
   arma_extra_debug_sigprint();
   
-  return (*this).operator=( (*this) / x.get_ref() );
+  const SpSubview<eT>& A = (*this);
+  
+  const quasi_unwrap<T1> U(x.get_ref());
+  const Mat<eT>& B     = U.M;
+  
+  arma_debug_assert_same_size(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "element-wise division");
+  
+  bool result_ok = true;
+  
+  constexpr eT zero = eT(0);
+  
+  const uword B_n_rows = B.n_rows;
+  const uword B_n_cols = B.n_cols;
+  
+  for(uword c=0; c < B_n_cols; ++c)
+    {
+    for(uword r=0; r < B_n_rows; ++r)
+      {
+      // a zero in B and A at the same location implies the division result is NaN;
+      // hence a zero in A (not stored) needs to be changed into a non-zero
+      
+      // for efficiency, an element in B is checked before checking the corresponding element in A
+      
+      if((B.at(r,c) == zero) && (A.at(r,c) == zero))  { result_ok = false; break; }
+      }
+    
+    if(result_ok == false)  { break; }
+    }
+  
+  if(result_ok)
+    {
+    const_iterator cit     = A.begin();
+    const_iterator cit_end = A.end();
+    
+    while(cit != cit_end)
+      {
+      const eT tmp = (*cit) / B.at(cit.row(), cit.col());
+      
+      if(tmp == zero)  { result_ok = false; break; }
+      
+      ++cit;
+      }
+    }
+  
+  if(result_ok)
+    {
+    iterator it     = (*this).begin();
+    iterator it_end = (*this).end();
+    
+    while(it != it_end)
+      {
+      (*it) /= B.at(it.row(), it.col());
+      
+      ++it;
+      }
+    }
+  else
+    {
+    (*this).operator=( (*this) / B );
+    }
+  
+  return (*this);
   }