From cc01e135c1c1f9335dc3c6fa81eacae795c19d92 Mon Sep 17 00:00:00 2001
From: conrad <conrad@localhost.localdomain>
Date: Tue, 22 Mar 2022 16:33:41 +1000
Subject: [PATCH] refactor handling of tiny matrices

---
 include/armadillo_bits/op_inv_gen_bones.hpp |   8 +-
 include/armadillo_bits/op_inv_gen_meat.hpp  | 276 +++++++++++---------
 include/armadillo_bits/op_inv_spd_bones.hpp |  10 +-
 include/armadillo_bits/op_inv_spd_meat.hpp  | 100 +++----
 4 files changed, 225 insertions(+), 169 deletions(-)

diff --git a/include/armadillo_bits/op_inv_gen_bones.hpp b/include/armadillo_bits/op_inv_gen_bones.hpp
index 42d2345b..0dda8b3c 100644
--- a/include/armadillo_bits/op_inv_gen_bones.hpp
+++ b/include/armadillo_bits/op_inv_gen_bones.hpp
@@ -55,7 +55,13 @@ class op_inv_gen_full
   inline static bool apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type,T1>& expr, const char* caller_sig, const uword flags);
   
   template<typename eT>
-  arma_cold inline static bool apply_tiny(Mat<eT>& X);
+  arma_cold inline static bool apply_tiny_2x2(Mat<eT>& X);
+  
+  template<typename eT>
+  arma_cold inline static bool apply_tiny_3x3(Mat<eT>& X);
+  
+  template<typename eT>
+  arma_cold inline static bool apply_tiny_4x4(Mat<eT>& X);
   };
 
 
