src/s_dc.cc

   1 /*                             -*- C++ -*-
   2  * Copyright (C) 2001 Albert Davis
   3  * Author: Albert Davis <aldavis@gnu.org>
   4  *
   5  * This file is part of "Gnucap", the Gnu Circuit Analysis Package
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  20  * 02110-1301, USA.
  21  *------------------------------------------------------------------
  22  * dc analysis top
  23  */
  24 #include "globals.h"
  25 #include "u_status.h"
  26 #include "u_prblst.h"
  27 #include "u_cardst.h"
  28 #include "e_elemnt.h"
  29 #include "s__.h"
  30 #include "io_matrix.h"
  31 /*--------------------------------------------------------------------------*/
  32 namespace { //
  33 /*--------------------------------------------------------------------------*/
  34 class DCOP : public SIM { //
  35 protected:
  36   enum STEP_CAUSE { //
  37     scUSER      =  1,   /* user requested                               */
  38     scSKIP      =  3,   /* effect of "skip" parameter                   */
  39     scITER_R    =  4,   /* iter count exceeds itl4 (reducing)           */
  40     scITER_A    =  5,   /* iter count exceeds itl3 (holding)            */
  41     scTE        =  6,   /* truncation error, or device stuff            */
  42     scINITIAL   =  9,   /* initial guess                                */
  43     scREJECT    = 10,   /* rejected previous time step                  */
  44     scZERO      = 11,   /* fixed zero time step                         */
  45     scSMALL     = 12,   /* time step too small                          */
  46     scNO_ADVANCE= 13,   /* after all that it still didn't advance       */
  47     scLAST      = 15    /* last step */
  48   };
  49   void  set_step_cause(STEP_CAUSE C) {::status.control = C;}
  50   STEP_CAUSE step_cause()const {return STEP_CAUSE(::status.control);}
  51 public:
  52   void  finish();
  53 protected:
  54   void  fix_args(int);
  55   void  options(CS&, int);
  56   std::string label()const{untested(); return "dc";}
  57 private:
  58   void  sweep();
  59   void  sweep_recursive(int);
  60   void  first(int);
  61   bool  next(int);
  62   explicit DCOP(const DCOP&): SIM() {unreachable(); incomplete();}
  63 protected:
  64   explicit DCOP();
  65   ~DCOP() {}
  66
  67 protected:
  68   enum {DCNEST = 4};
  69   int _n_sweeps;
  70   PARAMETER<double> _start[DCNEST];
  71   PARAMETER<double> _stop[DCNEST];
  72   PARAMETER<double> _step_in[DCNEST];
  73   double _step[DCNEST];
  74   double _val_by_user_request[DCNEST];
  75   double _sweepdamp[DCNEST];
  76   bool _linswp[DCNEST];
  77   double* (_sweepval[DCNEST]);  /* pointer to thing to sweep, dc command */
  78   ELEMENT* (_zap[DCNEST]);      /* to branch to zap, for re-expand */
  79   CARDSTASH _stash[DCNEST];     /* store std values of elements being swept */
  80   bool _loop[DCNEST];           /* flag: do it again backwards */
  81   bool _reverse_in[DCNEST];     /* flag: sweep backwards, input */
  82   bool _reverse[DCNEST];        /* flag: sweep backwards, working */
  83   bool _cont;                   /* flag: continue from previous run */
  84   TRACE _trace;                 /* enum: show extended diagnostics */
  85   bool _converged;
  86   bool _ever_converged;         /* don't try to step back otherwise... */
  87   enum {ONE_PT, LIN_STEP, LIN_PTS, TIMES, OCTAVE, DECADE} _stepmode[DCNEST];
  88 };
  89 /*--------------------------------------------------------------------------*/
  90 class DC : public DCOP { //
  91 public:
  92   explicit DC(): DCOP() {}
  93   ~DC() {}
  94   void  do_it(CS&, CARD_LIST*);
  95 private:
  96   void  setup(CS&);
  97   explicit DC(const DC&): DCOP() {unreachable(); incomplete();}
  98 };
  99 /*--------------------------------------------------------------------------*/
 100 class OP : public DCOP { //
 101 public:
 102   explicit OP(): DCOP() {}
 103   ~OP() {}
 104   void  do_it(CS&, CARD_LIST*);
 105 private:
 106   void  setup(CS&);
 107   explicit OP(const OP&): DCOP() {unreachable(); incomplete();}
 108 };
 109 /*--------------------------------------------------------------------------*/
 110 /*--------------------------------------------------------------------------*/
 111 void DC::do_it(CS& Cmd, CARD_LIST* Scope)
 112 {
 113   _scope = Scope;
 114   _sim->_time0 = 0.;
 115   _sim->set_command_dc();
 116   _sim->_phase = p_INIT_DC;
 117   ::status.dc.reset().start();
 118   command_base(Cmd);
 119   _sim->_has_op = s_DC;
 120   _scope = NULL;
 121   ::status.dc.stop();
 122 }
 123 /*--------------------------------------------------------------------------*/
 124 void OP::do_it(CS& Cmd, CARD_LIST* Scope)
 125 {
 126   _scope = Scope;
 127   _sim->_time0 = 0.;
 128   _sim->set_command_op();
 129   _sim->_phase = p_INIT_DC;
 130   ::status.op.reset().start();
 131   command_base(Cmd);
 132   _sim->_has_op = s_OP;
 133   _scope = NULL;
 134   ::status.op.stop();
 135 }
 136 /*--------------------------------------------------------------------------*/
 137 /*--------------------------------------------------------------------------*/
 138 DCOP::DCOP()
 139   :SIM(),
 140    _n_sweeps(1),
 141    _cont(false),
 142    _trace(tNONE),
 143    _converged(0)
 144 {
 145   for (int ii = 0; ii < DCNEST; ++ii) {
 146     _loop[ii] = false;
 147     _reverse_in[ii] = false;
 148     _reverse[ii] = false;
 149     _step[ii]=0.;
 150     _linswp[ii]=true;
 151     _sweepval[ii]=&_sim->_genout;
 152     _zap[ii]=NULL;
 153     _stepmode[ii] = ONE_PT;
 154   }
 155
 156   //BUG// in SIM.  should be initialized there.
 157   //_sim->_genout=0.;
 158   _out=IO::mstdout;
 159 }
 160 /*--------------------------------------------------------------------------*/
 161 void DCOP::finish(void)
 162 {
 163   for (int ii = 0;  ii < _n_sweeps;  ++ii) {
 164     if (_zap[ii]) { // component
 165       _stash[ii].restore();
 166       _zap[ii]->dec_probes();
 167       _zap[ii]->precalc_first();
 168       _zap[ii]->precalc_last();
 169     }else{
 170     }
 171   }
 172 }
 173 /*--------------------------------------------------------------------------*/
 174 void OP::setup(CS& Cmd)
 175 {
 176   _sim->_temp_c = OPT::temp_c;
 177   _cont = false;
 178   _trace = tNONE;
 179   _out = IO::mstdout;
 180   _out.reset(); //BUG// don't know why this is needed */
 181   bool ploton = IO::plotset  &&  plotlist().size() > 0;
 182
 183   _zap[0] = NULL;
 184   _sweepval[0] = &(_sim->_temp_c);
 185
 186   if (Cmd.match1("'\"({") || Cmd.is_float()) {
 187     Cmd >> _start[0];
 188     if (Cmd.match1("'\"({") || Cmd.is_float()) {
 189       Cmd >> _stop[0];
 190     }else{
 191       _stop[0] = _start[0];
 192     }
 193   }else{
 194   }
 195
 196   _step[0] = 0.;
 197   _sim->_genout = 0.;
 198
 199   options(Cmd,0);
 200 //  _sim->_temp_c = temp_c_in;
 201
 202   _n_sweeps = 1;
 203   Cmd.check(bWARNING, "what's this?");
 204   _sim->_freq = 0;
 205
 206   IO::plotout = (ploton) ? IO::mstdout : OMSTREAM();
 207   initio(_out);
 208
 209   _start[0].e_val(OPT::temp_c, _scope);
 210   fix_args(0);
 211 }
 212 /*--------------------------------------------------------------------------*/
 213 void DC::setup(CS& Cmd)
 214 {
 215   _sim->_temp_c = OPT::temp_c;
 216   _cont = false;
 217   _trace = tNONE;
 218   _out = IO::mstdout;
 219   _out.reset(); //BUG// don't know why this is needed */
 220   bool ploton = IO::plotset  &&  plotlist().size() > 0;
 221
 222   if (Cmd.more()) {
 223     for (_n_sweeps = 0; Cmd.more() && _n_sweeps < DCNEST; ++_n_sweeps) {
 224       CARD_LIST::fat_iterator ci = findbranch(Cmd, &CARD_LIST::card_list);
 225       if (!ci.is_end()) {                       // sweep a component
 226         if (ELEMENT* c = dynamic_cast<ELEMENT*>(*ci)) {
 227           _zap[_n_sweeps] = c;
 228         }else{ untested();
 229           throw Exception("dc/op: can't sweep " + (**ci).long_label() + '\n');
 230         }
 231       }else if (Cmd.is_float()) {               // sweep the generator
 232         _zap[_n_sweeps] = NULL;
 233       }else{
 234         // leave as it was .. repeat Cmd with no args
 235       }
 236
 237       if (Cmd.match1("'\"({") || Cmd.is_float()) {      // set up parameters
 238         _start[_n_sweeps] = "NA";
 239         _stop[_n_sweeps] = "NA";
 240         Cmd >> _start[_n_sweeps] >> _stop[_n_sweeps];
 241         _step[_n_sweeps] = 0.;
 242       }else{
 243         // leave it as it was .. repeat Cmd with no args
 244       }
 245
 246       _sim->_genout = 0.;
 247       options(Cmd,_n_sweeps);
 248 //      _sim->_temp_c = temp_c_in;
 249     }
 250   }else{
 251   }
 252   Cmd.check(bWARNING, "what's this?");
 253
 254   IO::plotout = (ploton) ? IO::mstdout : OMSTREAM();
 255   initio(_out);
 256
 257   assert(_n_sweeps > 0);
 258   for (int ii = 0;  ii < _n_sweeps;  ++ii) {
 259     _start[ii].e_val(0., _scope);
 260     fix_args(ii);
 261
 262     if (_zap[ii]) { // component
 263       _stash[ii] = _zap[ii];                    // stash the std value
 264       _zap[ii]->inc_probes();                   // we need to keep track of it
 265       _zap[ii]->set_value(_zap[ii]->value(),0); // zap out extensions
 266       _zap[ii]->set_constant(false);            // so it will be updated
 267       _sweepval[ii] = _zap[ii]->set__value();   // point to value to patch
 268     }else{ // generator
 269       _sweepval[ii] = &_sim->_genout;                   // point to value to patch
 270     }
 271   }
 272   _sim->_freq = 0;
 273 }
 274 /*--------------------------------------------------------------------------*/
 275 void DCOP::fix_args(int Nest)
 276 {
 277   _stop[Nest].e_val(_start[Nest], _scope);
 278   _step_in[Nest].e_val(0., _scope);
 279   _step[Nest] = _step_in[Nest];
 280
 281   switch (_stepmode[Nest]) { untested();
 282   case ONE_PT:
 283   case LIN_STEP:
 284     _linswp[Nest] = true;
 285     break;
 286   case LIN_PTS:untested();
 287     if (_step[Nest] <= 2.) { untested();
 288       _step[Nest] = 2.;
 289     }else{ untested();
 290     }
 291     _linswp[Nest] = true;
 292     break;
 293   case TIMES:
 294     if (_step[Nest] == 0.  &&  _start[Nest] != 0.) { untested();
 295       _step[Nest] = _stop[Nest] / _start[Nest];
 296     }else{
 297     }
 298     _linswp[Nest] = false;
 299     break;
 300   case OCTAVE:untested();
 301     if (_step[Nest] == 0.) {untested();
 302       _step[Nest] = 1.;
 303     }else{untested();
 304     }
 305     _step[Nest] = pow(2.00000001, 1./_step[Nest]);
 306     _linswp[Nest] = false;
 307     break;
 308   case DECADE:
 309     if (_step[Nest] == 0.) {untested();
 310       _step[Nest] = 1.;
 311     }else{
 312     }
 313     _step[Nest] = pow(10., 1./_step[Nest]);
 314     _linswp[Nest] = false;
 315     break;
 316   };
 317
 318   if (_step[Nest] == 0.) {      // prohibit log sweep from 0
 319     _step[Nest] = _stop[Nest] - _start[Nest];
 320     _linswp[Nest] = true;
 321   }else{
 322   }
 323 }
 324 /*--------------------------------------------------------------------------*/
 325 void DCOP::options(CS& Cmd, int Nest)
 326 {
 327   bool _dump_matrix=false;
 328   _sim->_uic = _loop[Nest] = _reverse_in[Nest] = false;
 329   unsigned here = Cmd.cursor();
 330   do{
 331     ONE_OF
 332       || (Cmd.match1("'\"({")   && ((Cmd >> _step_in[Nest]), (_stepmode[Nest] = LIN_STEP)))
 333       || (Cmd.is_float()        && ((Cmd >> _step_in[Nest]), (_stepmode[Nest] = LIN_STEP)))
 334       || (Get(Cmd, "*",           &_step_in[Nest]) && (_stepmode[Nest] = TIMES))
 335       || (Get(Cmd, "+",           &_step_in[Nest]) && (_stepmode[Nest] = LIN_STEP))
 336       || Get(Cmd, "dm",           &_dump_matrix)
 337       || (Get(Cmd, "by",          &_step_in[Nest]) && (_stepmode[Nest] = LIN_STEP))
 338       || (Get(Cmd, "step",        &_step_in[Nest]) && (_stepmode[Nest] = LIN_STEP))
 339       || (Get(Cmd, "d{ecade}",    &_step_in[Nest]) && (_stepmode[Nest] = DECADE))
 340       || (Get(Cmd, "ti{mes}",     &_step_in[Nest]) && (_stepmode[Nest] = TIMES))
 341       || (Get(Cmd, "lin",         &_step_in[Nest]) && (_stepmode[Nest] = LIN_PTS))
 342       || (Get(Cmd, "o{ctave}",    &_step_in[Nest]) && (_stepmode[Nest] = OCTAVE))
 343       || Get(Cmd, "c{ontinue}",   &_cont)
 344       || Get(Cmd, "uic",          &_sim->_uic)
 345       || Get(Cmd, "dt{emp}",      &(_sim->_temp_c),   mOFFSET, OPT::temp_c)
 346       || Get(Cmd, "lo{op}",       &_loop[Nest])
 347       || Get(Cmd, "re{verse}",    &_reverse_in[Nest])
 348       || Get(Cmd, "te{mperature}",&(_sim->_temp_c))
 349       || (Cmd.umatch("tr{ace} {=}") &&
 350           (ONE_OF
 351            || Set(Cmd, "n{one}",      &_trace, tNONE)
 352            || Set(Cmd, "o{ff}",       &_trace, tNONE)
 353            || Set(Cmd, "w{arnings}",  &_trace, tUNDER)
 354            || Set(Cmd, "a{lltime}",   &_trace, tALLTIME)
 355            || Set(Cmd, "r{ejected}",  &_trace, tREJECTED)
 356            || Set(Cmd, "i{terations}",&_trace, tITERATION)
 357            || Set(Cmd, "v{erbose}",   &_trace, tVERBOSE)
 358            || Cmd.warn(bWARNING,
 359                        "need none, off, warnings, iterations, verbose")
 360            )
 361           )
 362       || _out.outset(Cmd);
 363   }while (Cmd.more() && !Cmd.stuck(&here));
 364
 365   if(_dump_matrix) {
 366     _trace = (TRACE) (_trace | (int)tMATRIX);
 367   }else{
 368   }
 369 }
 370 /*--------------------------------------------------------------------------*/
 371 void DCOP::sweep()
 372 {
 373 //  _sim->_age= true; // hack?
 374   head(_start[0], _stop[0], " ");
 375   _sim->_bypass_ok = false;
 376   _sim->set_inc_mode_bad();
 377   if (_cont) {
 378     _sim->restore_voltages();
 379     CARD_LIST::card_list.tr_restore();
 380   }else{
 381     _sim->clear_limit();
 382     CARD_LIST::card_list.tr_begin();
 383   }
 384
 385   set_step_cause(scUSER);
 386   _converged = false;
 387   _ever_converged = false;
 388   for (int ii = 0; ii < _n_sweeps; ++ii) {
 389     if (!_zap[ii]) {
 390     }else if (_zap[ii]->is_constant()) {
 391       CARD_LIST::card_list.q_hack(_zap[ii]);
 392     }else{
 393     }
 394   }
 395   sweep_recursive(_n_sweeps-1);
 396   _sim->pop_voltages();
 397   _sim->_has_op = _sim->_mode;
 398   _sim->keep_voltages();
 399 }
 400 /*--------------------------------------------------------------------------*/
 401 static double mul(double a,double b){ return a*b; }
 402 static double sub(double a,double b){ return a-b; }
 403 static double div(double a,double b){ return a/b; }
 404 static double add(double a,double b){ return a+b; }
 405 static bool ge(double a,double b){ return a>=b; }
 406 static bool le(double a,double b){ return a<=b; }
 407 /*--------------------------------------------------------------------------*/
 408 void DCOP::sweep_recursive(int Nest)
 409 {
 410   for (int ii = int(_sim->_lu.size()); ii >= 0; --ii) {
 411     assert(!_sim->_nstat[ii].discont());
 412   }
 413   static unsigned extra_steps;
 414   assert(Nest >= 0);
 415   assert(Nest < DCNEST);
 416
 417   OPT::ITL itl = OPT::DCBIAS;
 418
 419   first(Nest);
 420
 421   double (*step)(double a, double b) = add;
 422   double (*back)(double a, double b) = sub;
 423   if (!_linswp[Nest]) {
 424     step = mul;
 425     back = div;
 426   }
 427   if (_reverse[Nest]) {
 428     std::swap(step,back);
 429   }
 430
 431   trace3("DCOP::sweep_recursive", Nest, *(_sweepval[Nest]), _step[Nest]);
 432
 433   bool firstloop=true;
 434   do {
 435 //     _sim->_temp_c = temp_c_in;
 436     if (Nest) {
 437       sweep_recursive(Nest-1);
 438       if(_converged || !_ever_converged){
 439         //if(!firstloop)
 440         _sim->pop_voltages();
 441         _sim->restore_voltages();
 442         if(firstloop) {
 443           _sim->_has_op = _sim->_mode;
 444           _sim->keep_voltages(true); // push values for sweep. hmmm
 445         }
 446       } else { untested();
 447         trace2("not converged II", Nest, *(_sweepval[Nest]));
 448         // step back...
 449         incomplete();
 450       }
 451     }else{ // leaf
 452       if (_sim->command_is_op()) {
 453         CARD_LIST::card_list.precalc_last();
 454       }else{
 455       }
 456       for (int ii = 0;  ii < _n_sweeps;  ++ii) {
 457         if (_zap[ii]) {
 458           _zap[ii]->do_tr();
 459         }
 460       }
 461       _converged = solve_with_homotopy(itl,_trace);
 462       _ever_converged |= _converged;
 463       ++::status.hidden_steps;
 464       bool printnow =
 465         (_trace >= tREJECTED)
 466         || firstloop
 467         || (_converged && ((_trace >= tALLTIME)
 468         || (step_cause() == scUSER )));
 469       trace4("outdata?", _trace, _converged, printnow, _ever_converged);
 470       if (!printnow) {
 471         ++extra_steps;
 472         if(extra_steps > 100){ untested();
 473           throw Exception("dc stepping did not succeed");
 474         }
 475       }
 476       if (!_converged) {
 477         trace3("DCOP::sweep_recursive noconv", Nest, *_sweepval[Nest], _ever_converged);
 478         error(bWARNING, "did not converge\n");
 479         _sim->restore_voltages();
 480       }else if (!_ever_converged) { incomplete();
 481         if (firstloop) {
 482           _sim->keep_voltages(true); // push values for sweep. hmmm
 483         }
 484       }else{
 485         ::status.accept.start();
 486         _sim->set_limit();
 487         CARD_LIST::card_list.tr_accept();
 488         ::status.accept.stop();
 489         _sim->keep_voltages();
 490         if (firstloop) {
 491           _sim->_has_op = _sim->_mode;
 492           _sim->keep_voltages(true); // push values for sweep. hmmm
 493         }
 494       }
 495       if (printnow) {
 496         extra_steps = 0;
 497         fixzero(_sweepval[Nest], _step[Nest]); // hack
 498         if (_converged){
 499           outdata(*_sweepval[Nest]);
 500         } else {
 501           outdata(- *_sweepval[Nest]);
 502         }
 503         ::status.hidden_steps = 0;
 504       }else{
 505       }
 506
 507       if (!_converged && firstloop && Nest) { untested();
 508         trace1("didnt converge in first", Nest);
 509         return;
 510       }
 511       // ::status.accept.start();
 512       // _sim->set_limit();
 513       // CARD_LIST::card_list.tr_accept();
 514       // ::status.accept.stop();
 515       // _sim->_has_op = _sim->_mode;
 516       // _sim->keep_voltages();
 517       // outdata(*_sweepval[Nest]); ??
 518       itl = OPT::DCXFER;
 519
 520     }
 521
 522     if(firstloop){
 523       if(step_cause() != scUSER){
 524         trace1("firststep nouser", Nest);
 525         return;
 526       }
 527     } else {
 528       // UGLY. next may have changed _reverse[Nest]
 529       step = add;
 530       back = sub;
 531       if (!_linswp[Nest]) {
 532         step=mul;
 533         back=div;
 534       }
 535       if (_reverse[Nest]) {
 536         std::swap(step,back);
 537       }
 538       // /UGLY
 539
 540     }
 541     if ((firstloop || _converged) && step_cause() == scUSER) {
 542       _val_by_user_request[Nest] = step(_val_by_user_request[Nest], _step[Nest]);
 543       trace2("ordered next step loop", Nest, _val_by_user_request[Nest]);
 544     }
 545     firstloop = false;
 546   } while (next(Nest));
 547
 548   // _sim->pop_voltages();
 549 }
 550 /*--------------------------------------------------------------------------*/
 551 void DCOP::first(int Nest)
 552 {
 553   assert(Nest >= 0);
 554   assert(Nest < DCNEST);
 555   assert(_start);
 556   assert(_sweepval);
 557   assert(_sweepval[Nest]);
 558
 559   if (ELEMENT* c = dynamic_cast<ELEMENT*>(_zap[Nest])) {
 560     c->set_constant(false);
 561     // because of extra precalc_last could set constant to true
 562     // will be obsolete, once pointer hack is fixed
 563   }else{
 564     // not needed if not sweeping an element
 565   }
 566
 567   *_sweepval[Nest] = _start[Nest];
 568   if(_converged){
 569     trace1("BUG", Nest);
 570 //    set_step_cause(scUSER);
 571   }else{
 572   }
 573   assert(step_cause());
 574
 575   _val_by_user_request[Nest] = _start[Nest];
 576   _sweepdamp[Nest] = 1;
 577   if (ELEMENT* c = dynamic_cast<ELEMENT*>(_zap[Nest])) {
 578     // because of extra precalc_last
 579     // obsolete, once pointer hack is fixed
 580     c->set_constant(false);
 581     trace1("zapq", _zap[Nest]->long_label());
 582 //    CARD_LIST::card_list.q_hack(_zap[Nest]);
 583   }else{
 584   }
 585   _reverse[Nest] = false;
 586   if (_reverse_in[Nest]) {
 587     _converged = true; //?
 588     double (*step)(double a, double b) = add;
 589     double (*back)(double a, double b) = sub;
 590     if (!_linswp[Nest]) {
 591       step=mul;
 592       back=div;
 593     }
 594     while (next(Nest)) {
 595       _val_by_user_request[Nest] = step(*(_sweepval[Nest]), _step[Nest]);
 596     }
 597     _val_by_user_request[Nest] = back(*(_sweepval[Nest]), _step[Nest]);
 598     _reverse[Nest] = true;
 599     next(Nest);
 600     _converged = false;
 601   }else{
 602   }
 603   _sim->_phase = p_INIT_DC;
 604 }
 605 /*--------------------------------------------------------------------------*/
 606 bool DCOP::next(int Nest)
 607 { itested();
 608   double sweepval = NOT_VALID;
 609   trace2("next", Nest, *(_sweepval[Nest]));
 610   bool ok = false;
 611   double nothing = 0;
 612   double (*step)(double a, double b) = add;
 613   double (*back)(double a, double b) = sub;
 614   bool (*further)(double a, double b) = ge;
 615   double (*scale)(double a, double b) = mul;
 616   if (!_linswp[Nest]) {
 617     step = mul;
 618     back = div;
 619     scale = pow;
 620     nothing = 1;
 621     scale = pow;
 622   }
 623   double fudge = scale(_step[Nest], 1e-6);
 624   if (_reverse[Nest]) {
 625     trace2("next, reverse", *_sweepval[Nest], _val_by_user_request[Nest]);
 626     std::swap(step,back);
 627     if (_step[Nest] < 0) {
 628       further = ge;
 629     }else{ itested();
 630       further = le;
 631     }
 632     fudge = scale(_step[Nest], -1e-6);
 633   } else if (_step[Nest] < 0) {
 634     further = le;
 635   }
 636   if (_step[Nest] == nothing) {
 637     ok = false;
 638     set_step_cause(scZERO);
 639   }else if (!_converged && _ever_converged) {
 640     if (_sweepdamp[Nest]<OPT::dtmin) { untested();
 641       throw Exception("step too small (does not converge)");
 642     }else{
 643     }
 644     _sweepdamp[Nest] /= 2.;
 645     trace2("reducing step by", _sweepdamp[Nest], Nest);
 646     sweepval = back(*_sweepval[Nest], scale(_step[Nest],_sweepdamp[Nest]));
 647     trace2("next at", sweepval, _val_by_user_request[Nest]);
 648     ok = true;
 649     set_step_cause(scREJECT);
 650   }else{
 651     if (_sweepdamp[Nest] != 1) {
 652       trace3("recovered from", _sweepdamp[Nest], Nest, sweepval);
 653       set_step_cause(scTE);
 654     }else{
 655     }
 656     _sweepdamp[Nest] *= 1.4;
 657     _sweepdamp[Nest] = std::min(_sweepdamp[Nest],1.);
 658     sweepval = step(*_sweepval[Nest], scale(_step[Nest],_sweepdamp[Nest]));
 659     fixzero(&sweepval, _step[Nest]);
 660     ok = in_order(back(_start[Nest],fudge), sweepval, step(_stop[Nest],fudge));
 661     if(ok){
 662     }else if (!_reverse[Nest] && !ok && _loop[Nest]) {
 663       _reverse[Nest] = true;
 664       if (_step[Nest] < 0) {
 665         further = le;
 666       }else{
 667         further = ge;
 668       }
 669       fudge = scale(_step[Nest], -.1);
 670       sweepval = back(*_sweepval[Nest], _step[Nest]);
 671       std::swap(step,back);
 672       ok = in_order(back(_start[Nest],fudge), sweepval, step(_stop[Nest],fudge));
 673       assert(ok);
 674       trace2("BUG?", sweepval, *_sweepval[Nest]);
 675       _val_by_user_request[Nest] = sweepval; // BUG: here?
 676     }else if(_reverse_in[Nest]){
 677       // hmm maybe _reverse_in && !_reverse?
 678       trace2("reverse end?", sweepval, *_sweepval[Nest]);
 679       *_sweepval[Nest] = sweepval;
 680     }else{
 681     }
 682   }
 683
 684   double v = _val_by_user_request[Nest];
 685   if(!ok){
 686   }else if (further(step(sweepval, scale(_step[Nest],1e-6) ), v)) {
 687     trace5("userstep at", v, sweepval, ok, _reverse[Nest], *_sweepval[Nest]);
 688     set_step_cause(scUSER); // here?!
 689     sweepval = v;
 690   }else{
 691   }
 692
 693   _sim->_phase = p_DC_SWEEP;
 694 //  *(_sweepval[Nest]) = sweepval; // ouch.
 695   if (ok) {
 696     assert(sweepval != NOT_VALID);
 697     *(_sweepval[Nest]) = sweepval;
 698     return true;
 699   }else{
 700     trace3("not ok at", v, sweepval, *(_sweepval[Nest]));
 701     //assert(sweepval == NOT_VALID);
 702     return false;
 703   }
 704 }
 705 /*--------------------------------------------------------------------------*/
 706 static DC p2;
 707 static OP p4;
 708 static DISPATCHER<CMD>::INSTALL d2(&command_dispatcher, "dc", &p2);
 709 static DISPATCHER<CMD>::INSTALL d4(&command_dispatcher, "op", &p4);
 710 }
 711 /*--------------------------------------------------------------------------*/
 712 /*--------------------------------------------------------------------------*/
 713 // vim:ts=8:sw=2:noet: