1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19 namespace newarp
20 {
21
22
23 template<typename eT>
24 inline
SparseGenRealShiftSolve(const SpMat<eT> & mat_obj,const eT shift)25 SparseGenRealShiftSolve<eT>::SparseGenRealShiftSolve(const SpMat<eT>& mat_obj, const eT shift)
26 #if defined(ARMA_USE_SUPERLU)
27 : perm_c(mat_obj.n_cols + 1)
28 , perm_r(mat_obj.n_rows + 1)
29 , n_rows(mat_obj.n_rows)
30 , n_cols(mat_obj.n_cols)
31 #endif
32 {
33 arma_extra_debug_sigprint();
34
35 #if defined(ARMA_USE_SUPERLU)
36 {
37 // Derived from sp_auxlib::run_aupd_shiftinvert()
38 superlu_opts superlu_opts_default;
39 superlu::superlu_options_t options;
40 sp_auxlib::set_superlu_opts(options, superlu_opts_default);
41
42 superlu::GlobalLU_t Glu;
43 arrayops::fill_zeros(reinterpret_cast<char*>(&Glu), sizeof(superlu::GlobalLU_t));
44
45 superlu_supermatrix_wrangler x;
46 superlu_supermatrix_wrangler xC;
47 superlu_array_wrangler<int> etree(mat_obj.n_cols+1);
48
49 // Copy A-shift*I to x
50 const bool status_x = sp_auxlib::copy_to_supermatrix_with_shift(x.get_ref(), mat_obj, shift);
51
52 if(status_x == false) { arma_stop_runtime_error("newarp::SparseGenRealShiftSolve::SparseGenRealShiftSolve(): could not construct SuperLU matrix"); return; }
53
54 int panel_size = superlu::sp_ispec_environ(1);
55 int relax = superlu::sp_ispec_environ(2);
56 int slu_info = 0; // Return code
57 int lwork = 0; // lwork = 0: allocate space internally by system malloc
58
59 superlu_stat_wrangler stat;
60
61 arma_extra_debug_print("superlu::gstrf()");
62 superlu::get_permutation_c(options.ColPerm, x.get_ptr(), perm_c.get_ptr());
63 superlu::sp_preorder_mat(&options, x.get_ptr(), perm_c.get_ptr(), etree.get_ptr(), xC.get_ptr());
64 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);
65
66 if(slu_info != 0)
67 {
68 arma_debug_warn_level(2, "matrix is singular to working precision");
69 return;
70 }
71
72 eT x_norm_val = sp_auxlib::norm1<eT>(x.get_ptr());
73 eT x_rcond = sp_auxlib::lu_rcond<eT>(l.get_ptr(), u.get_ptr(), x_norm_val);
74
75 if( (x_rcond < std::numeric_limits<eT>::epsilon()) || arma_isnan(x_rcond) )
76 {
77 if(x_rcond > eT(0)) { arma_debug_warn_level(2, "matrix is singular to working precision (rcond: ", x_rcond, ")"); }
78 else { arma_debug_warn_level(2, "matrix is singular to working precision"); }
79 return;
80 }
81
82 valid = true;
83 }
84 #else
85 {
86 arma_ignore(mat_obj);
87 arma_ignore(shift);
88 }
89 #endif
90 }
91
92
93
94 // Perform the shift-solve operation \f$y=(A-\sigma I)^{-1}x\f$.
95 // y_out = inv(A - sigma * I) * x_in
96 template<typename eT>
97 inline
98 void
perform_op(eT * x_in,eT * y_out) const99 SparseGenRealShiftSolve<eT>::perform_op(eT* x_in, eT* y_out) const
100 {
101 arma_extra_debug_sigprint();
102
103 #if defined(ARMA_USE_SUPERLU)
104 {
105 const Col<eT> x(x_in , n_cols, false, true);
106 Col<eT> y(y_out, n_rows, false, true);
107
108 // Derived from sp_auxlib::run_aupd_shiftinvert()
109 y = x;
110 superlu_supermatrix_wrangler out_slu;
111
112 const bool status_out_slu = sp_auxlib::wrap_to_supermatrix(out_slu.get_ref(), y);
113
114 if(status_out_slu == false) { arma_stop_runtime_error("newarp::SparseGenRealShiftSolve::perform_op(): could not construct SuperLU matrix"); return; }
115
116 superlu_stat_wrangler stat;
117 int info = 0;
118
119 arma_extra_debug_print("superlu::gstrs()");
120 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);
121
122 if(info != 0) { arma_stop_runtime_error("newarp::SparseGenRealShiftSolve::perform_op(): could not solve linear equation"); return; }
123
124 // No need to modify memory further since it was all done in-place.
125 }
126 #else
127 {
128 arma_ignore(x_in);
129 arma_ignore(y_out);
130 }
131 #endif
132 }
133
134
135 } // namespace newarp
136