From f7e45d2d67c6cb16efa058cafa4b67dd975bca9f Mon Sep 17 00:00:00 2001
From: conrad <conrad@localhost.localdomain>
Date: Wed, 13 Apr 2022 15:00:15 +1000
Subject: [PATCH] partial implementation of op_pinv_default::apply_solve()

---
 include/armadillo_bits/op_pinv_bones.hpp |  2 +
 include/armadillo_bits/op_pinv_meat.hpp  | 87 ++++++++++++++++++++++++
 2 files changed, 89 insertions(+)

diff --git a/include/armadillo_bits/op_pinv_bones.hpp b/include/armadillo_bits/op_pinv_bones.hpp
index bf83ddc9..194a534e 100644
--- a/include/armadillo_bits/op_pinv_bones.hpp
+++ b/include/armadillo_bits/op_pinv_bones.hpp
@@ -30,6 +30,8 @@ class op_pinv_default
   template<typename T1> inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_pinv_default>& in);
   
   template<typename T1> inline static bool apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type,T1>& expr);
+  
+  // template<typename eT> inline static bool apply_solve(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B);
   };
 
 
diff --git a/include/armadillo_bits/op_pinv_meat.hpp b/include/armadillo_bits/op_pinv_meat.hpp
index 283c369f..2c86fc88 100644
--- a/include/armadillo_bits/op_pinv_meat.hpp
+++ b/include/armadillo_bits/op_pinv_meat.hpp
@@ -57,6 +57,93 @@ op_pinv_default::apply_direct(Mat<typename T1::elem_type>& out, const Base<typen
 
 
 
+// template<typename T1>
+// inline
+// bool
+// op_pinv_default::apply_solve(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B)
+//   {
+//   arma_extra_debug_sigprint();
+//   
+//   // NOTE: assuming A is square-sized
+//   // NOTE: assuming sizes of A and B are conformant for matrix multiplication
+//   
+//   typedef typename T1::pod_type T;
+//   
+//   constexpr T     tol       = T(0);
+//   constexpr uword method_id = uword(0);
+//   
+//   if(A.is_diagmat())
+//     {
+//     arma_extra_debug_print("op_pinv: detected diagonal matrix");
+//     
+//     Mat<eT> tmp;
+//     
+//     const bool status = op_pinv::apply_diag(tmp, A, tol);
+//     
+//     if(status == false)  { return false; }
+//     
+//     out = diagmat(tmp)*B;
+//     
+//     return true;
+//     }
+//   
+//   bool do_sym   = false;
+//   bool do_sympd = false;
+//   
+//   const bool is_sym_size_ok = (A.n_rows > (is_cx<eT>::yes ? uword(20) : uword(40)));
+//   
+//   if( (arma_config::optimise_sympd) && (auxlib::crippled_lapack(A) == false) && (is_sym_size_ok) )
+//     {
+//     bool is_approx_sym   = false;
+//     bool is_approx_sympd = false;
+//     
+//     sympd_helper::analyse_matrix(is_approx_sym, is_approx_sympd, A);
+//     
+//     do_sym   = is_sym_size_ok && ((is_cx<eT>::no) ? (is_approx_sym) : (is_approx_sym && is_approx_sympd));
+//     do_sympd = is_arg_default && is_approx_sympd;
+//     }
+//   
+//   if(do_sympd)
+//     {
+//     arma_extra_debug_print("op_pinv: attempting sympd optimisation");
+//     
+//     Mat<eT> tmp = A;
+//     
+//           bool is_sympd_junk   = false;
+//           T    rcond_calc      = T(0);
+//     const T    rcond_threshold = T((std::max)(uword(100), uword(A.n_rows))) * std::numeric_limits<T>::epsilon();
+//     
+//     const bool status = auxlib::solve_sympd_rcond(tmp, is_sympd_junk, rcond_calc, A, B, false);
+//     
+//     if((status) && (rcond_calc >= rcond_threshold) && (arma_isfinite(rcond_calc)))
+//       {
+//       out.steal_mem(tmp);
+//       return true;
+//       }
+//     
+//     arma_extra_debug_print("op_pinv: sympd optimisation failed");
+//     }
+//   
+//   if(do_sym)
+//     {
+//     arma_extra_debug_print("op_pinv: symmetric/hermitian optimisation");
+//     
+//     Mat<eT> tmp;
+//     
+//     const bool status = op_pinv::apply_sym(tmp, A, tol, method_id);
+//     
+//     if(status == false)  { return false; }
+//     
+//     out = tmp*B;
+//     
+//     return true;
+//     }
+//   
+//   // TODO: call auxlib::solve_approx_svd()
+//   }
+
+
+
 //