summaryrefslogblamecommitdiffstatshomepage
path: root/src/lib/netlist/solver/nld_ms_sor_mat.h
blob: f53727cfe9bd32ae7a1bf605835fa7bfdbb38766 (plain) (tree)
1
2
3
4
5
6
7
8

                               
  




                                                        


   

                         
 

                    
                                 
                                     
                              
 



                         

                                                                            
 

                                     

       
                                                                                                                      
                                                                                                                   




                                                          
                 
                 
 
                                             
 
                                                                 
 
                                                                           

        
                                                 
 



                                       


                                                                                           
                                 

                                                                                           

                                                                                
 
                                                             

 

                           

                                                           
 




                                                        
                                                        




                                                                                                                                                 
                                                        







                                                                          


                                                                         
                                                        
                 
                                                                
                                                                                            



                                                                     


                                                                                      
                                                                                 
                    
                                                   



                                             
 
      
 

                                                                                          
 



                                                                                  
 
 
                                             
                                      
 


                                                                 



                             
 
     
                              
                 
                                   






                                            
                                                                  




                                              
                                                                  

                                                
                                                                  
                                            
                                           
                                              







                                                                                     
                                                                  


         
                                         

                                                         

                                
                                     
 
                                                 
                 
                                             
 
                                                                      


                                                                         
                                                                                
 
                                                                                                 
                                                               
                                          






                                                     
                                                                       
 
                                    
                                        
 



                                                                                                                
 
                                                                                                 

              
                                   

                                   
 

 

                             

                                   
// license:GPL-2.0+
// copyright-holders:Couriersud
/*
 * nld_ms_sor.h
 *
 * Generic successive over relaxation solver.
 *
 * Fow w==1 we will do the classic Gauss-Seidel approach
 *
 */

#ifndef NLD_MS_SOR_MAT_H_
#define NLD_MS_SOR_MAT_H_

#include <algorithm>

#include "solver/nld_ms_direct.h"
#include "solver/nld_matrix_solver.h"
#include "solver/nld_solver.h"

namespace netlist
{
	namespace devices
	{
template <unsigned m_N, unsigned storage_N>
class matrix_solver_SOR_mat_t: public matrix_solver_direct_t<m_N, storage_N>
{
	friend class matrix_solver_t;

public:

	matrix_solver_SOR_mat_t(netlist_t &anetlist, const pstring &name, const solver_parameters_t *params, int size)
		: matrix_solver_direct_t<m_N, storage_N>(anetlist, name, matrix_solver_t::DESCENDING, params, size)
		, m_Vdelta(*this, "m_Vdelta", 0.0)
		, m_omega(*this, "m_omega", params->m_sor)
		, m_lp_fact(*this, "m_lp_fact", 0)
		, m_gs_fail(*this, "m_gs_fail", 0)
		, m_gs_total(*this, "m_gs_total", 0)
		{
		}

	virtual ~matrix_solver_SOR_mat_t() {}

	virtual void vsetup(analog_net_t::list_t &nets) override;

	virtual int vsolve_non_dynamic(const bool newton_raphson) override;

private:
	state_var<nl_double[storage_N]> m_Vdelta;

