1 files changed, 138 insertions, 0 deletions
diff --git a/src/armadillo/include/armadillo_bits/newarp_SparseGenRealShiftSolve_meat.hpp b/src/armadillo/include/armadillo_bits/newarp_SparseGenRealShiftSolve_meat.hpp
new file mode 100644
index 0000000..ea20618
--- /dev/null
+++ b/src/armadillo/include/armadillo_bits/newarp_SparseGenRealShiftSolve_meat.hpp
@@ -0,0 +1,138 @@
+// SPDX-License-Identifier: Apache-2.0
+// 
+// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
+// Copyright 2008-2016 National ICT Australia (NICTA)
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+// ------------------------------------------------------------------------
+
+
+namespace newarp
+{
+
+
+template<typename eT>
+inline
+SparseGenRealShiftSolve<eT>::SparseGenRealShiftSolve(const SpMat<eT>& mat_obj, const eT shift)
+  #if defined(ARMA_USE_SUPERLU)
+    : perm_c(mat_obj.n_cols + 1)
+    , perm_r(mat_obj.n_rows + 1)
+    , n_rows(mat_obj.n_rows)
+    , n_cols(mat_obj.n_cols)
+  #else
+    : n_rows(0)
+    , n_cols(0)
+  #endif
+  {
+  arma_extra_debug_sigprint();
+  
+  #if defined(ARMA_USE_SUPERLU)
+    {
+    // Derived from sp_auxlib::run_aupd_shiftinvert()
+    superlu_opts superlu_opts_default;
+    superlu::superlu_options_t options;
+    sp_auxlib::set_superlu_opts(options, superlu_opts_default);
+    
+    superlu::GlobalLU_t Glu;
+    arrayops::fill_zeros(reinterpret_cast<char*>(&Glu), sizeof(superlu::GlobalLU_t));
+    
+    superlu_supermatrix_wrangler x;
+    superlu_supermatrix_wrangler xC;
+    superlu_array_wrangler<int> etree(mat_obj.n_cols+1);
+    
+    // Copy A-shift*I to x
+    const bool status_x = sp_auxlib::copy_to_supermatrix_with_shift(x.get_ref(), mat_obj, shift);
+    
+    if(status_x == false)  { arma_stop_runtime_error("newarp::SparseGenRealShiftSolve::SparseGenRealShiftSolve(): could not construct SuperLU matrix"); return; }
+    
+    int panel_size = superlu::sp_ispec_environ(1);
+    int relax      = superlu::sp_ispec_environ(2);
+    int slu_info   = 0; // Return code
+    int lwork      = 0; // lwork = 0: allocate space internally by system malloc
+    
+    superlu_stat_wrangler stat;
+    
+    arma_extra_debug_print("superlu::gstrf()");
+    superlu::get_permutation_c(options.ColPerm, x.get_ptr(), perm_c.get_ptr());
+    superlu::sp_preorder_mat(&options, x.get_ptr(), perm_c.get_ptr(), etree.get_ptr(), xC.get_ptr());
+    superlu::gstrf<eT>(&options, xC.get_ptr(), relax, panel_size, etree.get_ptr(), NULL, lwork, perm_c.get_ptr(), perm_r.get_ptr(), l.get_ptr(), u.get_ptr(), &Glu, stat.get_ptr(), &slu_info);
+    
+    if(slu_info != 0)
+      {
+      arma_debug_warn_level(2, "matrix is singular to working precision");
+      return;
+      }
+    
+    eT x_norm_val = sp_auxlib::norm1<eT>(x.get_ptr());
+    eT x_rcond    = sp_auxlib::lu_rcond<eT>(l.get_ptr(), u.get_ptr(), x_norm_val);
+    
+    if( (x_rcond < std::numeric_limits<eT>::epsilon()) || arma_isnan(x_rcond) )
+      {
+      if(x_rcond == eT(0))  { arma_debug_warn_level(2, "matrix is singular to working precision");                        }
+      else                  { arma_debug_warn_level(2, "matrix is singular to working precision (rcond: ", x_rcond, ")"); }
+      return;
+      }
+    
+    valid = true;
+    }
+  #else
+    {
+    arma_ignore(mat_obj);
+    arma_ignore(shift);
+    }
+  #endif
+  }
+
+
+
+// Perform the shift-solve operation \f$y=(A-\sigma I)^{-1}x\f$.
+// y_out = inv(A - sigma * I) * x_in
+template<typename eT>
+inline
+void
+SparseGenRealShiftSolve<eT>::perform_op(eT* x_in, eT* y_out) const
+  {
+  arma_extra_debug_sigprint();
+  
+  #if defined(ARMA_USE_SUPERLU)
+    {
+    const Col<eT> x(x_in , n_cols, false, true);
+          Col<eT> y(y_out, n_rows, false, true);
+    
+    // Derived from sp_auxlib::run_aupd_shiftinvert()
+    y = x;
+    superlu_supermatrix_wrangler out_slu;
+    
+    const bool status_out_slu = sp_auxlib::wrap_to_supermatrix(out_slu.get_ref(), y);
+    
+    if(status_out_slu == false)  { arma_stop_runtime_error("newarp::SparseGenRealShiftSolve::perform_op(): could not construct SuperLU matrix"); return; }
+    
+    superlu_stat_wrangler stat;
+    int info = 0;
+    
+    arma_extra_debug_print("superlu::gstrs()");
+    superlu::gstrs<eT>(superlu::NOTRANS, l.get_ptr(), u.get_ptr(), perm_c.get_ptr(), perm_r.get_ptr(), out_slu.get_ptr(), stat.get_ptr(), &info);
+    
+    if(info != 0)  { arma_stop_runtime_error("newarp::SparseGenRealShiftSolve::perform_op(): could not solve linear equation"); return; }
+    
+    // No need to modify memory further since it was all done in-place.
+    }
+  #else
+    {
+    arma_ignore(x_in);
+    arma_ignore(y_out);
+    }
+  #endif
+  }
+
+
+}  // namespace newarp