From b6aa7705fd419fcb50e2e1258179fc2e45a1bb1d Mon Sep 17 00:00:00 2001
From: conrad <conrad@localhost.localdomain>
Date: Tue, 15 Mar 2022 13:23:06 +1000
Subject: [PATCH] simplification

---
 include/armadillo_bits/op_expmat_meat.hpp | 162 +++++++++++-----------
 1 file changed, 81 insertions(+), 81 deletions(-)

diff --git a/include/armadillo_bits/op_expmat_meat.hpp b/include/armadillo_bits/op_expmat_meat.hpp
index c0b88d4e..178faa14 100644
--- a/include/armadillo_bits/op_expmat_meat.hpp
+++ b/include/armadillo_bits/op_expmat_meat.hpp
@@ -65,100 +65,100 @@ op_expmat::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1
     const uword N = (std::min)(out.n_rows, out.n_cols);
     
     for(uword i=0; i<N; ++i)  { out.at(i,i) = std::exp( out.at(i,i) ); }
+    
+    return true;
     }
-  else
+  
+  Mat<eT> A = expr.get_ref();
+  
+  arma_debug_check( (A.is_square() == false), "expmat(): given matrix must be square sized" );
+  
+  if(A.is_diagmat())
     {
-    Mat<eT> A = expr.get_ref();
+    arma_extra_debug_print("op_expmat: detected diagonal matrix");
     
-    arma_debug_check( (A.is_square() == false), "expmat(): given matrix must be square sized" );
+    const uword N = (std::min)(A.n_rows, A.n_cols);
     
-    if(A.is_diagmat())
+    out.zeros(N,N);
+    
+    for(uword i=0; i<N; ++i)  { out.at(i,i) = std::exp( A.at(i,i) ); }
+    
+    return true;
+    }
+  
+  #if defined(ARMA_OPTIMISE_SYMPD)
+    const bool try_sympd = sympd_helper::guess_sympd(A);
+  #else
+    const bool try_sympd = false;
+  #endif
+  
+  if(try_sympd)
+    {
+    arma_extra_debug_print("op_expmat: attempting sympd optimisation");
+    
+    // if matrix A is sympd, all its eigenvalues are positive
+    
+    Col< T> eigval;
+    Mat<eT> eigvec;
+    
+    const bool eig_status = eig_sym_helper(eigval, eigvec, A, 'd', "expmat()");
+    
+    if(eig_status)
       {
-      arma_extra_debug_print("op_expmat: detected diagonal matrix");
+      eigval = exp(eigval);
       
-      const uword N = (std::min)(A.n_rows, A.n_cols);
-      
-      out.zeros(N,N);
-      
-      for(uword i=0; i<N; ++i)  { out.at(i,i) = std::exp( A.at(i,i) ); }
+      out = eigvec * diagmat(eigval) * eigvec.t();
       
       return true;
       }
     
-    #if defined(ARMA_OPTIMISE_SYMPD)
-      const bool try_sympd = sympd_helper::guess_sympd(A);
-    #else
-      const bool try_sympd = false;
-    #endif
+    arma_extra_debug_print("op_expmat: sympd optimisation failed");
     
-    if(try_sympd)
-      {
-      arma_extra_debug_print("op_expmat: attempting sympd optimisation");
-      
-      // if matrix A is sympd, all its eigenvalues are positive
-      
-      Col< T> eigval;
-      Mat<eT> eigvec;
-      
-      const bool eig_status = eig_sym_helper(eigval, eigvec, A, 'd', "expmat()");
-      
-      if(eig_status)
-        {
-        eigval = exp(eigval);
-        
-        out = eigvec * diagmat(eigval) * eigvec.t();
-        
-        return true;
-        }
-      
-      arma_extra_debug_print("op_expmat: sympd optimisation failed");
-      
-      // fallthrough if eigen decomposition failed
-      }
-    
-    const T norm_val = arma::norm(A, "inf");
-    
-    const double log2_val = (norm_val > T(0)) ? double(eop_aux::log2(norm_val)) : double(0);
-    
-    int exponent = int(0);  std::frexp(log2_val, &exponent);
-    
-    const uword s = uword( (std::max)(int(0), exponent + int(1)) );
-    
-    A /= eT(eop_aux::pow(double(2), double(s)));
-    
-    T c = T(0.5);
-    
-    Mat<eT> E(A.n_rows, A.n_rows, fill::eye);  E += c * A;
-    Mat<eT> D(A.n_rows, A.n_rows, fill::eye);  D -= c * A;
-    
-    Mat<eT> X = A;
-    
-    bool positive = true;
-    
-    const uword N = 6;
-    
-    for(uword i = 2; i <= N; ++i)
-      {
-      c = c * T(N - i + 1) / T(i * (2*N - i + 1));
-      
-      X = A * X;
-      
-      E += c * X;
-      
-      if(positive)  { D += c * X; }  else  { D -= c * X; }
-      
-      positive = (positive) ? false : true;
-      }
-    
-    if( (D.is_finite() == false) || (E.is_finite() == false) )  { return false; }
-    
-    const bool status = solve(out, D, E, solve_opts::no_approx);
-    
-    if(status == false)  { return false; }
-    
-    for(uword i=0; i < s; ++i)  { out = out * out; }
+    // fallthrough if eigen decomposition failed
     }
   
+  const T norm_val = arma::norm(A, "inf");
+  
+  const double log2_val = (norm_val > T(0)) ? double(eop_aux::log2(norm_val)) : double(0);
+  
+  int exponent = int(0);  std::frexp(log2_val, &exponent);
+  
+  const uword s = uword( (std::max)(int(0), exponent + int(1)) );
+  
+  A /= eT(eop_aux::pow(double(2), double(s)));
+  
+  T c = T(0.5);
+  
+  Mat<eT> E(A.n_rows, A.n_rows, fill::eye);  E += c * A;
+  Mat<eT> D(A.n_rows, A.n_rows, fill::eye);  D -= c * A;
+  
+  Mat<eT> X = A;
+  
+  bool positive = true;
+  
+  const uword N = 6;
+  
+  for(uword i = 2; i <= N; ++i)
+    {
+    c = c * T(N - i + 1) / T(i * (2*N - i + 1));
+    
+    X = A * X;
+    
+    E += c * X;
+    
+    if(positive)  { D += c * X; }  else  { D -= c * X; }
+    
+    positive = (positive) ? false : true;
+    }
+  
+  if( (D.is_finite() == false) || (E.is_finite() == false) )  { return false; }
+  
+  const bool status = solve(out, D, E, solve_opts::no_approx);
+  
+  if(status == false)  { return false; }
+  
+  for(uword i=0; i < s; ++i)  { out = out * out; }
+  
   return true;
   }