	state_var<nl_double> m_omega;
	state_var<nl_double> m_lp_fact;
	state_var<int> m_gs_fail;
	state_var<int> m_gs_total;
};

// ----------------------------------------------------------------------------------------
// matrix_solver - Gauss - Seidel
// ----------------------------------------------------------------------------------------

template <unsigned m_N, unsigned storage_N>
void matrix_solver_SOR_mat_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
{
	matrix_solver_direct_t<m_N, storage_N>::vsetup(nets);
}

#if 0
//FIXME: move to solve_base
template <unsigned m_N, unsigned storage_N>
nl_double matrix_solver_SOR_mat_t<m_N, storage_N>::vsolve()
{
	/*
	 * enable linear prediction on first newton pass
	 */

	if (USE_LINEAR_PREDICTION)
		for (unsigned k = 0; k < this->N(); k++)
		{
			this->m_last_V[k] = this->m_nets[k]->m_cur_Analog;
			this->m_nets[k]->m_cur_Analog = this->m_nets[k]->m_cur_Analog + this->m_Vdelta[k] * this->current_timestep() * m_lp_fact;
		}
	else
		for (unsigned k = 0; k < this->N(); k++)
		{
			this->m_last_V[k] = this->m_nets[k]->m_cur_Analog;
		}

	this->solve_base(this);

	if (USE_LINEAR_PREDICTION)
	{
		nl_double sq = 0;
		nl_double sqo = 0;
		const nl_double rez_cts = 1.0 / this->current_timestep();
		for (unsigned k = 0; k < this->N(); k++)
		{
			const analog_net_t *n = this->m_nets[k];
			const nl_double nv = (n->Q_Analog() - this->m_last_V[k]) * rez_cts ;
			sq += nv * nv;
			sqo += this->m_Vdelta[k] * this->m_Vdelta[k];
			this->m_Vdelta[k] = nv;
		}

		// FIXME: used to be 1e90, but this would not be compatible with float
		if (sqo > NL_FCONST(1e-20))
			m_lp_fact = std::min(std::sqrt(sq/sqo), (nl_double) 2.0);
		else
			m_lp_fact = NL_FCONST(0.0);
	}


	return this->compute_next_timestep();
}
#endif

template <unsigned m_N, unsigned storage_N>
int matrix_solver_SOR_mat_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
{
	/* The matrix based code looks a lot nicer but actually is 30% slower than
	 * the optimized code which works directly on the data structures.
	 * Need something like that for gaussian elimination as well.
	 */


	nl_double new_v[storage_N] = { 0.0 };
	const unsigned iN = this->N();

	matrix_solver_t::build_LE_A<matrix_solver_SOR_mat_t>();
	matrix_solver_t::build_LE_RHS<matrix_solver_SOR_mat_t>();

	bool resched = false;

	int  resched_cnt = 0;


#if 0
	static int ws_cnt = 0;
	ws_cnt++;
	if (1 && ws_cnt % 200 == 0)
	{
		// update omega
		nl_double lambdaN = 0;
		nl_double lambda1 = 1e9;
		for (int k = 0; k < iN; k++)
		{
	#if 0
			nl_double akk = std::abs(this->m_A[k][k]);
			if ( akk > lambdaN)
				lambdaN = akk;
			if (akk < lambda1)
				lambda1 = akk;
	#else
			nl_double akk = std::abs(this->m_A[k][k]);
			nl_double s = 0.0;
			for (int i=0; i<iN; i++)
				s = s + std::abs(this->m_A[k][i]);
			akk = s / akk - 1.0;
			if ( akk > lambdaN)
				lambdaN = akk;
			if (akk < lambda1)
				lambda1 = akk;
	#endif
		}
		//printf("lambda: %f %f\n", lambda, 2.0 / (1.0 + 2 * sqrt(lambda)) );

		//ws = 2.0 / (2.0 - lambdaN - lambda1);
		m_omega = 2.0 / (2.0 - lambda1);
		//printf("%f %f %f\n", m_omega, lambda1, lambdaN);
	}
#endif

	for (unsigned k = 0; k < iN; k++)
		new_v[k] = this->m_nets[k]->m_cur_Analog;

	do {
		resched = false;
		nl_double cerr = 0.0;

		for (unsigned k = 0; k < iN; k++)
		{
			nl_double Idrive = 0;

			const auto *p = this->m_terms[k]->m_nz.data();
			const unsigned e = this->m_terms[k]->m_nz.size();

			for (unsigned i = 0; i < e; i++)
				Idrive = Idrive + this->A(k,p[i]) * new_v[p[i]];

			const nl_double delta = m_omega * (this->RHS(k) - Idrive) / this->A(k,k);
			cerr = std::max(cerr, std::abs(delta));
			new_v[k] += delta;
		}

		if (cerr > this->m_params.m_accuracy)
		{
			resched = true;
		}
		resched_cnt++;
	} while (resched && (resched_cnt < this->m_params.m_gs_loops));

	this->m_stat_calculations++;
	this->m_gs_total += resched_cnt;

	if (resched)
	{
		//this->netlist().warning("Falling back to direct solver .. Consider increasing RESCHED_LOOPS");
		this->m_gs_fail++;

		return matrix_solver_direct_t<m_N, storage_N>::solve_non_dynamic(newton_raphson);
	}
	else {
		this->store(new_v);
		return resched_cnt;
	}

}

	} //namespace devices
} // namespace netlist

#endif /* NLD_MS_GAUSS_SEIDEL_H_ */