diff --git a/include/armadillo_bits/op_inv_gen_meat.hpp b/include/armadillo_bits/op_inv_gen_meat.hpp
index 76f3fb57..f21eac1a 100644
--- a/include/armadillo_bits/op_inv_gen_meat.hpp
+++ b/include/armadillo_bits/op_inv_gen_meat.hpp
@@ -103,16 +103,42 @@ op_inv_gen_full::apply_direct(Mat<typename T1::elem_type>& out, const Base<typen
   
   arma_debug_check( (out.is_square() == false), caller_sig, ": given matrix must be square sized" );
   
-  if(tiny && (out.n_rows <= 4) && is_cx<eT>::no)
+  const uword N = out.n_rows;
+  
+  if(N == 0)  { return true; }
+  
+  if(is_cx<eT>::no)
     {
-    arma_extra_debug_print("op_inv_gen_full: attempting tinymatrix optimisation");
-    
-    const bool status = op_inv_gen_full::apply_tiny(out);
-    
-    if(status)  { return true; }
-    
-    arma_extra_debug_print("op_inv_gen_full: tinymatrix optimisation failed");
-    
+    if(N == 1)
+      {
+      const eT a = out[0];
+      
+      out[0] = eT(1) / a;
+      
+      return (a != eT(0));
+      }
+    else
+    if(N == 2)
+      {
+      const bool status = op_inv_gen_full::apply_tiny_2x2(out);
+      
+      if(status)  { return true; }
+      }
+    else
+    if((N == 3) && tiny)
+      {
+      const bool status = op_inv_gen_full::apply_tiny_3x3(out);
+      
+      if(status)  { return true; }
+      }
+    else
+    if((N == 4) && tiny)
+      {
+      const bool status = op_inv_gen_full::apply_tiny_4x4(out);
+      
+      if(status)  { return true; }
+      }
+      
     // fallthrough if optimisation failed
     }
   
@@ -120,8 +146,6 @@ op_inv_gen_full::apply_direct(Mat<typename T1::elem_type>& out, const Base<typen
     {
     arma_extra_debug_print("op_inv_gen_full: detected diagonal matrix");
     
-    const uword N = out.n_rows;
-    
     eT* colmem = out.memptr();
     
     for(uword i=0; i<N; ++i)
@@ -184,7 +208,7 @@ template<typename eT>
 arma_cold
 inline
 bool
-op_inv_gen_full::apply_tiny(Mat<eT>& X)
+op_inv_gen_full::apply_tiny_2x2(Mat<eT>& X)
   {
   arma_extra_debug_sigprint();
   
@@ -192,117 +216,133 @@ op_inv_gen_full::apply_tiny(Mat<eT>& X)
   
   // NOTE: assuming matrix X is square sized
   
-  const uword N = X.n_rows;
-  
   constexpr T det_min =        std::numeric_limits<T>::epsilon();
   constexpr T det_max = T(1) / std::numeric_limits<T>::epsilon();
   
   eT* Xm = X.memptr();
   
-       if(N == 0)  { return true; }
-  else if(N == 1)
-    {
-    const eT a = Xm[0];
-    
-    if(a == eT(0))  { return false; }
-    
-    Xm[0] = eT(1) / a;
-    
-    return true;
-    }
-  else if(N == 2)
-    {
-    const eT a = Xm[pos<0,0>::n2];
-    const eT b = Xm[pos<0,1>::n2];
-    const eT c = Xm[pos<1,0>::n2];
-    const eT d = Xm[pos<1,1>::n2];
-    
-    const eT     det_val = (a*d - b*c);
-    const  T abs_det_val = std::abs(det_val);
-    
-    if((abs_det_val < det_min) || (abs_det_val > det_max))  { return false; }
-    
-    Xm[pos<0,0>::n2] =  d / det_val;
-    Xm[pos<0,1>::n2] = -b / det_val;
-    Xm[pos<1,0>::n2] = -c / det_val;
-    Xm[pos<1,1>::n2] =  a / det_val;
-    
-    return true;
-    }
-  else if(N == 3)
-    {
-    Mat<eT> out(3, 3, arma_nozeros_indicator());
-    
-    eT* outm = out.memptr();
-    
-    const eT     det_val = op_det::apply_tiny(X);
-    const  T abs_det_val = std::abs(det_val);
-    
-    if((abs_det_val < det_min) || (abs_det_val > det_max))  { return false; }
-    
-    outm[pos<0,0>::n3] =  (Xm[pos<2,2>::n3]*Xm[pos<1,1>::n3] - Xm[pos<2,1>::n3]*Xm[pos<1,2>::n3]) / det_val;
-    outm[pos<1,0>::n3] = -(Xm[pos<2,2>::n3]*Xm[pos<1,0>::n3] - Xm[pos<2,0>::n3]*Xm[pos<1,2>::n3]) / det_val;
-    outm[pos<2,0>::n3] =  (Xm[pos<2,1>::n3]*Xm[pos<1,0>::n3] - Xm[pos<2,0>::n3]*Xm[pos<1,1>::n3]) / det_val;
-    
-    outm[pos<0,1>::n3] = -(Xm[pos<2,2>::n3]*Xm[pos<0,1>::n3] - Xm[pos<2,1>::n3]*Xm[pos<0,2>::n3]) / det_val;
-    outm[pos<1,1>::n3] =  (Xm[pos<2,2>::n3]*Xm[pos<0,0>::n3] - Xm[pos<2,0>::n3]*Xm[pos<0,2>::n3]) / det_val;
-    outm[pos<2,1>::n3] = -(Xm[pos<2,1>::n3]*Xm[pos<0,0>::n3] - Xm[pos<2,0>::n3]*Xm[pos<0,1>::n3]) / det_val;
-    
-    outm[pos<0,2>::n3] =  (Xm[pos<1,2>::n3]*Xm[pos<0,1>::n3] - Xm[pos<1,1>::n3]*Xm[pos<0,2>::n3]) / det_val;
-    outm[pos<1,2>::n3] = -(Xm[pos<1,2>::n3]*Xm[pos<0,0>::n3] - Xm[pos<1,0>::n3]*Xm[pos<0,2>::n3]) / det_val;
-    outm[pos<2,2>::n3] =  (Xm[pos<1,1>::n3]*Xm[pos<0,0>::n3] - Xm[pos<1,0>::n3]*Xm[pos<0,1>::n3]) / det_val;
-    
-    const eT check_val = Xm[pos<0,0>::n3]*outm[pos<0,0>::n3] + Xm[pos<0,1>::n3]*outm[pos<1,0>::n3] + Xm[pos<0,2>::n3]*outm[pos<2,0>::n3];
-    
-    const  T max_diff  = (is_float<T>::value) ? T(1e-4) : T(1e-10);  // empirically determined; may need tuning
-    
-    if(std::abs(T(1) - check_val) >= max_diff)  { return false; }
-    
-    arrayops::copy(X.memptr(), out.memptr(), out.n_elem);
-    }
-  else if(N == 4)
-    {
-    Mat<eT> out(4, 4, arma_nozeros_indicator());
-    
-    eT* outm = out.memptr();
-    
-    const eT     det_val = op_det::apply_tiny(X);
-    const  T abs_det_val = std::abs(det_val);
-    
-    if((abs_det_val < det_min) || (abs_det_val > det_max))  { return false; }
-    
-    outm[pos<0,0>::n4] = ( Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] + Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
-    outm[pos<1,0>::n4] = ( Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
-    outm[pos<2,0>::n4] = ( Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] ) / det_val;
-    outm[pos<3,0>::n4] = ( Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] ) / det_val;
-    
-    outm[pos<0,1>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
-    outm[pos<1,1>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
-    outm[pos<2,1>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] ) / det_val;
-    outm[pos<3,1>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] ) / det_val;
-    
-    outm[pos<0,2>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
-    outm[pos<1,2>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
-    outm[pos<2,2>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,3>::n4] ) / det_val;
-    outm[pos<3,2>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,2>::n4] ) / det_val;
-    
-    outm[pos<0,3>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4] ) / det_val;
-    outm[pos<1,3>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4] ) / det_val;
-    outm[pos<2,3>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4] ) / det_val;
-    outm[pos<3,3>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4] ) / det_val;
-    
-    const eT check_val = Xm[pos<0,0>::n4]*outm[pos<0,0>::n4] + Xm[pos<0,1>::n4]*outm[pos<1,0>::n4] + Xm[pos<0,2>::n4]*outm[pos<2,0>::n4] + Xm[pos<0,3>::n4]*outm[pos<3,0>::n4];
-    
-    const  T max_diff  = (is_float<T>::value) ? T(1e-4) : T(1e-10);  // empirically determined; may need tuning
-    
-    if(std::abs(T(1) - check_val) >= max_diff)  { return false; }
-    
-    arrayops::copy(X.memptr(), out.memptr(), out.n_elem);
-    }
-  else
-    {
-    return false;
-    }
+  const eT a = Xm[pos<0,0>::n2];
+  const eT b = Xm[pos<0,1>::n2];
+  const eT c = Xm[pos<1,0>::n2];
+  const eT d = Xm[pos<1,1>::n2];
+  
+  const eT     det_val = (a*d - b*c);
+  const  T abs_det_val = std::abs(det_val);
+  
+  if((abs_det_val < det_min) || (abs_det_val > det_max))  { return false; }
+  
+  Xm[pos<0,0>::n2] =  d / det_val;
+  Xm[pos<0,1>::n2] = -b / det_val;
+  Xm[pos<1,0>::n2] = -c / det_val;
+  Xm[pos<1,1>::n2] =  a / det_val;
+  
+  return true;
+  }
+
+
+
+template<typename eT>
+arma_cold
+inline
+bool
+op_inv_gen_full::apply_tiny_3x3(Mat<eT>& X)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename get_pod_type<eT>::result T;
+  
+  // NOTE: assuming matrix X is square sized
+  
+  constexpr T det_min =        std::numeric_limits<T>::epsilon();
+  constexpr T det_max = T(1) / std::numeric_limits<T>::epsilon();
+  
+  Mat<eT> Y(3, 3, arma_nozeros_indicator());
+  
+  eT* Xm = X.memptr();
+  eT* Ym = Y.memptr();
+  
+  const eT     det_val = op_det::apply_tiny(X);
+  const  T abs_det_val = std::abs(det_val);
+  
+  if((abs_det_val < det_min) || (abs_det_val > det_max))  { return false; }
+  
+  Ym[pos<0,0>::n3] =  (Xm[pos<2,2>::n3]*Xm[pos<1,1>::n3] - Xm[pos<2,1>::n3]*Xm[pos<1,2>::n3]) / det_val;
+  Ym[pos<1,0>::n3] = -(Xm[pos<2,2>::n3]*Xm[pos<1,0>::n3] - Xm[pos<2,0>::n3]*Xm[pos<1,2>::n3]) / det_val;
+  Ym[pos<2,0>::n3] =  (Xm[pos<2,1>::n3]*Xm[pos<1,0>::n3] - Xm[pos<2,0>::n3]*Xm[pos<1,1>::n3]) / det_val;
+  
+  Ym[pos<0,1>::n3] = -(Xm[pos<2,2>::n3]*Xm[pos<0,1>::n3] - Xm[pos<2,1>::n3]*Xm[pos<0,2>::n3]) / det_val;
+  Ym[pos<1,1>::n3] =  (Xm[pos<2,2>::n3]*Xm[pos<0,0>::n3] - Xm[pos<2,0>::n3]*Xm[pos<0,2>::n3]) / det_val;
+  Ym[pos<2,1>::n3] = -(Xm[pos<2,1>::n3]*Xm[pos<0,0>::n3] - Xm[pos<2,0>::n3]*Xm[pos<0,1>::n3]) / det_val;
+  
+  Ym[pos<0,2>::n3] =  (Xm[pos<1,2>::n3]*Xm[pos<0,1>::n3] - Xm[pos<1,1>::n3]*Xm[pos<0,2>::n3]) / det_val;
+  Ym[pos<1,2>::n3] = -(Xm[pos<1,2>::n3]*Xm[pos<0,0>::n3] - Xm[pos<1,0>::n3]*Xm[pos<0,2>::n3]) / det_val;
+  Ym[pos<2,2>::n3] =  (Xm[pos<1,1>::n3]*Xm[pos<0,0>::n3] - Xm[pos<1,0>::n3]*Xm[pos<0,1>::n3]) / det_val;
+  
+  const eT check_val = Xm[pos<0,0>::n3]*Ym[pos<0,0>::n3] + Xm[pos<0,1>::n3]*Ym[pos<1,0>::n3] + Xm[pos<0,2>::n3]*Ym[pos<2,0>::n3];
+  
+  const  T max_diff  = (is_float<T>::value) ? T(1e-4) : T(1e-10);  // empirically determined; may need tuning
+  
+  if(std::abs(T(1) - check_val) >= max_diff)  { return false; }
+  
+  arrayops::copy(Xm, Ym, uword(3*3));
+  
+  return true;
+  }
+
+
+
+template<typename eT>
+arma_cold
+inline
+bool
+op_inv_gen_full::apply_tiny_4x4(Mat<eT>& X)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename get_pod_type<eT>::result T;
+  
+  // NOTE: assuming matrix X is square sized
+  
+  constexpr T det_min =        std::numeric_limits<T>::epsilon();
+  constexpr T det_max = T(1) / std::numeric_limits<T>::epsilon();
+  
+  Mat<eT> Y(4, 4, arma_nozeros_indicator());
+  
+  eT* Xm = X.memptr();
+  eT* Ym = Y.memptr();
+  
+  const eT     det_val = op_det::apply_tiny(X);
+  const  T abs_det_val = std::abs(det_val);
+  
+  if((abs_det_val < det_min) || (abs_det_val > det_max))  { return false; }
+  
+  Ym[pos<0,0>::n4] = ( Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] + Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
+  Ym[pos<1,0>::n4] = ( Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
+  Ym[pos<2,0>::n4] = ( Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] ) / det_val;
+  Ym[pos<3,0>::n4] = ( Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] ) / det_val;
+  
+  Ym[pos<0,1>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
+  Ym[pos<1,1>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
+  Ym[pos<2,1>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] ) / det_val;
+  Ym[pos<3,1>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] ) / det_val;
+  
+  Ym[pos<0,2>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
+  Ym[pos<1,2>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,3>::n4] ) / det_val;
+  Ym[pos<2,2>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,3>::n4] ) / det_val;
+  Ym[pos<3,2>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,2>::n4] ) / det_val;
+  
+  Ym[pos<0,3>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4] ) / det_val;
+  Ym[pos<1,3>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4] ) / det_val;
+  Ym[pos<2,3>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4] ) / det_val;
+  Ym[pos<3,3>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4] ) / det_val;
+  
+  const eT check_val = Xm[pos<0,0>::n4]*Ym[pos<0,0>::n4] + Xm[pos<0,1>::n4]*Ym[pos<1,0>::n4] + Xm[pos<0,2>::n4]*Ym[pos<2,0>::n4] + Xm[pos<0,3>::n4]*Ym[pos<3,0>::n4];
+  
+  const  T max_diff  = (is_float<T>::value) ? T(1e-4) : T(1e-10);  // empirically determined; may need tuning
+  
+  if(std::abs(T(1) - check_val) >= max_diff)  { return false; }
+  
+  arrayops::copy(Xm, Ym, uword(4*4));
   
   return true;
   }
diff --git a/include/armadillo_bits/op_inv_spd_bones.hpp b/include/armadillo_bits/op_inv_spd_bones.hpp
index ac716922..693504f2 100644
--- a/include/armadillo_bits/op_inv_spd_bones.hpp
+++ b/include/armadillo_bits/op_inv_spd_bones.hpp
@@ -40,6 +40,14 @@ class op_inv_spd_full
   {
   public:
   
+  template<const uword row, const uword col>
+  struct pos
+    {
+    static constexpr uword n2 = row + col*2;
+    static constexpr uword n3 = row + col*3;
+    static constexpr uword n4 = row + col*4;
+    };
+  
   template<typename T1>
   inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_inv_spd_full>& in);
   
@@ -47,7 +55,7 @@ class op_inv_spd_full
   inline static bool apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type,T1>& expr, const uword flags);
   
   template<typename eT>
-  inline static bool apply_tiny(Mat<eT>& out);
+  arma_cold inline static bool apply_tiny_2x2(Mat<eT>& X);
   };
 
 
diff --git a/include/armadillo_bits/op_inv_spd_meat.hpp b/include/armadillo_bits/op_inv_spd_meat.hpp
index de935af6..b791ef36 100644
--- a/include/armadillo_bits/op_inv_spd_meat.hpp
+++ b/include/armadillo_bits/op_inv_spd_meat.hpp
@@ -119,17 +119,30 @@ op_inv_spd_full::apply_direct(Mat<typename T1::elem_type>& out, const Base<typen
       }
     }
   
