From 3ee77a6ded8b2590a59fa8dd0bebefd645fe308c Mon Sep 17 00:00:00 2001
From: conrad <conrad@localhost.localdomain>
Date: Thu, 18 Nov 2021 14:30:04 +1000
Subject: [PATCH] optimised handling of diagonal matrices

---
 include/armadillo_bits/op_inv_meat.hpp | 39 +++++++++++++++-----------
 1 file changed, 22 insertions(+), 17 deletions(-)

diff --git a/include/armadillo_bits/op_inv_meat.hpp b/include/armadillo_bits/op_inv_meat.hpp
index 69a46958..e8af684c 100644
--- a/include/armadillo_bits/op_inv_meat.hpp
+++ b/include/armadillo_bits/op_inv_meat.hpp
@@ -306,7 +306,7 @@ op_inv_sympd::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename
   arma_extra_debug_sigprint();
   
   typedef typename T1::elem_type eT;
-  // typedef typename T1::pod_type   T;
+  typedef typename T1::pod_type   T;
   
   out = expr.get_ref();
   
@@ -329,22 +329,27 @@ op_inv_sympd::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename
     // fallthrough if optimisation failed
     }
   
-  // if((is_cx<eT>::no) && (out.is_diagmat()))
-  //   {
-  //   const uword N = (std::min)(out.n_rows, out.n_cols);
-  //   
-  //   for(uword i=0; i<N; ++i)
-  //     {
-  //           eT&      out_ii = out.at(i,i);
-  //     const  T  real_out_ii = access::tmp_real(out_ii);
-  //     
-  //     if(real_out_ii <= T(0))  { return false; }
-  //     
-  //     out_ii = eT(T(1) / real_out_ii);
-  //     }
-  //     
-  //   return true;
-  //   }
+  if((is_cx<eT>::no) && (out.is_diagmat()))
+    {
+    // specialised handling of real matrices only;
+    // currently auxlib::inv_sympd() does not enforce that 
+    // imaginary components of diagonal elements must be zero;
+    // strictly enforcing this constraint may break existing user software.
+    
+    const uword N = (std::min)(out.n_rows, out.n_cols);
+    
+    for(uword i=0; i<N; ++i)
+      {
+            eT&      out_ii = out.at(i,i);
+      const  T  real_out_ii = access::tmp_real(out_ii);
+      
+      if(real_out_ii <= T(0))  { return false; }
+      
+      out_ii = eT(T(1) / real_out_ii);
+      }
+      
+    return true;
+    }
   
   return auxlib::inv_sympd(out);
   }