?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
fix 2 compiler warnings in modelgen
[gnucap-felix.git] / modules / s_ev.cc
blob5d762bb44a4c67274316b881066c74fcdb8762f5
1 /* Copyright (C) 2014 Felix Salfelder
2 * Author: Felix Salfelder
3 *
4 * This file is part of "Gnucap", the Gnu Circuit Analysis Package
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3, or (at your option)
9 * any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
19 * 02110-1301, USA.
20 *------------------------------------------------------------------
21 * eigenvalue stuff
22 */
23
24 #define ADD_VERSION
25
26 #include "config.h"
27
28 #ifdef HAVE_CLAPACK_H
29 // #include "u_atlas.h"
30 #endif
31
32
33 #ifndef HAVE_CLAPACK_H
34 #warning "untested"
35 #endif
36
37 #include "u_status.h"
38 #include "u_ssp.h"
39 #include "s_ev.h"
40 #include "u_gsl.h"
41 #include "u_prblst.h"
42 #include "u_cardst.h"
43 #include "e_elemnt.h"
44 #include "e_storag.h"
45 #include "s_ev.h"
46 #include "io_matrix.h"
47 #include "m_matrix_extra.h"
48 #include "e_aux.h"
49 using namespace std;
50 /*--------------------------------------------------------------------------*/
51 /*--------------------------------------------------------------------------*/
52 /*--------------------------------------------------------------------------*/
53 /*--------------------------------------------------------------------------*/
54 namespace { //
55 class EV : public EV_BASE { //
56 public:
57 explicit EV(): EV_BASE() {}
58 ~EV() {}
59 void do_it(CS&, CARD_LIST*);
60 private:
61 void setup(CS&);
62 explicit EV(const EV&): EV_BASE() {unreachable(); incomplete();}
63
64 protected:
65 void sweep_recursive(int);
66
67 void align_inf_eigenvectors(gsl_matrix_complex_view& EV_inf, size_t inputnumber=0);
68 unsigned sens_inf_vectors(gsl_matrix_complex_view& EV);
69
70 private:
71 };
72 /*--------------------------------------------------------------------------*/
73 /*--------------------------------------------------------------------------*/
74 // EV_DISC?
75 void EV::do_it(CS& Cmd, CARD_LIST* Scope)
76 {
77 trace0("EV::do_it(CS&, CARD_LIST*)");
78 _scope = Scope;
79 _sim->_time0 = 0.;
80 _sim->set_command_ddc();
81 _sim->_phase = p_INIT_DC;
82 ::status.ddc.reset().start();
83 _sim->_temp_c = temp_c_in;
84 _do_tran_step = 0;
85 _dump_matrix = 0;
86 command_base(Cmd);
87 ::status.ddc.stop();
88 }
89 /*--------------------------------------------------------------------------*/
90 /*--------------------------------------------------------------------------*/
91 void EV::setup(CS& Cmd)
92 {
93 _cont = false;
94 _trace = tNONE;
95 _out = IO::mstdout;
96 _out.reset(); //BUG// don't know why this is needed */
97 bool ploton = IO::plotset && plotlist().size() > 0;
98
99 _n_inputs = 0;
100
101 if (Cmd.more()) {
102 for (_n_sweeps = 0; Cmd.more() && _n_sweeps < DCNEST; ++_n_sweeps) {
103 if (Cmd.umatch("to")) {
104 break;
105 }
106 CARD_LIST::fat_iterator ci = findbranch(Cmd, &CARD_LIST::card_list);
107 if (!ci.is_end()) {
108 if (ELEMENT* c = dynamic_cast<ELEMENT*>(*ci)) {
109 _zap[_n_sweeps] = c;
110 trace1("EV::setup", _zap[_n_sweeps]->long_label());
111 } else { untested();
112 throw Exception("dc/op: can't sweep " + (**ci).long_label() + '\n');
113 }
114 _inputno[_n_sweeps] = _n_inputs;
115 ++_n_inputs;
116 }else if (Cmd.is_float()) { // sweep the generator
117 _zap[_n_sweeps] = NULL;
118 }else if (Cmd.is_alpha()) {
119 string pname;
120 unsigned here = Cmd.cursor();
121 Cmd >> pname;
122 try {
123 _param[_n_sweeps] = &(_scope->params()->find(pname));
124 } catch(Exception) {
125 Cmd.reset(here);
126 // throw Exception("ddc: can't sweep " + pname + '\n');
127 }
128 _zap[_n_sweeps] = NULL;
129 }else{ untested();
130 // leave as it was .. repeat Cmd with no args
131 }
132 if (Cmd.match1("'\"({") || Cmd.is_float()) {
133 _start[_n_sweeps] = "NA";
134 _stop[_n_sweeps] = "NA";
135 Cmd >> _start[_n_sweeps] >> _stop[_n_sweeps];
136 _step[_n_sweeps] = 0.;
137 }else{
138 // leave it as it was .. repeat Cmd with no args
139 }
140 _sim->_genout = 0.;
141 temp_c_in = OPT::temp_c;
142 _sim->_temp_c = temp_c_in;
143 options(Cmd,_n_sweeps);
144
145 }
146 unsigned here = Cmd.cursor();
147 string output;
148 bool newprobes = true;
149 do{ itested();
150 // FIXME: can we use the probe parser?!
151 if( Cmd >> "v" ){
152 output_t t;
153 t.brh[1] = 0;
154 trace1("SENS::setup, have v", Cmd.tail());
155
156 int paren = Cmd.skip1b('(');
157 Cmd >> output;
158 t.label = "v("+output;
159
160 CKT_NODE* node = dynamic_cast<CKT_NODE*>((*_scope).node(output));
161 if(node)
162 t.brh[0] = node;
163 else{
164 continue;
165 }
166
167 trace2("SENS::setup, have", output, Cmd.tail());
168 if(!Cmd.match1(')')){
169 output = Cmd.ctos(")");
170 trace1("SENS:skip1b:setup, have2", output);
171 t.label += "," + output;
172 node = dynamic_cast<CKT_NODE*>((*_scope).node(output));
173 if(node)
174 t.brh[1] = node;
175 else{
176 Cmd.warn(bWARNING, "probelist: what's this?");
177 }
178
179 }else{
180 trace0("SENS::setup no comma");
181 }
182
183 paren -= Cmd.skip1b(')');
184 if (paren != 0) {untested();
185 Cmd.warn(bWARNING, "need )");
186 }else if (output.empty()) {untested();
187 Cmd.warn(bWARNING, "probelist: what's this?");
188 }else{
189 }
190
191 t.label+=")";
192 if (newprobes){
193 _output.clear();
194 newprobes = false;
195 }
196 _output.push_back(t);
197 }
198 //try{
199 // _output.add_list(Cmd);
200 //}
201 //catch(Exception_Cant_Find)
202 //{}
203 ONE_OF
204 || Get(Cmd, "dm", &_dump_matrix)
205 || _out.outset(Cmd);
206 ;
207 }while (Cmd.more() && !Cmd.stuck(&here));
208 }else{ untested();
209 }
210 Cmd.check(bWARNING, "what's this?");
211
212 // hat peter auskommentiert
213 //_sim->_uic = _sim->_more_uic = true;
214
215 IO::plotout = (ploton) ? IO::mstdout : OMSTREAM();
216 initio(_out);
217
218 assert(_n_sweeps > 0);
219 for (int ii = 0; ii < _n_sweeps; ++ii) {
220 trace1("DDC_setup", ii);
221 _start[ii].e_val(0., _scope);
222 fix_args(ii);
223
224 if (_zap[ii]) {
225 trace2("zap inc_probes" + _zap[ii]->long_label(), ii, _zap[ii]->probes() );
226 _stash[ii] = _zap[ii]; // stash the std value
227 _zap[ii]->inc_probes(); // we need to keep track of it
228
229 STORAGE* s = dynamic_cast<STORAGE*>(_zap[ii]);
230 if (s) {
231 trace2("EV::setup ", _zap[ii]->long_label(), _zap[ii]->is_constant());
232 _zap[ii]->set_constant(false); // so it will be updated
233 _pushel[ii] = _zap[ii]; // point to value to patch
234 _uic_caplist.push_back(s);
235 _sweepval[ii] = s->set__ic();
236 }else{
237 trace1("EV::setup, no STORAGE", ii);
238 _zap[ii]->set_value(_zap[ii]->value(),0); // zap out extensions
239 _zap[ii]->set_constant(false); // so it will be updated
240 _pushel[ii] = _zap[ii]; // element to patch
241 _sweepval[ii] = _zap[ii]->set__value();
242 }
243 } else if (_param[ii]) {
244 _sweepval[ii] = _param[ii]->pointer_hack();
245 } else {
246 trace1("EV::setup does not exist", ii);
247 //_sweepval[ii] = 0;
248 _sweepval[ii] = &_sim->_genout;
249 _pushel[ii] = NULL; //&_sim->set_gen; // point to value to patch
250 // throw(Exception("something went wrong\n"));
251 }
252 }
253 _sim->_freq = 0;
254
255 }
256 /*--------------------------------------------------------------------------*/
257 static EV p2;
258 static DISPATCHER<CMD>::INSTALL d2(&command_dispatcher, "evs", &p2);
259 /*--------------------------------------------------------------------------*/
260 }// namespace
261 static void register_status()
262 {
263 static bool done;
264 if(!done) { itested();
265 new DISPATCHER<CKT_BASE>::INSTALL(&status_dispatcher, "ssp", &p2);
266 done = true;
267 }
268 }
269 /*--------------------------------------------------------------------------*/
270 static size_t n;
271 void EV_BASE::sweep()
272 { itested();
273 register_status();
274
275 trace0("EV_BASE::sweep()");
276
277 _sim->set_command_tran();
278 head(_start[0], _stop[0], "chart spl"); // here? hmmm
279 // _sim->_bypass_ok = false;
280 _sim->set_inc_mode_bad();
281 if (_cont) {untested();
282 _sim->restore_voltages();
283 CARD_LIST::card_list.tr_restore();
284 }else{
285 _sim->clear_limit();
286 CARD_LIST::card_list.tr_begin();
287 }
288
289 // allocate()?
290 n = _sim->_total_nodes;
291 // FIXME: allocate globally, where necessary.
292 G_ = gsl_matrix_alloc(n,n);
293 C_ = gsl_matrix_alloc(n,n);
294 Q_ = gsl_matrix_alloc(n,n);
295 Z_ = gsl_matrix_alloc(n,n);
296
297 sweep_recursive(_n_sweeps);
298 }
299 /*--------------------------------------------------------------------------*/
300 /*--------------------------------------------------------------------------*/
301 double EV_BASE::temp_c_in = 0.;
302 /*--------------------------------------------------------------------------*/
303 EV_BASE::EV_BASE()
304 :SIM(),
305 _n_sweeps(1),
306 _cont(false),
307 _trace(tNONE)
308 {
309
310 for (int ii = 0; ii < DCNEST; ++ii) {
311 _loop[ii] = false;
312 _reverse_in[ii] = false;
313 _reverse[ii] = false;
314 _step[ii]=0.;
315 _linswp[ii]=true;
316 _sweepval[ii]=&_sim->_genout;
317 _zap[ii]=NULL;
318 _stepmode[ii] = ONE_PT;
319 _param[ii]=NULL;
320 }
321
322 temp_c_in=OPT::temp_c;
323 _out=IO::mstdout;
324 }
325 /*--------------------------------------------------------------------------*/
326 void EV_BASE::finish(void)
327 {
328 trace0("EV_BASE::finish");
329 // _sim->_phase = p_;
330 for (int ii = 0; ii < _n_sweeps; ++ii) {
331 if (_zap[ii]) { // component
332 _stash[ii].restore();
333 _zap[ii]->dec_probes();
334 trace2("EV_BASE::finish dec_probes done", _zap[ii]->long_label(), _zap[ii]->probes());
335 _zap[ii]->precalc_first();
336 _zap[ii]->precalc_last();
337 }else{
338 }
339 }
340 _uic_caplist.clear();
341 }
342 /*--------------------------------------------------------------------------*/
343 void EV_BASE::fix_args(int Nest)
344 {
345 trace1("EV_BASE::fix_args(int Nest)", Nest);
346
347 _stop[Nest].e_val(_start[Nest], _scope);
348 _step_in[Nest].e_val(0., _scope);
349 _step[Nest] = _step_in[Nest];
350
351 switch (_stepmode[Nest]) { untested();
352 case ONE_PT:
353 case LIN_STEP:
354 _linswp[Nest] = true;
355 break;
356 case LIN_PTS:untested();
357 if (_step[Nest] <= 2.) {untested();
358 _step[Nest] = 2.;
359 }else{untested();
360 }
361 _linswp[Nest] = true;
362 break;
363 case TIMES:untested();
364 if (_step[Nest] == 0. && _start[Nest] != 0.) {untested();
365 _step[Nest] = _stop[Nest] / _start[Nest];
366 }else{untested();
367 }
368 _linswp[Nest] = false;
369 break;
370 case OCTAVE:
371 if (_step[Nest] == 0.) {untested();
372 _step[Nest] = 1.;
373 }else{ untested();
374 }
375 _step[Nest] = pow(2.00000001, 1./_step[Nest]);
376 _linswp[Nest] = false;
377 break;
378 case DECADE:
379 if (_step[Nest] == 0.) { untested();
380 _step[Nest] = 1.;
381 }else{ untested();
382 }
383 _step[Nest] = pow(10., 1./_step[Nest]);
384 _linswp[Nest] = false;
385 break;
386 };
387
388 if (_step[Nest] == 0.) { // prohibit log sweep from 0
389 _step[Nest] = _stop[Nest] - _start[Nest];
390 _linswp[Nest] = true;
391 }else{
392 }
393 }
394 /*--------------------------------------------------------------------------*/
395 // process tail vector.
396 // find maximum coeff in v0prime_tail at piv, swap to 0, and
397 // exchange, rescale first basis vector to get z0_tail = (1,0,0 ...)
398 // maybe unnecessary?
399 int EV_BASE::process_tail(gsl_matrix_complex_view& relevant_EV)
400 {
401 int err = 0;
402 size_t n = _sim->_total_nodes;
403 unsigned tailstart = unsigned(relevant_EV.matrix.size2);
404 unsigned tailsize = unsigned(n-tailstart);
405 trace2("preparing tail", tailstart, tailsize);
406 if(tailstart < n) {
407 gsl_matrix_complex_view EV_tail = gsl_matrix_complex_submatrix (EV_, 0, tailstart, n, tailsize);
408 gsl_vector_complex_view v0prime_tail = gsl_vector_complex_subvector(z0_, tailstart, tailsize);
409 size_t piv = gsl_blas_izamax(&v0prime_tail.vector);
410 if(piv != 0){
411 gsl_vector_complex_swap_elements(&v0prime_tail.vector,0,piv);
412 gsl_matrix_complex_swap_columns(&EV_tail.matrix,0,piv);
413 }
414 gsl_complex pivot = gsl_vector_complex_get(&v0prime_tail.vector,0);
415 if (0.!=GSL_REAL(pivot) || 0.!=GSL_IMAG(pivot)){
416 gsl_vector_complex_view first = gsl_matrix_complex_column(EV_, tailstart);
417
418 for(unsigned i=tailstart+1; i<n; ++i) {
419 gsl_complex s = gsl_vector_complex_get(z0_,i);
420 gsl_vector_complex_const_view b = gsl_matrix_complex_const_column(EV_, i);
421 gsl_blas_zaxpy(gsl_complex_div(s,pivot), &b.vector, &first.vector);
422
423 gsl_vector_complex_set(z0_,i,zero);
424 }
425
426 double norm = 1;
427 gsl_complex* p = gsl_vector_complex_ptr(&v0prime_tail.vector,0);
428 if(0){
429 norm = 1/gsl_blas_dznrm2(&first.vector);
430 }else if(1) {
431 norm = gsl_complex_abs(*p);
432 }
433 trace1("normalize tail[0]", norm);
434 err+= gsl_vector_complex_scale(&first.vector, gsl_complex_mul(one,gsl_complex_rect(norm,0.))); // inefficient.
435 GSL_REAL(*p)/= norm;
436 GSL_IMAG(*p)/= norm;
437 trace1("tail", first);
438 } else { itested();
439 trace2("no pivot. possible zero op", tailstart, tailsize);
440 }
441 }else{
442 trace2("not preparing tail. possible zero op", tailstart, tailsize);
443 }
444 return err;
445 }
446 /*--------------------------------------------------------------------------*/
447 // match _sweepval and dcsens coeffs.
448 inline void EV_BASE::compute_canonical_op(gsl_vector_complex_view EV_op, gsl_matrix_complex_const_view dcsens)
449 {
450 trace1("compute_canonical_op", *z0_);
451 for(unsigned i=0; i<_n_inputs; ++i) {
452 gsl_complex* k_i = gsl_vector_complex_ptr(z0_, _n_states + i);
453 gsl_complex delta = gsl_complex_sub(*k_i, gsl_complex_rect(*_sweepval[_inputno[i]],0.));
454 trace3("inputs vs inputs", *k_i, *_sweepval[_inputno[i]], delta);
455
456 *k_i = gsl_complex_rect(*_sweepval[_inputno[i]], 0.);
457
458 // op -= k_i * delta;
459 gsl_vector_complex_const_view s = gsl_matrix_complex_const_column(&dcsens.matrix,i);
460 gsl_blas_zaxpy( delta, &s.vector, &EV_op.vector);
461 }
462
463 // for(unsigned i=0; i<_n_inputs; ++i) {
464 // gsl_complex* op_k = gsl_vector_complex_ptr(z0_, _n_states + i);
465 // trace2("inputs vs inputs", *op_k, *_sweepval[_inputno[i]]);
466 // }
467 }
468 /*--------------------------------------------------------------------------*/
469 // foreach input
470 // compute node voltage DC sensitivities, stash in columns of dc_sens
471 // project to infspace
472 // fill inf vectors, preserving linear independence.
473 //
474 // append op to ensure v0 is in img. FIXME
475 //
476 inline unsigned EV_BASE::sens_inf_vectors(gsl_matrix_complex_view EV_inf, gsl_vector_view grid, gsl_matrix_complex* nodc_sens)
477 {
478 size_t n = _sim->_total_nodes;
479 assert(n==EV_inf.matrix.size1); // number of rows
480 size_t m = EV_inf.matrix.size2;
481 unsigned seek = 0;
482 trace4("sens_inf_vectors", EV_inf.matrix.size1, _n_sweeps, n, EV_inf.matrix.size2);
483 trace1("sens_inf_vectors", grid);
484
485 for(unsigned i=0; i<unsigned(_n_sweeps); ++i){ itested();
486 if(!_zap[i]){ untested();
487 // sweeping over parameter?
488 continue;
489 }
490
491 clear_sens();
492 _zap[i]->sens_load("dummy"); // sens input
493 complex_array_to_gsl_vector_complex(*sens_,_sim->_sens+1,n);
494 trace3("sensstamp", _zap[i]->long_label(), *sens_, _step[i]);
495
496 ::status.back.start();
497 _sim->_lu.fbsub(_sim->_sens);
498 ::status.back.stop();
499
500 trace1("input sens\n", _sim->_aa);
501
502 complex_array_to_gsl_vector_complex(*sens_,_sim->_sens+1,n);
503 gsl_vector_complex_view nodcsens = gsl_matrix_complex_column(nodc_sens, seek);
504
505 trace2("DC input sens", *sens_, *betabar_);
506
507 // gsl_vector_complex_memcpy(&dcsens.vector, sens_);
508
509 bool inf_sens = true;
510 if(inf_sens){ incomplete();
511 // hmm, no-DC sensitivity...
512 gsl_linalg_complex_LU_svx(EVLU_, LU_perm, sens_); // too many, only need left ones...
513 trace3("sensprime\n", *sens_, m,n);
514 gsl_vector_complex_view sens_right = gsl_vector_complex_subvector(sens_, n-m, m);
515
516 gsl_blas_zgemv( CblasNoTrans, one, &EV_inf.matrix, &sens_right.vector, zero, ztmp_ );
517
518 if(1){ // sanitycheck
519 mul(CblasTrans, *C_, *ztmp_, *sens_);
520 trace1("zero", *sens_);
521 }
522 trace1("sens projection", *ztmp_);
523 gsl_vector_complex_memcpy(&nodcsens.vector, ztmp_);
524 }else{
525 }
526
527 complex_array_to_gsl_vector_complex(*sens_,_sim->_sens+1,n);
528
529 if(1){ //emplace...
530 gsl_vector_complex_view seek_ = gsl_matrix_complex_column(&EV_inf.matrix, seek);
531 gsl_vector_complex_memcpy(&seek_.vector, sens_);
532 }else{ // rewrite basis... inefficient
533 // BUG. refresh dc_sens!!
534 gsl_vector_complex_const_view seek_ = gsl_matrix_complex_const_column(&EV_inf.matrix, seek);
535 gsl_vector_complex_sub(sens_, &seek_.vector);
536 if (0) for(size_t j=seek; j<m; ++j){
537 trace1("replacing", j);
538 gsl_vector_complex_view seek_ = gsl_matrix_complex_column(&EV_inf.matrix, j);
539 gsl_vector_complex_add(&seek_.vector, sens_);
540
541 // normalize
542 gsl_complex len;
543 GSL_IMAG(len) = 0.;
544 GSL_REAL(len) = 1./norm(seek_.vector);
545 gsl_vector_complex_scale(&seek_.vector, len);
546 }
547 // BUG. should use d z0dot/din as grid size!!
548 // FIXME: use chart!!1
549 }
550 {
551 double* g = gsl_vector_ptr(&grid.vector,seek);
552 if (fabs(_step[i]) < fabs(*g)){ itested();
553 *g = fabs(_step[i]);
554 // trace2("sens_inf_vectors", _step[i], c.grid(i)); // FIXME
555 }
556 }
557 ++seek;
558 }
559 trace1("done sens_inf_vectors\n", EV_inf.matrix);
560 return seek;
561 }
562 /*--------------------------------------------------------------------------*/
563 gsl_complex operator-(const gsl_complex&x){
564 return gsl_complex_negative(x);
565 }
566 /*--------------------------------------------------------------------------*/
567 SSP_STATE& EV_BASE::ev_solver(SSP_CHART*& chart)
568 {
569
570 // FIXME: fetch more information from chart.
571 // - number of inputs
572 // - number of states.
573 // - grid stuff
574
575 int err;
576 size_t n = _sim->_total_nodes;
577
578 trace2("========= ev_solver ========\n", _sim->_acx, SSP_TREE::num_charts());
579 assert(SSP_TREE::num_charts() < (1<<20));
580 // gsl_eigen_gen_workspace* w = gsl_eigen_gen_alloc(n);
581 gsl_eigen_genv_workspace* w = gsl_eigen_genv_alloc(n);
582 gsl_vector *tmp_;
583 gsl_vector_complex *v0_;
584 gsl_vector *grid_;
585 gsl_vector_view v0__;
586 gsl_vector_complex *v0back_;
587
588 z0_ = gsl_vector_complex_alloc(n);
589 z0dot_ = gsl_vector_complex_alloc(n);
590 v0back_ = gsl_vector_complex_alloc(n);
591
592 BSMATRIX<double> G = _sim->_acx.real();
593 BSMATRIX<double> C = _sim->_acx.imag();
594
595
596 // P_ = gsl_matrix_alloc(n,n);
597 EV_ = gsl_matrix_complex_alloc(n,n);
598 EVLU_ = gsl_matrix_complex_alloc(n,n);
599 LU_perm = gsl_permutation_alloc (n);
600
601 alpha_ = gsl_vector_complex_alloc(n);
602 beta_ = gsl_vector_alloc(n);
603 betabar_ = gsl_vector_calloc(n);
604 v0_ = gsl_vector_complex_alloc(n);
605 sens_ = gsl_vector_complex_alloc(n);
606 grid_ = gsl_vector_alloc(n);
607 tmp_ = gsl_vector_alloc(n);
608 ztmp_ = gsl_vector_complex_alloc(n);
609 lambda_ = gsl_vector_complex_alloc(n); // fixme: only need maximum _n_states
610 // _n_states may depend on operating point...
611 G >> *G_;
612 C >> *C_;
613 v0__ = gsl_vector_view_array(_sim->_v0+1,n);
614 trace1("v0", v0__.vector);
615
616 trace1("A = G =\n", *G_);
617 trace1("B = C =\n", *C_);
618
619 // want \lambda C x = G x
620 // or 1/\mu C x = G x
621 // solve beta G = alpha C
622 // lambda_i = alpha/beta
623 // mu_i = beta/alpha
624
625 err+= gsl_eigen_genv_QZ (G_, C_, alpha_, beta_, EV_, Q_, Z_, w);
626 trace1("raw genv\n", *EV_);
627 trace1("schur\n", *G_);
628 trace1("schur\n", *C_);
629
630 { // check.
631 err+= gsl_vector_complex_memcpy(lambda_,alpha_);
632 *lambda_ /= *beta_;
633 }
634 trace2("solved alpha G = beta C\n", *alpha_, *beta_);
635 trace1("lambda\n", *lambda_);
636
637 C >> *C_;
638
639 //gsl_eigen_gen_QZ (G_, C_, alpha_, beta_, Q_, Z_, w);
640
641 if(0){ // hmm check
642 gsl_matrix *test_;
643 test_ = gsl_matrix_alloc(n,n);
644 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1., Q_, Z_, 0., test_);
645 trace1("QZ^T?", *test_);
646 }
647 // { // hmm check
648 // gsl_matrix_complex *test_;
649 // test_ = gsl_matrix_complex_alloc(n,n);
650 // gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, EV_, EV_, zero, test_);
651 // trace1("identity?", *test_);
652 // }
653 //
654 to_norm(*grid_, *alpha_);
655 gsl_vector_div(grid_, beta_);
656 trace1("lambda tmp grid", *grid_);
657 trace1("...", *lambda_);
658
659 gsl_permutation *ev_perm = gsl_permutation_alloc (n);
660
661 // order_by_lambda(ev_perm, grid_ find_dim)
662 { // order eigenvalues by size (starting with small)
663 //
664 // weed out huge lambdas.
665 // replace 1st nan eigenvector by v0 ... hmm
666 // apply order to EV_
667 size_t fin_dim = (size_t)n;
668
669 err+= sort_vector_and_index (ev_perm, grid_, fin_dim);
670 trace2("evalue perm", *ev_perm, fin_dim);
671
672 for(unsigned i=0; i<fin_dim; ++i){ itested();
673 gsl_vector_set(betabar_,i,1);
674 }
675
676 trace1("beta", *betabar_);
677
678 for(unsigned i=0; i<n; ++i){
679 gsl_vector_complex_view v;
680 gsl_vector_view vr;
681 v = gsl_matrix_complex_row (EV_, i);
682 err+= gsl_permute_vector_complex (ev_perm, &v.vector);
683 vr = gsl_matrix_row (Q_, i);
684 err+= gsl_permute_vector (ev_perm, &vr.vector);
685 vr = gsl_matrix_row (Z_, i);
686 err+= gsl_permute_vector (ev_perm, &vr.vector);
687 }
688 trace1("ordered\n", *EV_);
689 if (0) {
690 gsl_matrix_complex *test_;
691 test_ = gsl_matrix_complex_alloc(n,n);
692 gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, EV_, EV_, zero, test_);
693 trace1("ordered, identity?", *test_);
694 }
695 // first n rows of EV_ (the non inf-EVs)
696
697
698 real_array_to_gsl_vector(*z0_,_sim->_v0+1,n);
699 trace1("v0", *z0_);
700 // {
701 // trace1("preimg1", *EV_);
702 // err+= gsl_blas_zgemv( CblasConjTrans, one, EV_, z0_, zero, ztmp_ );
703 // trace1("preimg1", *ztmp_);
704 // }
705
706 gsl_matrix_complex_memcpy (EVLU_, EV_);
707
708 int signum;
709 gsl_linalg_complex_LU_decomp(EVLU_, LU_perm, &signum);
710
711 // probably not necessary.
712 // compute EV^-1 v0 below, with modified EV
713 real_array_to_gsl_vector(*z0_,_sim->_v0+1,n);
714 gsl_linalg_complex_LU_svx(EVLU_, LU_perm, z0_);
715 trace2("preimg2 ", *z0_, signum);
716
717 // reverse test
718 err+= gsl_blas_zgemv( CblasNoTrans, one, EV_, z0_, zero, ztmp_ );
719 trace1("v0 again?", *ztmp_);
720 _n_states = (unsigned)fin_dim;
721 }
722
723 if(0){ //trying to figure out coeffs...
724 err+= gsl_blas_dgemv( CblasTrans, 1, Z_, &v0__.vector, 0, tmp_ );
725 err+= gsl_permute_vector (ev_perm, tmp_);
726 trace1("preimg3 ", *tmp_);
727 set_vector(*z0_,*tmp_);
728 err+= gsl_blas_zgemv( CblasNoTrans, one, EV_, z0_, zero, ztmp_ );
729 trace1("v0 again?", *ztmp_);
730 }
731
732 ::status.lud.start();
733 COMPLEX z = COMPLEX(OPT::gmin,0.);
734 _sim->_acx.dezero(z);
735 _sim->_acx.lu_decomp();
736 ::status.lud.stop();
737
738 trace2("hmm", n, _n_inputs);
739 real_array_to_gsl_vector(*v0_,_sim->_v0+1,n);
740
741 // fin_dim: dimension of noninf space.
742 // input_dim: dimension of constant Z input space
743 // == number of inputs (?)
744
745 unsigned input_dim;
746 gsl_matrix_complex* nodc_sens_ = gsl_matrix_complex_alloc(n, _n_inputs);
747 { itested();
748 trace1("pre input sens", *grid_);
749 // compute input sensitivities.
750 // add to basis ?
751 gsl_matrix_complex_view EV_inf = gsl_matrix_complex_submatrix (EV_, 0, _n_states, n, n-_n_states);
752 gsl_vector_view g = gsl_vector_subvector (grid_, _n_states, n-_n_states);
753 trace3("infspace\n", EV_inf.matrix, EV_inf.matrix.size1, EV_inf.matrix.size2);
754 input_dim = sens_inf_vectors(EV_inf, g, nodc_sens_);
755 trace1("modified EV", input_dim);
756 trace1("done sens_inf_vectors", *grid_);
757 trace1("done sens_inf_vectors", *nodc_sens_);
758 }
759 assert(input_dim==unsigned(_n_inputs));
760
761 trace1("tmp grid", *grid_);
762
763 { // recompute LU decomp. inefficient!
764 // why?
765 int signum;
766 gsl_matrix_complex_memcpy (EVLU_, EV_);
767 gsl_linalg_complex_LU_decomp(EVLU_, LU_perm, &signum);
768 }
769
770 { // v0' = EV^-1 v0
771 real_array_to_gsl_vector(*z0_,_sim->_v0+1,n);
772 trace1("v0", *z0_);
773 gsl_linalg_complex_LU_svx(EVLU_, LU_perm, z0_);
774 trace1("", *z0_);
775 }
776
777 {
778 gsl_matrix_complex_view relevant_EV;
779 relevant_EV = gsl_matrix_complex_submatrix (EV_, 0, 0, n, _n_states + input_dim);
780 trace1("unprocessed tail\n", *EV_);
781 trace1("unprocessed tail", *z0_);
782 process_tail(relevant_EV);
783 trace1("processed tail\n", *EV_);
784 trace1("processed tail", *z0_);
785 }
786
787 // gsl_vector_complex* z0limit_ = gsl_vector_complex_alloc(n);
788 #if 0
789 if(_n_states) { // new grid?
790 trace2("new grid", _n_states, input_dim);
791 gsl_vector_complex_memcpy(z0dot_, z0_);
792 gsl_vector_complex_view z0dot__ = gsl_vector_complex_subvector(z0dot_, 0, _n_states);
793 do_ddc(z0dot__);
794 gsl_vector_complex_memcpy(z0limit_, ztmp_);
795 }else{ // no movement, already at limit
796 gsl_vector_complex_memcpy(z0limit_, z0_);
797 }
798 #endif
799
800 trace1("precan op", *z0_);
801 trace1("precan op\n" , *EV_);
802 if(_n_states+_n_inputs < n) { itested();
803 gsl_vector_complex_view op__ = gsl_matrix_complex_column(EV_, _n_states+_n_inputs);
804 gsl_matrix_complex_const_view dcs = gsl_matrix_complex_const_submatrix (EV_, 0, _n_states, n, n -_n_states);
805 trace2("can op EV", _n_states, _n_inputs);
806 trace3("can op EV", _n_states, n, n-_n_states);
807 compute_canonical_op(op__, dcs);
808 trace1("can op", *z0_);
809 trace1("can op", op__);
810 trace1("can op\n", *EV_);
811 }
812
813 if(_n_states) { itested();
814 gsl_vector_complex_view z0dot = gsl_vector_complex_subvector(z0dot_,0,_n_states);
815 gsl_vector_complex_const_view z0 = gsl_vector_complex_const_subvector(z0_, 0, _n_states);
816 compute_zdot(z0, z0dot);
817 trace1("z0dot", z0dot);
818 } else { untested();
819 }
820
821 if(0 && _n_states){
822 gsl_vector_view g = gsl_vector_subvector (grid_, 0, _n_states);
823 trace1("scale sens", *grid_);
824 gsl_vector_scale(g, OPT::dtmin);
825 trace1("scale sens", *grid_);
826 }
827
828 {
829 gsl_matrix_complex_view relevant_EV;
830 trace1("pre output sens", *grid_);
831 trace2("pre output sens", _n_states, input_dim);
832 relevant_EV = gsl_matrix_complex_submatrix (EV_, 0, 0, n, min(unsigned(n), _n_states + input_dim + 1)); // +1??
833 output_sens_grid(relevant_EV, grid_, lambda_);
834 trace1("done output sens", *grid_);
835 }
836
837 { // recompute LU decomp for rhs. inefficient...
838 int signum;
839 gsl_matrix_complex_memcpy (EVLU_, EV_);
840 gsl_linalg_complex_LU_decomp(EVLU_, LU_perm, &signum);
841 }
842
843
844 trace1("unscaled\n", *EV_);
845 trace3("actual grid", *grid_, _n_states, input_dim);
846 to_grid(*grid_, OPT::reltol * OPT::reltol);
847 trace1("disc grid", *grid_);
848
849 if(_n_states+_n_inputs < n) {
850 // try to stabilize canonical op.
851 // doesn't work
852 gsl_vector_set(grid_,_n_states+_n_inputs,OPT::reltol);
853 }
854
855 chart = &SSP_TREE::_root.insert(*EV_, *grid_, _n_inputs);
856 assert(chart);
857
858 assert(chart->n_inputs() == _n_inputs);
859
860 if(chart->is_new()){ itested();
861 for(unsigned i=0; i<_n_inputs; ++i){
862 chart->_inputgrid[i] = gsl_vector_get(grid_,i+_n_states);
863 }
864 }else{ itested();
865 for(unsigned i=0; i<_n_inputs; ++i){
866 assert(chart->grid(i) == gsl_vector_get(grid_,i+_n_states));
867 }
868 }
869
870 assert(chart);
871 trace1("discretized\n", *EV_);
872 trace3("preminimize", _n_states, input_dim, *z0_);
873
874 trace1("v0 raw ", v0__.vector);
875 { // find z0, | v0 - EV_grid z0 | minimal.
876 size_t numcols = _n_states+input_dim;
877 bool have_op = numcols != n;
878 double rcond = 0;
879 integer rank = 0;
880 gsl_complex opcoeff;
881 integer info = 0;
882
883 if(have_op){ untested();
884 // hmm opcoeff could be one, always. but it's not (currently)
885 opcoeff = gsl_vector_complex_get(z0_, _n_states+input_dim);
886 }
887
888 real_array_to_gsl_vector(*z0_,_sim->_v0+1,n);
889
890 if(have_op){ untested();
891 gsl_vector_complex_view op__ = gsl_matrix_complex_column(EV_, _n_states+_n_inputs);
892 trace2("sub op", op__, opcoeff);
893 // gsl_vector_complex_sub(z0_, &op__.vector);
894 gsl_blas_zaxpy(-opcoeff, &op__.vector, z0_);
895 trace1("subd op", *z0_);
896 }
897
898 gsl_matrix_complex_view v0pm = gsl_matrix_complex_view_vector (z0_, 1, n);
899 trace1("zgelss", v0pm.matrix);
900
901 // use EVLU as scratch space
902 err+= gsl_matrix_complex_memcpy (EVLU_, EV_);
903 err+= gsl_matrix_complex_transpose(EVLU_);
904 gsl_matrix_complex_view zEVLU = gsl_matrix_complex_submatrix(EVLU_,0,0,numcols,n);
905 trace1("zgelss", zEVLU.matrix);
906 gsl_blas_zgelss (&zEVLU.matrix, &v0pm.matrix, tmp_, &rcond, &rank, &info );
907 trace3("zgelss result", info, rank, *z0_);
908 trace1("zgelss result", v0pm.matrix);
909 if(have_op){ untested();
910 gsl_vector_complex_set(z0_,_n_states+input_dim, opcoeff);
911 }
912 trace1("zgelss result+op", *z0_);
913 err+= gsl_blas_zgemv( CblasNoTrans, one, EV_, z0_, zero, v0back_ );
914 trace2("zgelss EV_*z0_", *v0back_, info);
915 }
916
917 trace1("hmmm", *z0_);
918 trace1("hmmm", *grid_);
919
920 for(unsigned i=0; i<_n_states; ++i){
921 gsl_complex dotgrid = gsl_vector_complex_get(z0dot_,i);
922 double on = 1./norm(dotgrid);
923 double ln = gsl_vector_get(grid_, i);
924 if(on<ln){ incomplete();
925 // gsl_vector_set(grid_,i,on);
926 }
927 }
928
929 { // refine state grid. using nodc-sens
930 for(unsigned i=0; i<_n_states; ++i){
931 gsl_vector_complex_const_view state = gsl_matrix_complex_const_column(EV_, i);
932 for(unsigned j=0; j<_n_inputs; ++j){
933 gsl_vector_complex_const_view input = gsl_matrix_complex_const_column(EV_, _n_states + j);
934 gsl_complex cpl;
935 gsl_blas_zdotc (&state.vector, &input.vector, &cpl);
936
937 trace2("nodc cpl", input, state);
938 trace3("nodc cpl", i, j , cpl);
939 trace1("nodc cpl", _step[j]);
940
941 double cpln = norm(cpl);
942
943 // set state grid to <= inputgrid/cpln
944 if(cpln * gsl_vector_get(grid_,i) > gsl_vector_get(grid_,j+_n_inputs)){ itested();
945 trace2("nodc reducing grid", i, gsl_vector_get(grid_,j+_n_inputs)/cpln);
946 gsl_vector_set(grid_,i,gsl_vector_get(grid_,j+_n_inputs)/cpln);
947 }else{
948 trace0("nodc not reducing grid");
949 }
950 }
951 }
952 }
953
954
955 // unsigned gridbase=2; not yet.
956 index_t state_index[_n_states];
957 if(_n_states){ itested();
958 // gsl_vector_complex_view z0dot = gsl_vector_complex_subvector(z0dot_,0,_n_states);
959 gsl_vector_const_view tol = gsl_vector_const_subvector(grid_,0,_n_states);
960 gsl_vector_complex_view z0 = gsl_vector_complex_subvector(z0_, 0, _n_states);
961 gsl_vector_complex_const_view lambda = gsl_vector_complex_const_subvector(lambda_, 0, _n_states);
962 trace2("quant tol", tol, lambda);
963 trace1("quant input", z0);
964 if (_quantize_states) { itested();
965 quantize(z0, lambda, tol, state_index);
966 }else{ untested();
967 to_grid(z0, tol.vector, state_index);
968 }
969 }
970
971 trace1("refined", *grid_);
972 to_grid(*z0_, *grid_);
973 trace1("discrete", *z0_);
974 trace1("discretized\n", *EV_);
975 err+= gsl_blas_zgemv( CblasNoTrans, one, EV_, z0_, zero, v0back_ );
976 trace1("v0back", *v0back_);
977
978 for(unsigned i=0; i<n; ++i){
979 _sim->_v0[i+1] = GSL_REAL(*gsl_vector_complex_const_ptr(v0back_,i));
980 }
981
982 // FIXME: merge insert/discretize. use unit grid
983 SSP_STATE* state;
984 if(_n_states) { untested();
985 gsl_vector_view r = gsl_vector_complex_real(z0_);
986 gsl_vector_view i = gsl_vector_complex_imag(z0_);
987 gsl_vector_div(&r.vector, grid_);
988 gsl_vector_div(&i.vector, grid_);
989
990 state = &chart->insert_state(state_index, _n_states);
991
992 }else{ itested();
993 state = &chart->insert();
994 }
995 USE(state);
996 trace1("inserted", chart->size());
997
998 trace1("========= ev_solver done ========\n", _sim->_acx);
999
1000 // BUG use global alloc/free the other stuff!
1001 gsl_eigen_genv_free(w);
1002 return *state;
1003 }
1004 /*--------------------------------------------------------------------------*/
1005 void EV_BASE::options(CS& Cmd, int Nest)
1006 {
1007 trace1("EV_BASE::options(CS&, int)", Nest);
1008
1009 _quantize_states = true;
1010 _loop[Nest] = _reverse_in[Nest] = false;
1011 _sim->_uic = false;
1012 _tran_step = OPT::dtddc;
1013 unsigned here = Cmd.cursor();
1014 do{
1015 ONE_OF
1016 || (Cmd.match1("'\"({") && ((Cmd >> _step_in[Nest]), (_stepmode[Nest] = LIN_STEP)))
1017 || (Cmd.is_float() && ((Cmd >> _step_in[Nest]), (_stepmode[Nest] = LIN_STEP)))
1018 || (Get(Cmd, "*", &_step_in[Nest]) && (_stepmode[Nest] = TIMES))
1019 || (Get(Cmd, "+", &_step_in[Nest]) && (_stepmode[Nest] = LIN_STEP))
1020 || (Get(Cmd, "by", &_step_in[Nest]) && (_stepmode[Nest] = LIN_STEP))
1021 || (Get(Cmd, "trstep", &_tran_step))
1022 || (Get(Cmd, "uic", &_uic))
1023 || (Get(Cmd, "qu{antize}", &_quantize_states))
1024 || (Get(Cmd, "step", &_step_in[Nest]) && (_stepmode[Nest] = LIN_STEP))
1025 || (Get(Cmd, "d{ecade}", &_step_in[Nest]) && (_stepmode[Nest] = DECADE))
1026 || (Get(Cmd, "ti{mes}", &_step_in[Nest]) && (_stepmode[Nest] = TIMES))
1027 || (Get(Cmd, "lin", &_step_in[Nest]) && (_stepmode[Nest] = LIN_PTS))
1028 || (Get(Cmd, "o{ctave}", &_step_in[Nest]) && (_stepmode[Nest] = OCTAVE))
1029 || Get(Cmd, "c{ontinue}", &_cont)
1030 || Get(Cmd, "tr{s}", &_do_tran_step)
1031 || Get(Cmd, "sr", &_slew_rate)
1032 || Get(Cmd, "dm", &_dump_matrix)
1033 || Get(Cmd, "dt{emp}", &temp_c_in, mOFFSET, OPT::temp_c)
1034 || Get(Cmd, "lo{op}", &_loop[Nest])
1035 || Get(Cmd, "re{verse}", &_reverse_in[Nest])
1036 || Get(Cmd, "te{mperature}",&temp_c_in)
1037 || (Cmd.umatch("tr{ace} {=}") &&
1038 (ONE_OF
1039 || Set(Cmd, "n{one}", &_trace, tNONE)
1040 || Set(Cmd, "o{ff}", &_trace, tNONE)
1041 || Set(Cmd, "w{arnings}", &_trace, tUNDER)
1042 || Set(Cmd, "i{terations}",&_trace, tITERATION)
1043 || Set(Cmd, "v{erbose}", &_trace, tVERBOSE)
1044 || Cmd.warn(bWARNING,
1045 "need none, off, warnings, iterations, verbose")
1046 )
1047 )
1048 || _out.outset(Cmd);
1049 }while (Cmd.more() && !Cmd.stuck(&here));
1050
1051 }
1052 /*--------------------------------------------------------------------------*/
1053 void EV_BASE::set_uic_caps_constant(bool x)
1054 {
1055 for( vector<STORAGE*>::iterator i=_uic_caplist.begin(); i!=_uic_caplist.end(); ++i)
1056 {
1057 (*i)->set_constant(x);
1058 }
1059 }
1060 /*--------------------------------------------------------------------------*/
1061 // void EV_BASE::compute_op(OPT::ITL itl)
1062 void EV_BASE::get_op(OPT::ITL itl)
1063 { itested();
1064 _sim->_time0 = _sim->_dt0 = 0.0;
1065 _sim->tr_reset();
1066 _sim->_phase = p_INIT_DC; // BUG?
1067
1068 _sim->_loadq.clear();
1069 _sim->_evalq->clear();
1070 _sim->_evalq_uc->clear();
1071 _sim->set_inc_mode_bad();
1072
1073 TRACE oldtrace=_trace;
1074 _trace=tNONE;
1075 int converged = solve_with_homotopy(itl,_trace);
1076 _trace=oldtrace;
1077
1078 size_t d = _sim->_total_nodes;
1079
1080
1081 if (!converged) {itested();
1082 error(bWARNING, "did not converge\n");
1083 throw(Exception("foobar"));
1084 }
1085 ::status.accept.start();
1086 _sim->set_limit();
1087
1088 if(_dump_matrix){
1089 _out << " ======================== \n";
1090 _out << "_i ( " << _sim->_i[1];
1091 for(unsigned a=2; a <= d; ++a){
1092 _out << " " << _sim->_i[a];
1093 }
1094 _out << ") \n";
1095 }
1096
1097 _sim->_loadq.clear();
1098
1099 _sim->_uic = false;
1100 CARD_LIST::card_list.tr_accept(); // hmmm.
1101 CARD_LIST::card_list.do_tr(); // no-uic evaluates caps
1102 _sim->_phase = p_INIT_DC;
1103 _sim->_uic = false;
1104 _sim->_evalq->clear();
1105 _sim->_evalq_uc->clear();
1106 _sim->count_iterations(iTOTAL);
1107 ::status.accept.stop();
1108 rhs_ = gsl_vector_alloc(d);
1109 array_to_gsl_vector(*rhs_,_sim->_i+1,d);
1110
1111 #if 0 // HACK, somehow use inc mode?!
1112 _sim->set_inc_mode_no();
1113 clear_arrays();
1114 CARD_LIST::card_list.tr_load();
1115 #else
1116 while (!_sim->_loadq.empty()) {
1117 trace1("loading from q", _sim->_loadq.back()->long_label());
1118 _sim->_loadq.back()->tr_load();
1119 _sim->_loadq.pop_back();
1120 }
1121 #endif
1122 rhs_ = gsl_vector_alloc(d);
1123 array_to_gsl_vector(*rhs_,_sim->_i+1,d);
1124 ::status.lud.start();
1125 _sim->_lu.lu_decomp(_sim->_aa); // , bool(OPT::lubypass && _sim->is_inc_mode()));
1126 ::status.lud.stop();
1127
1128
1129 trace1("fresh", *rhs_);
1130 trace1("fresh\n", _sim->_aa);
1131 _sim->_lu.fbsub(_sim->_i);
1132 array_to_gsl_vector(*rhs_,_sim->_i+1,d);
1133 trace1("G^-1 rhs", *rhs_);
1134
1135 // ???
1136 finish_building_evalq();
1137 evaluate_models();
1138 // CARD_LIST::card_list.tr_load();
1139 // array_to_gsl_vector(*rhs_,_sim->_i+1,d);
1140
1141 _sim->keep_voltages(); // vdc = v0
1142
1143 _sim->_acx.reallocate(); // ?
1144 _sim->_jomega = COMPLEX(0., 1.);
1145 _sim->_phase = p_AC;
1146 CARD_LIST::card_list.ac_begin();
1147
1148 if (_dump_matrix) {
1149 _out << "solved w/ht\n";
1150 _out << "i ( " << _sim->_i[1]; // K-put
1151
1152 for(unsigned a=2; a <= _sim->_total_nodes; ++a){
1153 _out << " " << _sim->_i[a];
1154 }
1155 _out << ") \n";
1156 _out << "v0 = ( " << _sim->_v0[1];
1157 for(unsigned a=2;a <= _sim->_total_nodes; ++a){
1158 _out << " " << _sim->_v0[a];
1159 }
1160 _out << ") \n";
1161 }
1162
1163 _sim->_uic = false;
1164 _sim->set_inc_mode_yes();
1165 while (!_sim->_loadq.empty()) {
1166 trace1("loading from q", _sim->_loadq.back()->long_label());
1167 _sim->_loadq.back()->tr_load();
1168 _sim->_loadq.pop_back();
1169 }
1170 }
1171 /*--------------------------------------------------------------------------*/
1172 void EV_BASE::get_op_()
1173 {
1174 _sim->_time0 = 0.;
1175 _sim->set_inc_mode_no(); // uaargh
1176 _sim->_uic = false;
1177 _sim->_phase = p_RESTORE;
1178 size_t n = _sim->_total_nodes;
1179
1180 gsl_vector_view rhs__ = gsl_vector_view_array(_sim->_i+1,n);
1181 trace1("get_op b4", rhs__);
1182 clear_arrays();
1183
1184 // CARD_LIST::card_list.tr_accept(); // hmmm.
1185 CARD_LIST::card_list.tr_begin(); // hmmm.
1186 // finish_building_evalq();
1187 CARD_LIST::card_list.do_tr();
1188
1189 #if 1 // hack
1190 CARD_LIST::card_list.tr_load();
1191 #else
1192
1193 while (!_sim->_loadq.empty()) {
1194 trace1("get_op loading from q", _sim->_loadq.back()->long_label());
1195 _sim->_loadq.back()->tr_load();
1196 _sim->_loadq.pop_back();
1197 }
1198 #endif
1199
1200 assert(rhs_->size == n);
1201 array_to_gsl_vector(*rhs_,_sim->_i+1,n);
1202 ::status.lud.start();
1203 _sim->_lu.lu_decomp(_sim->_aa); // , bool(OPT::lubypass && _sim->is_inc_mode()));
1204 ::status.lud.stop();
1205
1206 trace1("get_op fresh", *rhs_);
1207 trace1("get_op fresh\n", _sim->_aa);
1208 _sim->_lu.fbsub(_sim->_i);
1209 array_to_gsl_vector(*rhs_,_sim->_i+1,n);
1210 trace1("G^-1 rhs", *rhs_);
1211
1212 { // also fetch AC stuff
1213 _sim->_acx.zero();
1214 std::fill_n(_sim->_ac, n+1, 0.);
1215 ::status.load.start();
1216 _sim->count_iterations(iTOTAL);
1217 CARD_LIST::card_list.do_ac();
1218 CARD_LIST::card_list.ac_load();
1219 ::status.load.stop();
1220 }
1221 }
1222 /*--------------------------------------------------------------------------*/
1223 void EV::sweep_recursive(int Nest)
1224 {
1225 trace1("EV_BASE::sweep_recursive(int)", Nest);
1226 size_t d = _sim->_total_nodes;
1227 unsigned n = _sim->_total_nodes;
1228 --Nest;
1229 assert(Nest >= 0);
1230 assert(Nest < DCNEST);
1231
1232 // double iddc[d];
1233
1234 OPT::ITL itl = OPT::DCBIAS;
1235
1236 first(Nest);
1237 do {
1238 trace1("EV_BASE::sweep_recursive loopstart", Nest);
1239 _sim->_temp_c = temp_c_in;
1240 if (Nest == 0) {
1241 _sim->_time0 = _sim->_dt0 = 0.0;
1242 _sim->tr_reset();
1243 // _sim->zero_currents();
1244 //
1245 // why not?
1246 _sim->_uic = true;
1247 _sim->_phase = p_INIT_DC;
1248
1249 _sim->_loadq.clear();
1250 _sim->_evalq->clear();
1251 _sim->_evalq_uc->clear();
1252 trace0("hot");
1253 _sim->set_inc_mode_bad();
1254
1255 _sim->_uic = true; // ?!
1256 get_op(itl);
1257
1258 { // fetch ACX
1259 trace0("AC::solve");
1260 _sim->_acx.zero();
1261 std::fill_n(_sim->_ac, _sim->_total_nodes+1, 0.);
1262 ::status.load.start();
1263 _sim->count_iterations(iTOTAL);
1264 CARD_LIST::card_list.do_ac();
1265 CARD_LIST::card_list.ac_load();
1266 ::status.load.stop();
1267 }
1268
1269 // if verbose outdata(1./0.);
1270 SSP_CHART* chart = NULL;
1271 SSP_STATE& state = ev_solver(chart);
1272 gsl_vector_view v0__ = gsl_vector_view_array(_sim->_v0+1,n);
1273 SSP_SPL& spl = state.insert_spl(v0__, chart);
1274 USE(spl);
1275 USE(state);
1276
1277 gsl_matrix_complex_const_view dcs = gsl_matrix_complex_const_submatrix (EV_, 0, _n_states, d, d -_n_states);
1278 USE(dcs);
1279 // state.insert_op(dcs);
1280
1281 { // some more AC stuff
1282
1283 if(_dump_matrix){ untested();
1284 _out << "RS ( " << _Gu[0];
1285 for(unsigned a=1; a < d; ++a){
1286 _out << " " << _Gu[a];
1287 }
1288 _out << ") \n";
1289 }
1290
1291 // irgendwoher di/du holen...
1292 // C.fbsub( dv, Gu , dv );
1293
1294 }
1295
1296 if (spl.is_new()){
1297 EV_BASE::_out << chart->id();
1298 outdata(spl.id());
1299 }
1300
1301 // here??
1302 //CARD_LIST::card_list.tr_regress(); incomplete(); // => do_tran_step ?
1303 itl = OPT::DCXFER;
1304 }else{
1305 sweep_recursive(Nest);
1306 }
1307 } while (next(Nest));
1308 }
1309 /*--------------------------------------------------------------------------*/
1310 void EV_BASE::first(int Nest)
1311 {
1312 trace2("EV_BASE::first", Nest, _start[Nest]);
1313 assert(Nest >= 0);
1314 assert(Nest < DCNEST);
1315 assert(_start);
1316 assert(_sweepval);
1317 // assert(_pushel[Nest]);
1318
1319 (*_sweepval[Nest]) = _start[Nest];
1320 CARD_LIST::card_list.precalc_last();
1321 if (ELEMENT* c = dynamic_cast<ELEMENT*>(_zap[Nest])) {
1322 c->set_constant(false); // because of extra precalc_last
1323 // obsolete, once pointer hack is fixed
1324 }
1325
1326 // here? (hack...)
1327 //if(_pushel[Nest])
1328 // _pushel[Nest]->set_ic(_start[Nest]);
1329 _reverse[Nest] = false;
1330 if (_reverse_in[Nest]) {itested();
1331 while (next(Nest)) {itested();
1332 /* nothing */;
1333 }
1334 _reverse[Nest] = true;
1335 next(Nest);
1336 }else{
1337 }
1338 _sim->_phase = p_INIT_DC;
1339 }
1340 /*--------------------------------------------------------------------------*/
1341 // EV_DISC
1342 bool EV_BASE::next(int Nest)
1343 {
1344 trace1("EV_BASE::next(int)", Nest);
1345
1346 bool ok = false;
1347 if (_linswp[Nest]) {
1348 double fudge = _step[Nest] / 10.;
1349 if (_step[Nest] == 0.) {
1350 ok = false;
1351 }else{
1352 if (!_reverse[Nest]) {
1353 *(_sweepval[Nest]) += _step[Nest];
1354 fixzero(_sweepval[Nest], _step[Nest]);
1355 CARD_LIST::card_list.precalc_last();
1356 if (ELEMENT* c = dynamic_cast<ELEMENT*>(_zap[Nest])) {
1357 c->set_constant(false); // because of extra precalc_last
1358 // obsolete, once pointer hack is fixed
1359 }
1360 ok=in_order(_start[Nest]-fudge,(*_sweepval[Nest]),_stop[Nest]+fudge);
1361 //_pushel[Nest]->set_ic(_sweepval[Nest]);
1362 if (!ok && _loop[Nest]) { untested();
1363 _reverse[Nest] = true;
1364 }else{
1365 }
1366 }else{ untested();
1367 }
1368 if (_reverse[Nest]) { untested();
1369 *(_sweepval[Nest]) -= _step[Nest];
1370 CARD_LIST::card_list.precalc_last();
1371 if (ELEMENT* c = dynamic_cast<ELEMENT*>(_zap[Nest])) {
1372 c->set_constant(false); // because of extra precalc_last
1373 // obsolete, once pointer hack is fixed
1374 }
1375 fixzero(_sweepval[Nest], _step[Nest]);
1376 ok=in_order(_start[Nest]-fudge,(*_sweepval[Nest]),_stop[Nest]+fudge);
1377 //_pushel[Nest]->set_ic(_sweepval[Nest]);
1378 }else{
1379 }
1380 }
1381 }else{ untested();
1382 double fudge = pow(_step[Nest], .1);
1383 if (_step[Nest] == 1.) {untested();
1384 ok = false;
1385 }else{ untested();
1386 if (!_reverse[Nest]) { untested();
1387 *(_sweepval[Nest]) *= _step[Nest];
1388 ok=in_order(_start[Nest]/fudge,*(_sweepval[Nest]),_stop[Nest]*fudge);
1389 //_pushel[Nest]->set_ic(_sweepval[Nest]);
1390 if (!ok && _loop[Nest]) {untested();
1391 _reverse[Nest] = true;
1392 }else{ untested();
1393 }
1394 }else{ untested();
1395 }
1396 if (_reverse[Nest]) {untested();
1397 *(_sweepval[Nest]) /= _step[Nest];
1398 ok=in_order(_start[Nest]/fudge,*(_sweepval[Nest]),_stop[Nest]*fudge);
1399 //_pushel[Nest]->set_ic(_sweepval[Nest]);
1400 }else{ untested();
1401 }
1402 }
1403 }
1404 _sim->_phase = p_DC_SWEEP;
1405 return ok;
1406 }
1407 /*-----------------------------------------------------------*/
1408 void EV_BASE::ac_snapshot()
1409 { itested(); // in sock
1410 trace0("EV_BASE::ac_snapshot");
1411
1412 _sim->_acx.reallocate();
1413 _sim->_jomega = COMPLEX(0., 1.0);
1414
1415 _sim->_phase = p_AC;
1416
1417 // _sim->_acx.set_min_pivot(OPT::pivtol);
1418 //
1419 CARD_LIST::card_list.ac_begin();
1420
1421 if(_dump_matrix){ untested();
1422 _out << "i ( " << _sim->_i[1];
1423 for(unsigned a=2; a <= _sim->_total_nodes; ++a){ untested();
1424 _out << " " << _sim->_i[a];
1425 }
1426 _out << ") \n";
1427 _out << "v0 = ( " << _sim->_v0[1];
1428 for(unsigned a=2;a <= _sim->_total_nodes; ++a){ untested();
1429 _out << " " << _sim->_v0[a];
1430 }
1431 _out << ") \n";
1432 }
1433
1434 // if verbose
1435 _sim->_uic = false;
1436
1437 _sim->_acx.zero();
1438 std::fill_n(_sim->_ac, _sim->_total_nodes+1, 0.);
1439
1440 ::status.load.start();
1441 _sim->count_iterations(iTOTAL);
1442 CARD_LIST::card_list.do_ac();
1443 CARD_LIST::card_list.ac_load();
1444 ::status.load.stop();
1445 }
1446 /*--------------------------------------------------------------------------*/
1447 /*--------------------------------------------------------------------------*/
1448 // take some eigenvectors (alpha!=0, beta==0)
1449 // rearrange considering input sensitivities.
1450 // obsolete?
1451 #if 0
1452 void EV::align_inf_eigenvectors(gsl_matrix_complex_view& EV_inf, size_t i)
1453 {
1454 if (!_zap[i]){ untested();
1455 return;
1456 }
1457 size_t n = _sim->_total_nodes;
1458 trace3("fixing\n", EV_inf.matrix, EV_inf.matrix.size1, EV_inf.matrix.size2);
1459
1460
1461
1462 trace1("lu", _sim->_acx);
1463
1464 assert(_zap[i]); // incomplete
1465 clear_sens();
1466 _zap[i]->sens_load("dummy");
1467 complex_array_to_gsl_vector_complex(*sens_,_sim->_sens+1,n);
1468 trace1("2 sensstamp", *sens_);
1469
1470 ::status.back.start();
1471 _sim->_acx.fbsub(_sim->_sens);
1472 ::status.back.stop();
1473
1474 complex_array_to_gsl_vector_complex(*sens_,_sim->_sens+1,n);
1475 trace2("sens", _zap[i]->long_label(), *sens_);
1476
1477 gsl_vector_complex_view first_EV = firstcol(EV_inf);
1478 gsl_vector_complex_memcpy(&first_EV.vector, sens_);
1479
1480
1481 if(EV_inf.matrix.size2 > 1) { untested();
1482 EV_inf = gsl_matrix_complex_submatrix (&EV_inf.matrix, 0, 1, n, EV_inf.matrix.size2-1);
1483 align_inf_eigenvectors(EV_inf, i+1);
1484 }
1485 }
1486 #endif
1487 /*--------------------------------------------------------------------------*/
1488 // output sensitivities. compute grid.
1489 // foreach output, use eigenvector -> output sensitivity
1490 // FIXME: use actual probes
1491 void EV_BASE::output_sens_grid(gsl_matrix_complex_view& relevant_EV, gsl_vector* grid_, const gsl_vector_complex* lambda)
1492 {
1493 int err = 0;
1494 trace1("output sens\n", relevant_EV.matrix);
1495 trace1("output sens", *grid_);
1496 size_t n = _sim->_total_nodes;
1497
1498 _output_iter = _output.begin();
1499
1500 for(; _output_iter!=_output.end(); ++_output_iter){
1501 set_sens(_output_iter);
1502 gsl_vector_complex_view ztmp_left = gsl_vector_complex_subvector(ztmp_, 0, relevant_EV.matrix.size2 );
1503
1504 complex_array_to_gsl_vector_complex(*sens_,_sim->_sens+1,n);
1505 trace2("sens stamp", _output_iter->label, *sens_);
1506 err+= gsl_blas_zgemv( CblasTrans, one, &relevant_EV.matrix, sens_, zero, &ztmp_left.vector );
1507 trace1("EV senstivities", *ztmp_);
1508
1509 for(unsigned i=0; i<relevant_EV.matrix.size2; ++i){
1510
1511 double ln = gsl_vector_get(grid_, i);
1512 gsl_complex lam = gsl_vector_complex_get(lambda, i);
1513 gsl_complex output_sens = gsl_vector_complex_get(ztmp_,i);
1514 trace3("grid", output_sens, lam, i);
1515
1516 double n = 1;
1517 if(i<_n_states){
1518 // n = .1/sqrt(OPT::dtmin)/norm(lam);
1519 // trace3("grid", output_sens, lam, n);
1520 gsl_vector_set(grid_,i,1./0.); // reset. grid has been abused above
1521 }
1522
1523 double on = norm(output_sens) * n;
1524
1525 if (0.!=on){ itested();
1526 double outgrid = sqrt(OPT::reltol)/on;
1527 // outgrid = (OPT::reltol)/on;
1528 if(outgrid < fabs(ln)){ itested();
1529 trace3("reducing grid", i, outgrid, fabs(ln));
1530 gsl_vector_set(grid_,i,outgrid);
1531 }else{
1532 trace2("not reducing grid", outgrid, fabs(ln));
1533 }
1534 }
1535 }
1536 }
1537 trace1("output sens done", *grid_);
1538 }
1539 /*--------------------------------------------------------------------------*/
1540 // quantize z0, stash index into index
1541 // z0= \pm (2**(\pm index) - 1) * tol
1542 // index = floor ( ln (z0 + reltol/ reltol) / ln(2) )
1543 void EV_BASE::quantize(gsl_vector_complex_view z0, gsl_vector_complex_const_view /*lambda*/, gsl_vector_const_view tol, index_t* index)
1544 {
1545 trace1("inp", z0);
1546 for(unsigned i=0; i<_n_states; ++i) {
1547 gsl_complex z = gsl_vector_complex_get(z0,i);
1548 double t_ = gsl_vector_get(tol,i);
1549 gsl_complex t = gsl_complex_rect(t_,0); // hmmm...
1550 int sign = 1-2* int(bool(signbit(GSL_REAL(z))));
1551 trace2("log discretized", z, sign);
1552 z = gsl_complex_mul_real(z, double(sign));
1553
1554 z = gsl_complex_div(z, t);
1555 z = gsl_complex_add(z,gsl_complex_rect(1,0));
1556 z = gsl_complex_log(z);
1557
1558 double ln2 = logf(2.);
1559 z = gsl_complex_add( gsl_complex_rect(ln2 * .5, 0. ),z) ;
1560 z = gsl_complex_div( z, gsl_complex_rect(ln2 ,0.));
1561
1562 if (GSL_REAL(z) < 0) { untested();
1563 sign = 0;
1564 }
1565
1566 index[i] = sign*int(floor(GSL_REAL(z)));
1567 GSL_SET_REAL(&z, sign * GSL_REAL(z));
1568
1569 if (index[i]){
1570 z = gsl_complex_rect( sign*t_ * (pow(2,sign*index[i]-1)), 0 );
1571 } else {
1572 z = zero;
1573 }
1574 assert(fabs(norm(z))<10); // for now.
1575 gsl_vector_complex_set(&z0.vector,i,z);
1576 }
1577 trace2("quantized", z0, index[0]);
1578 // if err; throw...
1579 }
1580 /*--------------------------------------------------------------------------*/
1581 string EV_BASE::status()const
1582 {
1583 return string("state space: charts=")
1584 + to_string(SSP_TREE::num_charts())
1585 + ", states=" + to_string(SSP_TREE::num_states())
1586 + ", samples=" + to_string(SSP_TREE::num_spls()) + "\n";
1587 }
1588 /*--------------------------------------------------------------------------*/
1589 /*--------------------------------------------------------------------------*/