1 // license:GPL-2.0+ 2 // copyright-holders:Couriersud 3 4 #ifndef NLD_MS_GMRES_H_ 5 #define NLD_MS_GMRES_H_ 6 7 /// 8 /// \file nld_ms_gmres.h 9 /// 10 11 #include "nld_matrix_solver_ext.h" 12 #include "nld_ms_direct.h" 13 #include "nld_solver.h" 14 #include "plib/gmres.h" 15 #include "plib/parray.h" 16 #include "plib/pmatrix_cr.h" 17 #include "plib/vector_ops.h" 18 19 #include <algorithm> 20 21 namespace netlist 22 { 23 namespace solver 24 { 25 26 template <typename FT, int SIZE> 27 class matrix_solver_GMRES_t: public matrix_solver_direct_t<FT, SIZE> 28 { 29 public: 30 31 using float_type = FT; 32 33 // Sort rows in ascending order. This should minimize fill-in and thus 34 // maximize the efficiency of the incomplete LUT. 35 // This is already preconditioning. 36 matrix_solver_GMRES_t(devices::nld_solver & main_solver,const pstring & name,matrix_solver_t::net_list_t & nets,const solver::solver_parameters_t * params,const std::size_t size)37 matrix_solver_GMRES_t(devices::nld_solver &main_solver, const pstring &name, 38 matrix_solver_t::net_list_t &nets, 39 const solver::solver_parameters_t *params, 40 const std::size_t size) 41 : matrix_solver_direct_t<FT, SIZE>(main_solver, name, nets, params, size) 42 , m_ops(size, 0) 43 , m_gmres(size) 44 { 45 const std::size_t iN = this->size(); 46 47 std::vector<std::vector<unsigned>> fill(iN); 48 49 for (std::size_t k=0; k<iN; k++) 50 { 51 fill[k].resize(iN, decltype(m_ops.m_mat)::FILL_INFINITY); 52 terms_for_net_t & row = this->m_terms[k]; 53 for (const auto &nz_j : row.m_nz) 54 { 55 fill[k][static_cast<mattype>(nz_j)] = 0; 56 } 57 } 58 59 m_ops.build(fill); 60 this->log_fill(fill, m_ops.m_mat); 61 62 // build pointers into the compressed row format matrix for each terminal 63 64 for (std::size_t k=0; k<iN; k++) 65 { 66 std::size_t cnt = 0; 67 for (std::size_t j=0; j< this->m_terms[k].railstart();j++) 68 { 69 for (std::size_t i = m_ops.m_mat.row_idx[k]; i<m_ops.m_mat.row_idx[k+1]; i++) 70 if (this->m_terms[k].m_connected_net_idx[j] == static_cast<int>(m_ops.m_mat.col_idx[i])) 71 { 72 this->m_mat_ptr[k][j] = &m_ops.m_mat.A[i]; 73 cnt++; 74 break; 75 } 76 } 77 nl_assert(cnt == this->m_terms[k].railstart()); 78 this->m_mat_ptr[k][this->m_terms[k].railstart()] = &m_ops.m_mat.A[m_ops.m_mat.diag[k]]; 79 } 80 } 81 82 void vsolve_non_dynamic() override; 83 84 private: 85 86 using mattype = typename plib::pmatrix_cr<FT, SIZE>::index_type; 87 88 //plib::mat_precondition_none<FT, SIZE> m_ops; 89 plib::mat_precondition_ILU<FT, SIZE> m_ops; 90 //plib::mat_precondition_diag<FT, SIZE> m_ops; 91 plib::gmres_t<FT, SIZE> m_gmres; 92 }; 93 94 // ---------------------------------------------------------------------------------------- 95 // matrix_solver - GMRES 96 // ---------------------------------------------------------------------------------------- 97 98 template <typename FT, int SIZE> vsolve_non_dynamic()99 void matrix_solver_GMRES_t<FT, SIZE>::vsolve_non_dynamic() 100 { 101 const std::size_t iN = this->size(); 102 103 m_ops.m_mat.set_scalar(plib::constants<FT>::zero()); 104 105 // populate matrix and V for first estimate 106 this->fill_matrix_and_rhs(); 107 108 for (std::size_t k = 0; k < iN; k++) 109 { 110 this->m_new_V[k] = static_cast<float_type>(this->m_terms[k].getV()); 111 } 112 113 const auto accuracy(static_cast<float_type>(this->m_params.m_accuracy)); 114 115 auto iter = std::max(plib::constants<std::size_t>::one(), this->m_params.m_gs_loops()); 116 auto gsl = m_gmres.solve(m_ops, this->m_new_V, this->m_RHS, iter, accuracy); 117 118 this->m_iterative_total += gsl; 119 120 if (gsl > iter) 121 { 122 this->m_iterative_fail++; 123 matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(); 124 } 125 126 } 127 128 129 130 131 } // namespace solver 132 } // namespace netlist 133 134 #endif // NLD_MS_GMRES_H_ 135