-  const uword N = (std::min)(out.n_rows, out.n_cols);
+  const uword N = out.n_rows;
   
-  if(tiny && (N <= 4) && (is_cx<eT>::no))
+  if(is_cx<eT>::no)
     {
-    arma_extra_debug_print("op_inv_spd_full: attempting tinymatrix optimisation");
+    if(N == 1)
+      {
+      const T a = access::tmp_real(out[0]);
+      
+      out[0] = eT(T(1) / a);
+      
+      return (a > T(0));
+      }
+    else
+    if(N == 2)
+      {
+      const bool status = op_inv_spd_full::apply_tiny_2x2(out);
+      
+      if(status)  { return true; }
+      }
     
-    const bool status = op_inv_spd_full::apply_tiny(out);
-    
-    if(status)  { return true; }
-    
-    arma_extra_debug_print("op_inv_spd_full: tinymatrix optimisation failed");
+    // NOTE: currently the "tiny" option is not used,
+    // NOTE: as it's cumbersome to implement optimisations for 3x3 and 4x4 matrices;
+    // NOTE: need to ensure that all leading principal minors are positive,
+    // NOTE: which would involve a lot of dedicated code
     
     // fallthrough if optimisation failed
     }
@@ -163,56 +176,45 @@ op_inv_spd_full::apply_direct(Mat<typename T1::elem_type>& out, const Base<typen
 
 
 template<typename eT>
+arma_cold
 inline
 bool
-op_inv_spd_full::apply_tiny(Mat<eT>& out)
+op_inv_spd_full::apply_tiny_2x2(Mat<eT>& X)
   {
   arma_extra_debug_sigprint();
   
   typedef typename get_pod_type<eT>::result T;
   
-  const uword N = out.n_rows;
+  // NOTE: assuming matrix X is square sized
+  // NOTE: assuming matrix X is symmetric
+  // NOTE: assuming matrix X is real
   
-  if((arma_config::debug) && (arma_config::warn_level > 0))
-    {
-    bool print_warning = false;
-    
-    T max_diag = T(0);
-    
-    const eT* colmem = out.memptr();
-    
-    for(uword i=0; i<N; ++i)
-      {
-      const eT& out_ii      = colmem[i];
-      const  T  out_ii_real = access::tmp_real(out_ii);
-      
-      // NOTE: inv_opts::tiny is also used as a workaround for broken user software
-        
-      print_warning = (out_ii_real <= T(0)) ? true : print_warning;
-      
-      max_diag = (out_ii_real > max_diag) ? out_ii_real : max_diag;
-      
-      colmem += N;
-      }
-    
-    colmem = out.memptr();
-    
-    for(uword c=0; c < N; ++c)
-      {
-      for(uword r=(c+1); r < N; ++r)
-        {
-        const T abs_val = std::abs(colmem[r]);
-        
-        print_warning = (abs_val > max_diag) ? true : print_warning;
-        }
-      
-      colmem += N;
-      }
-    
-    if(print_warning)  { arma_debug_warn_level(1, "inv_sympd(): given matrix is not positive definite"); }
-    }
+  constexpr T det_min =        std::numeric_limits<T>::epsilon();
+  constexpr T det_max = T(1) / std::numeric_limits<T>::epsilon();
   
-  return op_inv_gen_full::apply_tiny(out);
+  eT* Xm = X.memptr();
+  
+  const T a = access::tmp_real(Xm[pos<0,0>::n2]);
+  const T c = access::tmp_real(Xm[pos<1,0>::n2]);
+  const T d = access::tmp_real(Xm[pos<1,1>::n2]);
+  
+  const T det_val = (a*d - c*c);
+  
+  // positive definite iff all leading principal minors are positive
+  // a       = first  leading principal minor (top-left 1x1 submatrix)
+  // det_val = second leading principal minor (top-left 2x2 submatrix)
+  
+  if(a <= T(0))  { return false; }
+  
+  // NOTE: since det_min is positive, this also checks whether det_val is positive
+  if((det_val < det_min) || (det_val > det_max))  { return false; }
+  
+  Xm[pos<0,0>::n2] =  d / det_val;
+  Xm[pos<0,1>::n2] = -c / det_val;
+  Xm[pos<1,0>::n2] = -c / det_val;
+  Xm[pos<1,1>::n2] =  a / det_val;
+  
+  return true;
   }