From ad23146544ebb3e5f8477ffca4e37c032c26a929 Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Sat, 9 Sep 2023 18:37:13 +0200 Subject: [PATCH 01/13] Don't check continuously for autostop, only when option flag is set This speeds up simulations with simple device evaluation, but many time stepps. --- src/frontend/measure.c | 2 +- src/spicelib/analysis/dctran.c | 18 ++++++++++++++++-- 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/src/frontend/measure.c b/src/frontend/measure.c index c91715a11..b7b6680cc 100644 --- a/src/frontend/measure.c +++ b/src/frontend/measure.c @@ -532,7 +532,7 @@ check_autostop(char* what) { bool flag = FALSE; - if (cp_getvar("autostop", CP_BOOL, NULL, 0)) +// if (cp_getvar("autostop", CP_BOOL, NULL, 0)) flag = do_measure(what, TRUE); return flag; diff --git a/src/spicelib/analysis/dctran.c b/src/spicelib/analysis/dctran.c index 5f7d66c60..d3f6f8c12 100644 --- a/src/spicelib/analysis/dctran.c +++ b/src/spicelib/analysis/dctran.c @@ -94,6 +94,8 @@ DCtran(CKTcircuit *ckt, int numNames; double maxstepsize = 0.0; + bool have_autostop = FALSE, flag_autostop = FALSE; + int ltra_num; CKTnode *node; #ifdef XSPICE @@ -184,6 +186,8 @@ DCtran(CKTcircuit *ckt, save_mode = (ckt->CKTmode&MODEUIC) | MODETRANOP | MODEINITJCT; save_order = ckt->CKTorder; + have_autostop = cp_getvar("autostop", CP_BOOL, NULL, 0); + /* Add breakpoints here which have been requested by the user setting the stop command as 'stop when time = xx'. Get data from the global dbs data base. @@ -466,8 +470,14 @@ DCtran(CKTcircuit *ckt, /* gtri - end - wbk - Update event queues/data for accepted timepoint */ #endif ckt->CKTstat->STAToldIter = ckt->CKTstat->STATnumIter; - if (check_autostop("tran") || - ckt->CKTfinalTime - ckt->CKTtime < ckt->CKTminBreak) { + /* check for the end of the tran simulation, either by< stop time given, + or final time has been reached. */ + if (have_autostop) + /* time consuming autostop check only, when variable 'autostop' has been set + before tran is started.*/ + flag_autostop = check_autostop("tran"); + /* If CKTtime and CKTfinalTime are almost equal, then finish */ + if (flag_autostop || AlmostEqualUlps(ckt->CKTtime, ckt->CKTfinalTime, 100)) { #ifdef STEPDEBUG printf(" done: time is %g, final time is %g, and tol is %g\n", ckt->CKTtime, ckt->CKTfinalTime, ckt->CKTminBreak); @@ -480,6 +490,10 @@ DCtran(CKTcircuit *ckt, ckt->CKTsenInfo->SENmode = save; } #endif + if (flag_autostop) + fprintf(stdout, "\nNote: Autostop after %e s, all measurement conditions are fulfilled\n", ckt->CKTtime); + + /* Final return from tran*/ return(OK); } if(SPfrontEnd->IFpauseTest()) { From 7a646c0a12313b346de90b365673381f9e6729c9 Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Sat, 9 Sep 2023 18:38:09 +0200 Subject: [PATCH 02/13] If 'strict_errorhandling' is set, bail out if operating point is not found --- src/spicelib/analysis/cktop.c | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/spicelib/analysis/cktop.c b/src/spicelib/analysis/cktop.c index 1db078981..6bebb681a 100644 --- a/src/spicelib/analysis/cktop.c +++ b/src/spicelib/analysis/cktop.c @@ -10,6 +10,7 @@ Modified: 2005 Paolo Nenzi - Restructured #include "ngspice/devdefs.h" #include "ngspice/sperror.h" #include "ngspice/cpextern.h" +#include "ngspice/fteext.h" #ifdef XSPICE #include "ngspice/enh.h" @@ -108,6 +109,8 @@ CKTop (CKTcircuit *ckt, long int firstmode, long int continuemode, #endif fprintf(cp_err, "\nError: The operating point could not be simulated successfully.\n"); + if (ft_stricterror) + controlled_exit(1); fprintf(cp_err, " Any of the following steps may fail.!\n\n"); return converged; From 492bb64d9213fdec249580e6767f6e8b1f693281 Mon Sep 17 00:00:00 2001 From: Brian Taylor Date: Sun, 3 Sep 2023 13:16:44 -0700 Subject: [PATCH 03/13] By default, use the shortest typical delay estimate. This makes the digi_74LS90_74LS42.cir testcase for bug641 behave almost the same as MicroCap 12. In ngspice and MicroCap, the only signal with a glitch is not_y8. The other not_* signals look the same. Setting ps_use_mntymx in .spiceinit will change the delay estimates. See the function set_u_devices_info in src/frontend/udevices.c for the various settings of ps_use_mntymx. --- src/frontend/logicexp.c | 80 +++++++++++--- src/frontend/udevices.c | 185 +++++++++++++++++++++++---------- src/include/ngspice/udevices.h | 6 ++ 3 files changed, 201 insertions(+), 70 deletions(-) diff --git a/src/frontend/logicexp.c b/src/frontend/logicexp.c index 5911da317..32626be0e 100644 --- a/src/frontend/logicexp.c +++ b/src/frontend/logicexp.c @@ -1912,7 +1912,38 @@ static char *get_typ_estimate(char *min, char *typ, char *max, DSTRING *pds) return NULL; } -static char *typical_estimate(char *delay_str, DSTRING *pds) +static char *get_one_estimate(char *s, DSTRING *pds) +{ + ds_clear(pds); + if (s && strlen(s) > 0 && s[0] != '-') { + ds_cat_str(pds, s); + return ds_get_buf(pds); + } else { + return NULL; + } +} + +static char *get_delay_estimate(char *min, char *typ, char *max, DSTRING *pds) +{ + char *one = NULL; + struct udevices_info info = u_get_udevices_info(); + int delay_type = info.mntymx; + if (delay_type == 1) { // min + one = get_one_estimate(min, pds); + if (one) { + return one; + } + } else if (delay_type == 2) { // max + one = get_one_estimate(max, pds); + if (one) { + return one; + } + } + // typ + return get_typ_estimate(min, typ, max, pds); +} + +static char *mntymx_estimate(char *delay_str, DSTRING *pds) { /* Input string (t1,t2,t2) */ int which = 0; @@ -1945,7 +1976,7 @@ static char *typical_estimate(char *delay_str, DSTRING *pds) break; } } - s = get_typ_estimate(ds_get_buf(&dmin), ds_get_buf(&dtyp), + s = get_delay_estimate(ds_get_buf(&dmin), ds_get_buf(&dtyp), ds_get_buf(&dmax), pds); ds_free(&dmin); ds_free(&dtyp); @@ -1967,22 +1998,25 @@ static BOOL extract_delay( BOOL in_delay = FALSE, ret_val = TRUE; int i; BOOL prit = PRINT_ALL; - float typ_max_val = 0.0, typ_val = 0.0; + BOOL shorter = FALSE, update_val = FALSE; + struct udevices_info info = u_get_udevices_info(); + float del_max_val = 0.0, del_val = 0.0, del_min_val = FLT_MAX; char *units; + shorter = info.shorter_delays; DS_CREATE(dly, 64); - DS_CREATE(dtyp_max_str, 16); + DS_CREATE(ddel_str, 16); DS_CREATE(tmp_ds, 128); if (val != '=') { ds_free(&dly); - ds_free(&dtyp_max_str); + ds_free(&ddel_str); ds_free(&tmp_ds); return FALSE; } val = lexer_scan(lx); if (val != '{') { ds_free(&dly); - ds_free(&dtyp_max_str); + ds_free(&ddel_str); ds_free(&tmp_ds); return FALSE; } @@ -2005,7 +2039,7 @@ static BOOL extract_delay( char *tmps; ds_clear(&tmp_ds); in_delay = FALSE; - tmps = typical_estimate(ds_get_buf(&dly), &tmp_ds); + tmps = mntymx_estimate(ds_get_buf(&dly), &tmp_ds); if (!tmps) { ret_val = FALSE; ds_clear(&tmp_ds); @@ -2015,22 +2049,34 @@ static BOOL extract_delay( printf("%s\n", ds_get_buf(&dly)); printf("estimate \"%s\"\n", tmps); } - typ_val = strtof(tmps, &units); - if (typ_val > typ_max_val) { + del_val = strtof(tmps, &units); + update_val = FALSE; + if (shorter) { + if (del_val < del_min_val) { + update_val = TRUE; + } + } else if (del_val > del_max_val) { + update_val = TRUE; + } + if (update_val) { ds_clear(&delay_string); - ds_clear(&dtyp_max_str); - ds_cat_str(&dtyp_max_str, tmps); - typ_max_val = typ_val; - if (ds_get_length(&dtyp_max_str) > 0) { + ds_clear(&ddel_str); + ds_cat_str(&ddel_str, tmps); + if (shorter) { + del_min_val = del_val; + } else { + del_max_val = del_val; + } + if (ds_get_length(&ddel_str) > 0) { if (tri) { ds_cat_printf(&delay_string, "(inertial_delay=true delay=%s)", - ds_get_buf(&dtyp_max_str)); + ds_get_buf(&ddel_str)); } else { ds_cat_printf(&delay_string, "(inertial_delay=true rise_delay=%s fall_delay=%s)", - ds_get_buf(&dtyp_max_str), - ds_get_buf(&dtyp_max_str)); + ds_get_buf(&ddel_str), + ds_get_buf(&ddel_str)); } } else { printf("WARNING pindly DELAY not found\n"); @@ -2055,7 +2101,7 @@ static BOOL extract_delay( val = lexer_scan(lx); } // end while != '}' ds_free(&dly); - ds_free(&dtyp_max_str); + ds_free(&ddel_str); ds_free(&tmp_ds); return ret_val; } diff --git a/src/frontend/udevices.c b/src/frontend/udevices.c index c06184a01..f4433e70a 100644 --- a/src/frontend/udevices.c +++ b/src/frontend/udevices.c @@ -307,6 +307,10 @@ static int ps_udevice_exit = 0; static int ps_tpz_delays = 0; // For tristate delays static int ps_with_inverters = 0; // For ff/latch control inputs static int ps_with_tri_inverters = 0; // For inv3/inv3a data inputs + +static int ps_use_mntymx = 0; +static struct udevices_info mntymx_shorter = {0, FALSE}; + static NAME_ENTRY new_names_list = NULL; static NAME_ENTRY input_names_list = NULL; static NAME_ENTRY output_names_list = NULL; @@ -319,6 +323,46 @@ static BOOL add_drive_hilo = FALSE; static char *current_subckt = NULL; static unsigned int subckt_msg_count = 0; +static void set_u_devices_info(int mntymx_duration) +{ + switch (mntymx_duration) { + case 0: // typ + longer delays + mntymx_shorter.mntymx = 0; + mntymx_shorter.shorter_delays = FALSE; + break; + case 1: // min + longer delays + mntymx_shorter.mntymx = 1; + mntymx_shorter.shorter_delays = FALSE; + break; + case 2: // max + longer delays + mntymx_shorter.mntymx = 2; + mntymx_shorter.shorter_delays = FALSE; + break; + case 4: // typ + shorter delays + mntymx_shorter.mntymx = 0; + mntymx_shorter.shorter_delays = TRUE; + break; + case 5: // min + shorter delays + mntymx_shorter.mntymx = 1; + mntymx_shorter.shorter_delays = TRUE; + break; + case 6: // max + shorter delays + mntymx_shorter.mntymx = 2; + mntymx_shorter.shorter_delays = TRUE; + break; + default: + mntymx_shorter.mntymx = 0; + mntymx_shorter.shorter_delays = FALSE; + break; + } + return; +} + +struct udevices_info u_get_udevices_info(void) +{ + return mntymx_shorter; +} + static void check_name_unused(char *name) { if (find_name_entry(name, new_names_list)) { @@ -897,7 +941,7 @@ void initialize_udevice(char *subckt_line) if they are the only delays. */ if (!cp_getvar("ps_tpz_delays", CP_NUM, &ps_tpz_delays, 0)) { - ps_tpz_delays = 0; + ps_tpz_delays = 1; // default: use tpz... delays if necessary } /* If non-zero use inverters with ff/latch control inputs */ if (!cp_getvar("ps_with_inverters", CP_NUM, &ps_with_inverters, 0)) { @@ -907,6 +951,11 @@ void initialize_udevice(char *subckt_line) if (!cp_getvar("ps_with_tri_inverters", CP_NUM, &ps_with_tri_inverters, 0)) { ps_with_tri_inverters = 0; } + if (!cp_getvar("ps_use_mntymx", CP_NUM, &ps_use_mntymx, 0)) { + ps_use_mntymx = 4; // default: typ + shorter delays + } + set_u_devices_info(ps_use_mntymx); + if (subckt_line && strncmp(subckt_line, ".subckt", 7) == 0) { add_all_port_names(subckt_line); current_subckt = TMALLOC(char, strlen(subckt_line) + 1); @@ -2759,10 +2808,32 @@ static void estimate_typ(struct timing_data *tdp) return; } +static void estimate_delay(struct timing_data *tdp) +{ + char *del = NULL; + if (!tdp) { return; } + if (mntymx_shorter.mntymx == 1) { // use min delay + del = tdp->min; + if (del && strlen(del) > 0 && del[0] != '-') { + tdp->estimate = EST_MIN; + return; + } + } else if (mntymx_shorter.mntymx == 2) { // use max delay + del = tdp->max; + if (del && strlen(del) > 0 && del[0] != '-') { + tdp->estimate = EST_MAX; + return; + } + } + // use typ delay + estimate_typ(tdp); + return; +} + static char *get_estimate(struct timing_data *tdp) { /* - Call after estimate_typ. + Call after estimate_delay. Don't call delete_timing_data until you have copied or finished with this return value. */ @@ -2777,8 +2848,11 @@ static char *get_estimate(struct timing_data *tdp) return NULL; } -static char *larger_delay(char *delay1, char *delay2) +static char *select_delay(char *delay1, char *delay2) { + /* Return the shorter or longer delay + depending on the mntymx_shorter.shorter_delays setting + */ float val1, val2; char *units1, *units2; @@ -2787,11 +2861,16 @@ static char *larger_delay(char *delay1, char *delay2) if (!eq(units1, units2)) { printf("WARNING units do not match\n"); } - if (val1 >= val2) { - return delay1; + if (mntymx_shorter.shorter_delays) { + if (val1 <= val2) { + return delay1; + } } else { - return delay2; + if (val1 >= val2) { + return delay1; + } } + return delay2; } /* NOTE @@ -2813,10 +2892,10 @@ static char *get_delays_ugate(char *rem) BOOL has_rising = FALSE, has_falling = FALSE; tdp1 = create_min_typ_max("tplh", rem); - estimate_typ(tdp1); + estimate_delay(tdp1); rising = get_estimate(tdp1); tdp2 = create_min_typ_max("tphl", rem); - estimate_typ(tdp2); + estimate_delay(tdp2); falling = get_estimate(tdp2); has_rising = (rising && strlen(rising) > 0); has_falling = (falling && strlen(falling) > 0); @@ -2845,7 +2924,7 @@ static char *get_delays_udly(char *rem) struct timing_data *tdp1; tdp1 = create_min_typ_max("dly", rem); - estimate_typ(tdp1); + estimate_delay(tdp1); udelay = get_estimate(tdp1); if (udelay) { delays = tprintf( @@ -2865,22 +2944,22 @@ static char *get_delays_utgate(char *rem) char *rising, *falling, *delays = NULL; struct timing_data *tdp1, *tdp2; struct timing_data *tdp3, *tdp4, *tdp5, *tdp6; - char *tplz, *tphz, *tpzl, *tpzh, *larger, *larger1, *larger2, *larger3; + char *tplz, *tphz, *tpzl, *tpzh, *select0, *select1, *select2, *select3; BOOL use_zdelays = FALSE; BOOL has_rising = FALSE, has_falling = FALSE; tdp1 = create_min_typ_max("tplh", rem); - estimate_typ(tdp1); + estimate_delay(tdp1); rising = get_estimate(tdp1); tdp2 = create_min_typ_max("tphl", rem); - estimate_typ(tdp2); + estimate_delay(tdp2); falling = get_estimate(tdp2); has_rising = (rising && strlen(rising) > 0); has_falling = (falling && strlen(falling) > 0); if (has_rising) { if (has_falling) { - larger = larger_delay(rising, falling); - delays = tprintf("(inertial_delay=true delay = %s)", larger); + select0 = select_delay(rising, falling); + delays = tprintf("(inertial_delay=true delay = %s)", select0); } else { delays = tprintf("(inertial_delay=true delay = %s)", rising); } @@ -2889,49 +2968,49 @@ static char *get_delays_utgate(char *rem) } else if (use_zdelays || (ps_tpz_delays & 1)) { /* No lh/hl delays, so try the largest lz/hz/zl/zh delay */ tdp3 = create_min_typ_max("tplz", rem); - estimate_typ(tdp3); + estimate_delay(tdp3); tplz = get_estimate(tdp3); tdp4 = create_min_typ_max("tphz", rem); - estimate_typ(tdp4); + estimate_delay(tdp4); tphz = get_estimate(tdp4); - larger1 = NULL; + select1 = NULL; if (tplz && strlen(tplz) > 0) { if (tphz && strlen(tphz) > 0) { - larger1 = larger_delay(tplz, tphz); + select1 = select_delay(tplz, tphz); } else { - larger1 = tplz; + select1 = tplz; } } else if (tphz && strlen(tphz) > 0) { - larger1 = tphz; + select1 = tphz; } tdp5 = create_min_typ_max("tpzl", rem); - estimate_typ(tdp5); + estimate_delay(tdp5); tpzl = get_estimate(tdp5); tdp6 = create_min_typ_max("tpzh", rem); - estimate_typ(tdp6); + estimate_delay(tdp6); tpzh = get_estimate(tdp6); - larger2 = NULL; + select2 = NULL; if (tpzl && strlen(tpzl) > 0) { if (tpzh && strlen(tpzh) > 0) { - larger2 = larger_delay(tpzl, tpzh); + select2 = select_delay(tpzl, tpzh); } else { - larger2 = tpzl; + select2 = tpzl; } } else if (tpzh && strlen(tpzh) > 0) { - larger2 = tpzh; + select2 = tpzh; } - larger3 = NULL; - if (larger1) { - if (larger2) { - larger3 = larger_delay(larger1, larger2); + select3 = NULL; + if (select1) { + if (select2) { + select3 = select_delay(select1, select2); } else { - larger3 = larger1; + select3 = select1; } - } else if (larger2) { - larger3 = larger2; + } else if (select2) { + select3 = select2; } - if (larger3) { - delays = tprintf("(inertial_delay=true delay = %s)", larger3); + if (select3) { + delays = tprintf("(inertial_delay=true delay = %s)", select3); } else { delays = tprintf("(inertial_delay=true delay=1.0e-12)"); } @@ -2951,26 +3030,26 @@ static char *get_delays_ueff(char *rem) { char *delays = NULL; char *clkqrise, *clkqfall, *pcqrise, *pcqfall; - char *clkd, *setd, *resetd, *larger; + char *clkd, *setd, *resetd, *select; struct timing_data *tdp1, *tdp2, *tdp3, *tdp4; tdp1 = create_min_typ_max("tpclkqlh", rem); - estimate_typ(tdp1); + estimate_delay(tdp1); clkqrise = get_estimate(tdp1); tdp2 = create_min_typ_max("tpclkqhl", rem); - estimate_typ(tdp2); + estimate_delay(tdp2); clkqfall = get_estimate(tdp2); tdp3 = create_min_typ_max("tppcqlh", rem); - estimate_typ(tdp3); + estimate_delay(tdp3); pcqrise = get_estimate(tdp3); tdp4 = create_min_typ_max("tppcqhl", rem); - estimate_typ(tdp4); + estimate_delay(tdp4); pcqfall = get_estimate(tdp4); clkd = NULL; if (clkqrise && strlen(clkqrise) > 0) { if (clkqfall && strlen(clkqfall) > 0) { - larger = larger_delay(clkqrise, clkqfall); - clkd = larger; + select = select_delay(clkqrise, clkqfall); + clkd = select; } else { clkd = clkqrise; } @@ -3018,7 +3097,7 @@ static char *get_delays_ugff(char *rem, char *d_name) char *delays = NULL, *dname; char *tpdqlh, *tpdqhl, *tpgqlh, *tpgqhl, *tppcqlh, *tppcqhl; char *d_delay, *enab, *setd, *resetd; - char *s1, *s2, *larger; + char *s1, *s2, *select; struct timing_data *tdp1, *tdp2, *tdp3, *tdp4, *tdp5, *tdp6; if (eq(d_name, "d_dlatch")) { @@ -3029,28 +3108,28 @@ static char *get_delays_ugff(char *rem, char *d_name) return NULL; } tdp1 = create_min_typ_max("tpdqlh", rem); - estimate_typ(tdp1); + estimate_delay(tdp1); tpdqlh = get_estimate(tdp1); tdp2 = create_min_typ_max("tpdqhl", rem); - estimate_typ(tdp2); + estimate_delay(tdp2); tpdqhl = get_estimate(tdp2); tdp3 = create_min_typ_max("tpgqlh", rem); - estimate_typ(tdp3); + estimate_delay(tdp3); tpgqlh = get_estimate(tdp3); tdp4 = create_min_typ_max("tpgqhl", rem); - estimate_typ(tdp4); + estimate_delay(tdp4); tpgqhl = get_estimate(tdp4); tdp5 = create_min_typ_max("tppcqlh", rem); - estimate_typ(tdp5); + estimate_delay(tdp5); tppcqlh = get_estimate(tdp5); tdp6 = create_min_typ_max("tppcqhl", rem); - estimate_typ(tdp6); + estimate_delay(tdp6); tppcqhl = get_estimate(tdp6); d_delay = NULL; if (tpdqlh && strlen(tpdqlh) > 0) { if (tpdqhl && strlen(tpdqhl) > 0) { - larger = larger_delay(tpdqlh, tpdqhl); - d_delay = larger; + select = select_delay(tpdqlh, tpdqhl); + d_delay = select; } else { d_delay = tpdqlh; } @@ -3060,8 +3139,8 @@ static char *get_delays_ugff(char *rem, char *d_name) enab = NULL; if (tpgqlh && strlen(tpgqlh) > 0) { if (tpgqhl && strlen(tpgqhl) > 0) { - larger = larger_delay(tpgqlh, tpgqhl); - enab = larger; + select = select_delay(tpgqlh, tpgqhl); + enab = select; } else { enab = tpgqlh; } diff --git a/src/include/ngspice/udevices.h b/src/include/ngspice/udevices.h index 5ee5a5fb1..7fe19a881 100644 --- a/src/include/ngspice/udevices.h +++ b/src/include/ngspice/udevices.h @@ -1,6 +1,11 @@ #ifndef INCLUDED_UDEVICES #define INCLUDED_UDEVICES +struct udevices_info { + int mntymx; + BOOL shorter_delays; +}; + BOOL u_process_instance(char *line); BOOL u_process_model_line(char *line); BOOL u_check_instance(char *line); @@ -10,5 +15,6 @@ void cleanup_udevice(void); void u_add_instance(char *str); void u_add_logicexp_model(char *tmodel, char *xspice_gate, char *model_name); void u_remember_pin(char *name, int type); +struct udevices_info u_get_udevices_info(void); #endif From 3b89410b8f222892065c6888adabb149fa54eb47 Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Sat, 9 Sep 2023 23:34:49 +0200 Subject: [PATCH 04/13] Slight cosmetics in comment. --- src/spicelib/analysis/dctran.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spicelib/analysis/dctran.c b/src/spicelib/analysis/dctran.c index d3f6f8c12..fcd3f5240 100644 --- a/src/spicelib/analysis/dctran.c +++ b/src/spicelib/analysis/dctran.c @@ -491,7 +491,7 @@ DCtran(CKTcircuit *ckt, } #endif if (flag_autostop) - fprintf(stdout, "\nNote: Autostop after %e s, all measurement conditions are fulfilled\n", ckt->CKTtime); + fprintf(stdout, "\nNote: Autostop after %e s, all measurement conditions are fulfilled.\n", ckt->CKTtime); /* Final return from tran*/ return(OK); From 9d841382169f90f542120eade1f10e5d5d8baecc Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Sat, 9 Sep 2023 23:35:54 +0200 Subject: [PATCH 05/13] Remove memory leak by not mallocing unused node_ids. Remove some compiler warnings. --- src/osdi/osdisetup.c | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/osdi/osdisetup.c b/src/osdi/osdisetup.c index 08bae9942..f1a868761 100644 --- a/src/osdi/osdisetup.c +++ b/src/osdi/osdisetup.c @@ -478,8 +478,10 @@ int OSDIbindCSC(GENmodel *inModel, CKTcircuit *ckt) { GENmodel *gen_model; GENinstance *gen_inst; - /* setup a temporary buffer */ - uint32_t *node_ids = TMALLOC(uint32_t, descr->num_nodes); + NG_IGNORE(ckt); + + /* setup a temporary buffer + uint32_t *node_ids = TMALLOC(uint32_t, descr->num_nodes);*/ for (gen_model = inModel; gen_model; gen_model = gen_model->GENnextModel) { void *model = osdi_model_data(gen_model); @@ -503,8 +505,10 @@ int OSDIupdateCSC(GENmodel *inModel, CKTcircuit *ckt, bool complex) { GENmodel *gen_model; GENinstance *gen_inst; - /* setup a temporary buffer */ - uint32_t *node_ids = TMALLOC(uint32_t, descr->num_nodes); + NG_IGNORE(ckt); + + /* setup a temporary buffer + uint32_t *node_ids = TMALLOC(uint32_t, descr->num_nodes);*/ for (gen_model = inModel; gen_model; gen_model = gen_model->GENnextModel) { void *model = osdi_model_data(gen_model); @@ -522,10 +526,10 @@ int OSDIupdateCSC(GENmodel *inModel, CKTcircuit *ckt, bool complex) { return (OK); } int OSDIbindCSCComplexToReal(GENmodel *inModel, CKTcircuit *ckt) { - OSDIupdateCSC(inModel, ckt, false); + return OSDIupdateCSC(inModel, ckt, false); } int OSDIbindCSCComplex(GENmodel *inModel, CKTcircuit *ckt) { - OSDIupdateCSC(inModel, ckt, true); + return OSDIupdateCSC(inModel, ckt, true); } #endif From f2247a3c6fc2b1d9179b67b8fad714be7a7d5b10 Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Sun, 10 Sep 2023 13:43:09 +0200 Subject: [PATCH 06/13] Use sparse, as KLU will fail. --- examples/p-to-n-examples/op-test-adi.cir | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/p-to-n-examples/op-test-adi.cir b/examples/p-to-n-examples/op-test-adi.cir index a82eca9aa..cdde2ccf6 100644 --- a/examples/p-to-n-examples/op-test-adi.cir +++ b/examples/p-to-n-examples/op-test-adi.cir @@ -30,6 +30,7 @@ vin in 0 DC 0 PULSE(0.1 0.2 200uS 200uS 200uS 5m 10m) .tran 10u 10m .control +option sparse ; KLU will fail run plot dc1.v(outo) vs dc1.v(in) plot v(in) v(a1) v(outo) From 64307ba907e438ce879407f40e48b0b2cda4efdf Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Sun, 10 Sep 2023 13:44:18 +0200 Subject: [PATCH 07/13] Add a title line to 'listing r', so the resulting output may be re-loaded by the 'source' command. --- src/frontend/inp.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/frontend/inp.c b/src/frontend/inp.c index c4d6844d5..9c5f0f405 100644 --- a/src/frontend/inp.c +++ b/src/frontend/inp.c @@ -179,6 +179,8 @@ com_listing(wordlist *wl) } else { if (type != LS_DECK && type != LS_RUNNABLE) fprintf(cp_out, "\t%s\n\n", ft_curckt->ci_name); + else if (type == LS_RUNNABLE) + fprintf(cp_out, "* expanded deck of %s\n", ft_curckt->ci_name); inp_list(cp_out, expand ? ft_curckt->ci_deck : ft_curckt->ci_origdeck, ft_curckt->ci_options, type); From a6d5ce8ea60054393a54d42786efffa6cae1680f Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Sun, 10 Sep 2023 14:41:16 +0200 Subject: [PATCH 08/13] Some KLU warnings are useless for the normal user, as she or he does not have any means to further analyze or repair the issue: Warning: KLU ReFactor failed. Factoring again... Warning (ReFactor Complex): KLU Matrix is SINGULAR Numerical Rank: %d\n Singular Node: %d\n So print these messages only in debug mode. --- src/maths/KLU/klusmp.c | 58 +++++++++++++++++++++++++++--------------- src/maths/ni/niiter.c | 4 +-- 2 files changed, 39 insertions(+), 23 deletions(-) diff --git a/src/maths/KLU/klusmp.c b/src/maths/KLU/klusmp.c index 8985d440c..174dee419 100644 --- a/src/maths/KLU/klusmp.c +++ b/src/maths/KLU/klusmp.c @@ -10,6 +10,7 @@ #include "ngspice/spmatrix.h" #include "../sparse/spdefs.h" #include "ngspice/smpdefs.h" +#include "ngspice/fteext.h" #if defined (_MSC_VER) extern double scalbn(double, int); @@ -535,9 +536,11 @@ SMPcLUfac (SMPmatrix *Matrix, double PivTol) if (ret == 0) { if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) { - fprintf (stderr, "Warning (ReFactor Complex): KLU Matrix is SINGULAR\n") ; - fprintf (stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank) ; - fprintf (stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1) ; + if (ft_ngdebug) { + fprintf(stderr, "Warning (ReFactor Complex): KLU Matrix is SINGULAR\n"); + fprintf(stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank); + fprintf(stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1); + } return E_SINGULAR ; } if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) { @@ -590,9 +593,11 @@ SMPluFac (SMPmatrix *Matrix, double PivTol, double Gmin) if (ret == 0) { if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) { - fprintf (stderr, "Warning (ReFactor): KLU Matrix is SINGULAR\n") ; - fprintf (stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank) ; - fprintf (stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1) ; + if (ft_ngdebug) { + fprintf(stderr, "Warning (ReFactor): KLU Matrix is SINGULAR\n"); + fprintf(stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank); + fprintf(stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1); + } return E_SINGULAR ; } if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) { @@ -655,9 +660,11 @@ SMPluFacKLUforCIDER (SMPmatrix *Matrix) if (ret == 0) { if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) { - fprintf (stderr, "Warning (ReFactor for CIDER): KLU Matrix is SINGULAR\n") ; - fprintf (stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank) ; - fprintf (stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1) ; + if (ft_ngdebug) { + fprintf(stderr, "Warning (ReFactor for CIDER): KLU Matrix is SINGULAR\n"); + fprintf(stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank); + fprintf(stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1); + } return E_SINGULAR ; } if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) { @@ -707,9 +714,11 @@ SMPcReorder (SMPmatrix *Matrix, double PivTol, double PivRel, int *NumSwaps) if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL) { if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) { - fprintf (stderr, "Warning (Factor Complex): KLU Matrix is SINGULAR\n") ; - fprintf (stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank) ; - fprintf (stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1) ; + if (ft_ngdebug) { + fprintf(stderr, "Warning (Factor Complex): KLU Matrix is SINGULAR\n"); + fprintf(stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank); + fprintf(stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1); + } return E_SINGULAR ; } if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) { @@ -765,9 +774,11 @@ SMPreorder (SMPmatrix *Matrix, double PivTol, double PivRel, double Gmin) if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL) { if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) { - fprintf (stderr, "Warning (Factor): KLU Matrix is SINGULAR\n") ; - fprintf (stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank) ; - fprintf (stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1) ; + if (ft_ngdebug) { + fprintf(stderr, "Warning (Factor): KLU Matrix is SINGULAR\n"); + fprintf(stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank); + fprintf(stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1); + } return E_SINGULAR ; } if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) { @@ -835,9 +846,11 @@ SMPreorderKLUforCIDER (SMPmatrix *Matrix) if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL) { if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) { - fprintf (stderr, "Warning (Factor for CIDER): KLU Matrix is SINGULAR\n") ; - fprintf (stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank) ; - fprintf (stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1) ; + if (ft_ngdebug) { + fprintf(stderr, "Warning (Factor for CIDER): KLU Matrix is SINGULAR\n"); + fprintf(stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank); + fprintf(stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1); + } return E_SINGULAR ; } if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) { @@ -960,9 +973,12 @@ SMPsolve (SMPmatrix *Matrix, double RHS[], double Spare[]) if (ret == 0) { if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) { - fprintf (stderr, "Warning (Solve): KLU Matrix is SINGULAR\n") ; - fprintf (stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank) ; - fprintf (stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1) ; + if (ft_ngdebug) { + fprintf(stderr, "Warning (Solve): KLU Matrix is SINGULAR\n"); + fprintf(stderr, " Numerical Rank: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->numerical_rank); + fprintf(stderr, " Singular Node: %d\n", Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1); + } + /* FIXME: Do we need a 'return E_SINGULAR' here? */ } if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) { fprintf (stderr, "Error (Solve): KLUcommon object is NULL. A problem occurred\n") ; diff --git a/src/maths/ni/niiter.c b/src/maths/ni/niiter.c index bf945f3f6..941859b4e 100644 --- a/src/maths/ni/niiter.c +++ b/src/maths/ni/niiter.c @@ -166,8 +166,8 @@ NIiter(CKTcircuit *ckt, int maxIter) * This is my mod with KLU. It saves run-time, but also the system at the next iteration may be different. * How do we guarantee that the system is the same at the next iteration? So, the original SPARSE version below sounds like a bug. */ - - fprintf (stderr, "Warning: KLU ReFactor failed. Factoring again...\n") ; + if (ft_ngdebug) + fprintf (stderr, "Warning: KLU ReFactor failed. Factoring again...\n") ; ckt->CKTniState |= NISHOULDREORDER; ckt->CKTmatrix->SMPkluMatrix->KLUloadDiagGmin = 0 ; error = SMPreorder(ckt->CKTmatrix, ckt->CKTpivotAbsTol, ckt->CKTpivotRelTol, ckt->CKTdiagGmin); From 64c29e667ec2557811d585fbd86fd01c0bb980dc Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Mon, 11 Sep 2023 12:02:59 +0200 Subject: [PATCH 09/13] Improve error message on unimplemented dot command --- src/spicelib/parser/inp2dot.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/spicelib/parser/inp2dot.c b/src/spicelib/parser/inp2dot.c index bf035e52d..f1f5ee510 100644 --- a/src/spicelib/parser/inp2dot.c +++ b/src/spicelib/parser/inp2dot.c @@ -877,7 +877,9 @@ INP2dot(CKTcircuit *ckt, INPtables *tab, struct card *current, TSKtask *task, CK rtn = 0; goto quit; } - LITERR(" unimplemented control card - error \n"); + char *token2 = tprintf(" unimplemented dot command '%s'\n", token); + LITERR(token2); + tfree(token2); quit: tfree(token); return rtn; From dfeb0bdb4e221362964c255ebbe2035442917bf5 Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Mon, 11 Sep 2023 14:28:20 +0200 Subject: [PATCH 10/13] Improve error message on obsolete dot command --- src/spicelib/parser/inp2dot.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/spicelib/parser/inp2dot.c b/src/spicelib/parser/inp2dot.c index f1f5ee510..e1e545a2e 100644 --- a/src/spicelib/parser/inp2dot.c +++ b/src/spicelib/parser/inp2dot.c @@ -770,7 +770,9 @@ INP2dot(CKTcircuit *ckt, INPtables *tab, struct card *current, TSKtask *task, CK } else if ((strcmp(token, ".width") == 0) || strcmp(token, ".print") == 0 || strcmp(token, ".plot") == 0) { /* obsolete - ignore */ - LITERR(" Warning: obsolete control card - ignored \n"); + char* token2 = tprintf(" obsolete dot command '%s' - ignored \n", token); + LITERR(token2); + tfree(token2); goto quit; } else if ((strcmp(token, ".temp") == 0)) { /* .temp temp1 temp2 temp3 temp4 ..... */ From b7993bb53019ff742e7bcb74866e0a58085a4781 Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Mon, 11 Sep 2023 14:28:49 +0200 Subject: [PATCH 11/13] Formatting cktpzstr.c --- src/spicelib/analysis/cktpzstr.c | 1475 +++++++++++++++--------------- 1 file changed, 762 insertions(+), 713 deletions(-) diff --git a/src/spicelib/analysis/cktpzstr.c b/src/spicelib/analysis/cktpzstr.c index a18097278..ae784ef21 100644 --- a/src/spicelib/analysis/cktpzstr.c +++ b/src/spicelib/analysis/cktpzstr.c @@ -23,21 +23,21 @@ static unsigned int Debug = 3; # endif #endif -/* static function definitions */ + /* static function definitions */ -static int CKTpzStrat(PZtrial **set); -static int CKTpzStep(int strat, PZtrial **set); -static int CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set); -static int CKTpzVerify(PZtrial **set, PZtrial *new_trial); +static int CKTpzStrat(PZtrial** set); +static int CKTpzStep(int strat, PZtrial** set); +static int CKTpzRunTrial(CKTcircuit* ckt, PZtrial** new_trialp, PZtrial** set); +static int CKTpzVerify(PZtrial** set, PZtrial* new_trial); static void clear_trials(int mode); -static void check_flat(PZtrial *a, PZtrial *b); -void CKTpzUpdateSet(PZtrial **set, PZtrial *new); -void zaddeq(double *a, int *amag, double x, int xmag, double y, int ymag); -void CKTpzReset(PZtrial **set); +static void check_flat(PZtrial* a, PZtrial* b); +void CKTpzUpdateSet(PZtrial** set, PZtrial* new); +void zaddeq(double* a, int* amag, double x, int xmag, double y, int ymag); +void CKTpzReset(PZtrial** set); #ifdef PZDEBUG -static void show_trial(PZtrial *new_trial, char x); +static void show_trial(PZtrial* new_trial, char x); #endif #define NITER_LIM 200 @@ -67,25 +67,25 @@ static void show_trial(PZtrial *new_trial, char x); #define MID_RIGHT 9 #ifdef PZDEBUG -static char *snames[ ] = { - "none", - "none", - "shift left", - "shift right", - "skip left", - "skip right", - "init", - "guess", - "split left", - "split right", - "Muller", - "sym 1", - "sym 2", - "complex_init", - "complex_guess", - "quit", - "none" - }; +static char* snames[] = { + "none", + "none", + "shift left", + "shift right", + "skip left", + "skip right", + "init", + "guess", + "split left", + "split right", + "Muller", + "sym 1", + "sym 2", + "complex_init", + "complex_guess", + "quit", + "none" +}; #endif #define sgn(X) ((X) < 0 ? -1 : (X) == 0 ? 0 : 1) @@ -101,7 +101,7 @@ extern int NIpzK_mag; int CKTpzTrapped; static int NZeros, NFlat, Max_Zeros; -static PZtrial *ZeroTrial, *Trials; +static PZtrial* ZeroTrial, * Trials; static int Seq_Num; static double Guess_Param; static double High_Guess, Low_Guess; @@ -109,15 +109,15 @@ static int Last_Move, Consec_Moves; static int NIter, NTrials; static int Aberr_Num; -int PZeval(int strat, PZtrial **set, PZtrial **new_trial_p); -static PZtrial *pzseek(PZtrial *t, int dir); -static int alter(PZtrial *new, PZtrial *nearto, double abstol, double reltol); +int PZeval(int strat, PZtrial** set, PZtrial** new_trial_p); +static PZtrial* pzseek(PZtrial* t, int dir); +static int alter(PZtrial* new, PZtrial* nearto, double abstol, double reltol); int -CKTpzFindZeros(CKTcircuit *ckt, PZtrial **rootinfo, int *rootcount) +CKTpzFindZeros(CKTcircuit* ckt, PZtrial** rootinfo, int* rootcount) { - PZtrial *new_trial; - PZtrial *neighborhood[3]; + PZtrial* new_trial; + PZtrial* neighborhood[3]; int strat; int error; @@ -143,82 +143,87 @@ CKTpzFindZeros(CKTcircuit *ckt, PZtrial **rootinfo, int *rootcount) do { - while ((strat = CKTpzStrat(neighborhood)) < GUESS && !CKTpzTrapped) - if (!CKTpzStep(strat, neighborhood)) { - strat = GUESS; + while ((strat = CKTpzStrat(neighborhood)) < GUESS && !CKTpzTrapped) + if (!CKTpzStep(strat, neighborhood)) { + strat = GUESS; #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, "\t\tGuess\n"); + DEBUG(1) fprintf(stderr, "\t\tGuess\n"); #endif - break; - } + break; + } - NIter += 1; - - /* Evaluate current strategy */ - error = PZeval(strat, neighborhood, &new_trial); - if (error != OK) - return error; + NIter += 1; - error = CKTpzRunTrial(ckt, &new_trial, neighborhood); - if (error != OK) - return error; + /* Evaluate current strategy */ + error = PZeval(strat, neighborhood, &new_trial); + if (error != OK) + return error; - if (new_trial->flags & ISAROOT) { - if (CKTpzVerify(neighborhood, new_trial)) { - NIter = 0; - CKTpzReset(neighborhood); - } else - /* XXX Verify fails ?!? */ - CKTpzUpdateSet(neighborhood, new_trial); - } else if (new_trial->flags & ISANABERRATION) { - CKTpzReset(neighborhood); - Aberr_Num += 1; - tfree(new_trial); - } else if (new_trial->flags & ISAMINIMA) { - neighborhood[0] = NULL; - neighborhood[1] = new_trial; - neighborhood[2] = NULL; - } else { - CKTpzUpdateSet(neighborhood, new_trial); /* Replace a value */ - } + error = CKTpzRunTrial(ckt, &new_trial, neighborhood); + if (error != OK) + return error; - if (SPfrontEnd->IFpauseTest()) { - SPfrontEnd->IFerrorf (ERR_WARNING, "Pole-Zero analysis interrupted; %d trials, %d roots\n", Seq_Num, NZeros); - error = E_PAUSE; - break; - } + if (new_trial->flags & ISAROOT) { + if (CKTpzVerify(neighborhood, new_trial)) { + NIter = 0; + CKTpzReset(neighborhood); + } + else + /* XXX Verify fails ?!? */ + CKTpzUpdateSet(neighborhood, new_trial); + } + else if (new_trial->flags & ISANABERRATION) { + CKTpzReset(neighborhood); + Aberr_Num += 1; + tfree(new_trial); + } + else if (new_trial->flags & ISAMINIMA) { + neighborhood[0] = NULL; + neighborhood[1] = new_trial; + neighborhood[2] = NULL; + } + else { + CKTpzUpdateSet(neighborhood, new_trial); /* Replace a value */ + } + + if (SPfrontEnd->IFpauseTest()) { + SPfrontEnd->IFerrorf(ERR_WARNING, "Pole-Zero analysis interrupted; %d trials, %d roots\n", Seq_Num, NZeros); + error = E_PAUSE; + break; + } } while (High_Guess - Low_Guess < 1e40 - && NZeros < Max_Zeros - && NIter < NITER_LIM && Aberr_Num < 3 - && High_Guess - Low_Guess < 1e35 /* XXX Should use mach const */ - && (!neighborhood[0] || !neighborhood[2] || CKTpzTrapped - || neighborhood[2]->s.real - neighborhood[0]->s.real < 1e22)); - /* XXX ZZZ */ + && NZeros < Max_Zeros + && NIter < NITER_LIM && Aberr_Num < 3 + && High_Guess - Low_Guess < 1e35 /* XXX Should use mach const */ + && (!neighborhood[0] || !neighborhood[2] || CKTpzTrapped + || neighborhood[2]->s.real - neighborhood[0]->s.real < 1e22)); + /* XXX ZZZ */ #ifdef PZDEBUG DEBUG(1) fprintf(stderr, - "Finished: NFlat %d, NZeros: %d, NTrials %d, Guess %g to %g, aber %d\n", - NFlat, NZeros, NTrials, Low_Guess, High_Guess, Aberr_Num); + "Finished: NFlat %d, NZeros: %d, NTrials %d, Guess %g to %g, aber %d\n", + NFlat, NZeros, NTrials, Low_Guess, High_Guess, Aberr_Num); #endif if (NZeros >= Seq_Num - 1) { - /* Short */ - clear_trials(ISAROOT); - *rootinfo = NULL; - *rootcount = 0; - MERROR(E_SHORT, "The input signal is shorted on the way to the output"); - } else - clear_trials(0); + /* Short */ + clear_trials(ISAROOT); + *rootinfo = NULL; + *rootcount = 0; + MERROR(E_SHORT, "The input signal is shorted on the way to the output"); + } + else + clear_trials(0); *rootinfo = Trials; *rootcount = NZeros; if (Aberr_Num > 2) { - SPfrontEnd->IFerrorf (ERR_WARNING, "Pole-zero converging to numerical aberrations; giving up after %d trials", Seq_Num); + SPfrontEnd->IFerrorf(ERR_WARNING, "Pole-zero converging to numerical aberrations; giving up after %d trials", Seq_Num); } if (NIter >= NITER_LIM) { - SPfrontEnd->IFerrorf (ERR_WARNING, "Pole-zero iteration limit reached; giving up after %d trials", Seq_Num); + SPfrontEnd->IFerrorf(ERR_WARNING, "Pole-zero iteration limit reached; giving up after %d trials", Seq_Num); } return error; @@ -227,12 +232,12 @@ CKTpzFindZeros(CKTcircuit *ckt, PZtrial **rootinfo, int *rootcount) /* PZeval: evaluate an estimation function (given by 'strat') for the next guess (returned in a PZtrial) */ -/* XXX ZZZ */ + /* XXX ZZZ */ int -PZeval(int strat, PZtrial **set, PZtrial **new_trial_p) +PZeval(int strat, PZtrial** set, PZtrial** new_trial_p) { int error; - PZtrial *new_trial; + PZtrial* new_trial; new_trial = TMALLOC(PZtrial, 1); new_trial->multiplicity = 0; @@ -241,153 +246,158 @@ PZeval(int strat, PZtrial **set, PZtrial **new_trial_p) switch (strat) { case GUESS: - if (High_Guess < Low_Guess) - Guess_Param = 0.0; - else if (Guess_Param > 0.0) { - if (High_Guess > 0.0) - Guess_Param = High_Guess * 10.0; - else - Guess_Param = 1.0; - } else { - if (Low_Guess < 0.0) - Guess_Param = Low_Guess * 10.0; - else - Guess_Param = -1.0; - } - if (High_Guess < Guess_Param) - High_Guess = Guess_Param; - if (Low_Guess > Guess_Param) - Low_Guess = Guess_Param; - new_trial->s.real = Guess_Param; - if (set[1]) - new_trial->s.imag = set[1]->s.imag; - else - new_trial->s.imag = 0.0; - error = OK; - break; + if (High_Guess < Low_Guess) + Guess_Param = 0.0; + else if (Guess_Param > 0.0) { + if (High_Guess > 0.0) + Guess_Param = High_Guess * 10.0; + else + Guess_Param = 1.0; + } + else { + if (Low_Guess < 0.0) + Guess_Param = Low_Guess * 10.0; + else + Guess_Param = -1.0; + } + if (High_Guess < Guess_Param) + High_Guess = Guess_Param; + if (Low_Guess > Guess_Param) + Low_Guess = Guess_Param; + new_trial->s.real = Guess_Param; + if (set[1]) + new_trial->s.imag = set[1]->s.imag; + else + new_trial->s.imag = 0.0; + error = OK; + break; case SYM: case SYM2: - error = NIpzSym(set, new_trial); + error = NIpzSym(set, new_trial); - if (CKTpzTrapped == 1) { - if (new_trial->s.real < set[0]->s.real - || new_trial->s.real > set[1]->s.real) { + if (CKTpzTrapped == 1) { + if (new_trial->s.real < set[0]->s.real + || new_trial->s.real > set[1]->s.real) { #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, - "FIXED UP BAD Strat: %s (%d) was (%.15g,%.15g)\n", - snames[strat], CKTpzTrapped, - new_trial->s.real, new_trial->s.imag); + DEBUG(1) fprintf(stderr, + "FIXED UP BAD Strat: %s (%d) was (%.15g,%.15g)\n", + snames[strat], CKTpzTrapped, + new_trial->s.real, new_trial->s.imag); #endif - new_trial->s.real = (set[0]->s.real + set[1]->s.real) / 2.0; - } - } else if (CKTpzTrapped == 2) { - if (new_trial->s.real < set[1]->s.real - || new_trial->s.real > set[2]->s.real) { + new_trial->s.real = (set[0]->s.real + set[1]->s.real) / 2.0; + } + } + else if (CKTpzTrapped == 2) { + if (new_trial->s.real < set[1]->s.real + || new_trial->s.real > set[2]->s.real) { #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, - "FIXED UP BAD Strat: %s (%d) was (%.15g,%.15g)\n", - snames[strat], CKTpzTrapped, - new_trial->s.real, new_trial->s.imag); + DEBUG(1) fprintf(stderr, + "FIXED UP BAD Strat: %s (%d) was (%.15g,%.15g)\n", + snames[strat], CKTpzTrapped, + new_trial->s.real, new_trial->s.imag); #endif - new_trial->s.real = (set[1]->s.real + set[2]->s.real) / 2.0; - } - } else if (CKTpzTrapped == 3) { - if (new_trial->s.real <= set[0]->s.real - || (new_trial->s.real == set[1]->s.real - && new_trial->s.imag == set[1]->s.imag) - || new_trial->s.real >= set[2]->s.real) { + new_trial->s.real = (set[1]->s.real + set[2]->s.real) / 2.0; + } + } + else if (CKTpzTrapped == 3) { + if (new_trial->s.real <= set[0]->s.real + || (new_trial->s.real == set[1]->s.real + && new_trial->s.imag == set[1]->s.imag) + || new_trial->s.real >= set[2]->s.real) { #ifdef PZDEBUG - DEBUG(1) - fprintf(stderr, - "FIXED UP BAD Strat: %s (%d), was (%.15g %.15g)\n", - snames[strat], CKTpzTrapped, - new_trial->s.real, new_trial->s.imag); + DEBUG(1) + fprintf(stderr, + "FIXED UP BAD Strat: %s (%d), was (%.15g %.15g)\n", + snames[strat], CKTpzTrapped, + new_trial->s.real, new_trial->s.imag); #endif - new_trial->s.real = (set[0]->s.real + set[2]->s.real) / 2.0; - if (new_trial->s.real == set[1]->s.real) { + new_trial->s.real = (set[0]->s.real + set[2]->s.real) / 2.0; + if (new_trial->s.real == set[1]->s.real) { #ifdef PZDEBUG - DEBUG(1) - fprintf(stderr, "Still off!"); + DEBUG(1) + fprintf(stderr, "Still off!"); #endif - if (Last_Move == MID_LEFT || Last_Move == NEAR_RIGHT) - new_trial->s.real = (set[0]->s.real + set[1]->s.real) - / 2.0; - else - new_trial->s.real = (set[1]->s.real + set[2]->s.real) - / 2.0; - } - } - } + if (Last_Move == MID_LEFT || Last_Move == NEAR_RIGHT) + new_trial->s.real = (set[0]->s.real + set[1]->s.real) + / 2.0; + else + new_trial->s.real = (set[1]->s.real + set[2]->s.real) + / 2.0; + } + } + } - break; + break; case COMPLEX_INIT: - /* Not automatic */ + /* Not automatic */ #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, "\tPZ minima at: %-30g %d\n", - NIpzK, NIpzK_mag); + DEBUG(1) fprintf(stderr, "\tPZ minima at: %-30g %d\n", + NIpzK, NIpzK_mag); #endif - new_trial->s.real = set[1]->s.real; + new_trial->s.real = set[1]->s.real; - /* NIpzK is a good idea, but the value gets trashed - * due to the numerics when zooming in on a minima. - * The key is to know when to stop taking new values for NIpzK - * (which I don't). For now I take the first value indicated - * by the NIpzSym2 routine. A "hack". - */ + /* NIpzK is a good idea, but the value gets trashed + * due to the numerics when zooming in on a minima. + * The key is to know when to stop taking new values for NIpzK + * (which I don't). For now I take the first value indicated + * by the NIpzSym2 routine. A "hack". + */ - if (NIpzK != 0.0 && NIpzK_mag > -10) { - while (NIpzK_mag > 0) { - NIpzK *= 2.0; - NIpzK_mag -= 1; - } - while (NIpzK_mag < 0) { - NIpzK /= 2.0; - NIpzK_mag += 1; - } - new_trial->s.imag = NIpzK; - } else - new_trial->s.imag = 10000.0; + if (NIpzK != 0.0 && NIpzK_mag > -10) { + while (NIpzK_mag > 0) { + NIpzK *= 2.0; + NIpzK_mag -= 1; + } + while (NIpzK_mag < 0) { + NIpzK /= 2.0; + NIpzK_mag += 1; + } + new_trial->s.imag = NIpzK; + } + else + new_trial->s.imag = 10000.0; - /* - * Reset NIpzK so the same value doesn't get used again. - */ + /* + * Reset NIpzK so the same value doesn't get used again. + */ - NIpzK = 0.0; - NIpzK_mag = 0; - error = OK; - break; + NIpzK = 0.0; + NIpzK_mag = 0; + error = OK; + break; case COMPLEX_GUESS: - if (!set[2]) { - new_trial->s.real = set[0]->s.real; - new_trial->s.imag = 1.0e8; - } else { - new_trial->s.real = set[0]->s.real; - new_trial->s.imag = 1.0e12; - } - error = OK; - break; + if (!set[2]) { + new_trial->s.real = set[0]->s.real; + new_trial->s.imag = 1.0e8; + } + else { + new_trial->s.real = set[0]->s.real; + new_trial->s.imag = 1.0e12; + } + error = OK; + break; case MULLER: - error = NIpzMuller(set, new_trial); - break; + error = NIpzMuller(set, new_trial); + break; case SPLIT_LEFT: - new_trial->s.real = (set[0]->s.real + 2 * set[1]->s.real) / 3.0; - error = OK; - break; + new_trial->s.real = (set[0]->s.real + 2 * set[1]->s.real) / 3.0; + error = OK; + break; case SPLIT_RIGHT: - new_trial->s.real = (set[2]->s.real + 2 * set[1]->s.real) / 3.0; - error = OK; - break; + new_trial->s.real = (set[2]->s.real + 2 * set[1]->s.real) / 3.0; + error = OK; + break; default: - MERROR(E_PANIC, "Step type unknown"); - break; + MERROR(E_PANIC, "Step type unknown"); + break; } *new_trial_p = new_trial; @@ -397,8 +407,8 @@ PZeval(int strat, PZtrial **set, PZtrial **new_trial_p) /* CKTpzStrat: given three points, determine a good direction or method for guessing the next zero */ -/* XXX ZZZ what is a strategy for complex hunting? */ -int CKTpzStrat(PZtrial **set) + /* XXX ZZZ what is a strategy for complex hunting? */ +int CKTpzStrat(PZtrial** set) { int suggestion; double a, b; @@ -409,95 +419,105 @@ int CKTpzStrat(PZtrial **set) new_trap = 0; if (set[1] && (set[1]->flags & ISAMINIMA)) { - suggestion = COMPLEX_INIT; - } else if (set[0] && set[0]->s.imag != 0.0) { - if (!set[1] || !set[2]) - suggestion = COMPLEX_GUESS; - else - suggestion = MULLER; - } else if (!set[0] || !set[1] || !set[2]) { - suggestion = INIT; - } else { - if (sgn(set[0]->f_def.real) != sgn(set[1]->f_def.real)) { - /* Zero crossing between s[0] and s[1] */ - new_trap = 1; - suggestion = SYM2; - } else if (sgn(set[1]->f_def.real) != sgn(set[2]->f_def.real)) { - /* Zero crossing between s[1] and s[2] */ - new_trap = 2; - suggestion = SYM2; - } else { + suggestion = COMPLEX_INIT; + } + else if (set[0] && set[0]->s.imag != 0.0) { + if (!set[1] || !set[2]) + suggestion = COMPLEX_GUESS; + else + suggestion = MULLER; + } + else if (!set[0] || !set[1] || !set[2]) { + suggestion = INIT; + } + else { + if (sgn(set[0]->f_def.real) != sgn(set[1]->f_def.real)) { + /* Zero crossing between s[0] and s[1] */ + new_trap = 1; + suggestion = SYM2; + } + else if (sgn(set[1]->f_def.real) != sgn(set[2]->f_def.real)) { + /* Zero crossing between s[1] and s[2] */ + new_trap = 2; + suggestion = SYM2; + } + else { - zaddeq(&a, &a_mag, set[1]->f_def.real, set[1]->mag_def, - -set[0]->f_def.real, set[0]->mag_def); - zaddeq(&b, &b_mag, set[2]->f_def.real, set[2]->mag_def, - -set[1]->f_def.real, set[1]->mag_def); + zaddeq(&a, &a_mag, set[1]->f_def.real, set[1]->mag_def, + -set[0]->f_def.real, set[0]->mag_def); + zaddeq(&b, &b_mag, set[2]->f_def.real, set[2]->mag_def, + -set[1]->f_def.real, set[1]->mag_def); - if (!CKTpzTrapped) { + if (!CKTpzTrapped) { - k1 = set[1]->s.real - set[0]->s.real; - k2 = set[2]->s.real - set[1]->s.real; - if (a_mag + 10 < set[0]->mag_def - && a_mag + 10 < set[1]->mag_def - && b_mag + 10 < set[1]->mag_def - && b_mag + 10 < set[2]->mag_def) { - if (k1 > k2) - suggestion = SKIP_RIGHT; - else - suggestion = SKIP_LEFT; - } else if (sgn(a) != -sgn(b)) { - if (a == 0.0) - suggestion = SKIP_LEFT; - else if (b == 0.0) - suggestion = SKIP_RIGHT; - else if (sgn(a) == sgn(set[1]->f_def.real)) - suggestion = SHIFT_LEFT; - else - suggestion = SHIFT_RIGHT; - } else if (sgn(a) == -sgn(set[1]->f_def.real)) { - new_trap = 3; - /* minima in magnitude above the x axis */ - /* Search for exact mag. minima, look for complex pair */ - suggestion = SYM; - } else if (k1 > k2) - suggestion = SKIP_RIGHT; - else - suggestion = SKIP_LEFT; - } else { - new_trap = 3; /* still */ - /* XXX ? Are these tests needed or is SYM safe all the time? */ - if (sgn(a) != sgn(b)) { - /* minima in magnitude */ - /* Search for exact mag. minima, look for complex pair */ - suggestion = SYM; - } else if (a_mag > b_mag || (a_mag == b_mag - && fabs(a) > fabs(b))) - suggestion = SPLIT_LEFT; - else - suggestion = SPLIT_RIGHT; - } - } - if (Consec_Moves >= 3 && CKTpzTrapped == new_trap) { - new_trap = CKTpzTrapped; - if (Last_Move == MID_LEFT || Last_Move == NEAR_RIGHT) - suggestion = SPLIT_LEFT; - else if (Last_Move == MID_RIGHT || Last_Move == NEAR_LEFT) - suggestion = SPLIT_RIGHT; - else - abort( ); /* XXX */ - Consec_Moves = 0; - } + k1 = set[1]->s.real - set[0]->s.real; + k2 = set[2]->s.real - set[1]->s.real; + if (a_mag + 10 < set[0]->mag_def + && a_mag + 10 < set[1]->mag_def + && b_mag + 10 < set[1]->mag_def + && b_mag + 10 < set[2]->mag_def) { + if (k1 > k2) + suggestion = SKIP_RIGHT; + else + suggestion = SKIP_LEFT; + } + else if (sgn(a) != -sgn(b)) { + if (a == 0.0) + suggestion = SKIP_LEFT; + else if (b == 0.0) + suggestion = SKIP_RIGHT; + else if (sgn(a) == sgn(set[1]->f_def.real)) + suggestion = SHIFT_LEFT; + else + suggestion = SHIFT_RIGHT; + } + else if (sgn(a) == -sgn(set[1]->f_def.real)) { + new_trap = 3; + /* minima in magnitude above the x axis */ + /* Search for exact mag. minima, look for complex pair */ + suggestion = SYM; + } + else if (k1 > k2) + suggestion = SKIP_RIGHT; + else + suggestion = SKIP_LEFT; + } + else { + new_trap = 3; /* still */ + /* XXX ? Are these tests needed or is SYM safe all the time? */ + if (sgn(a) != sgn(b)) { + /* minima in magnitude */ + /* Search for exact mag. minima, look for complex pair */ + suggestion = SYM; + } + else if (a_mag > b_mag || (a_mag == b_mag + && fabs(a) > fabs(b))) + suggestion = SPLIT_LEFT; + else + suggestion = SPLIT_RIGHT; + } + } + if (Consec_Moves >= 3 && CKTpzTrapped == new_trap) { + new_trap = CKTpzTrapped; + if (Last_Move == MID_LEFT || Last_Move == NEAR_RIGHT) + suggestion = SPLIT_LEFT; + else if (Last_Move == MID_RIGHT || Last_Move == NEAR_LEFT) + suggestion = SPLIT_RIGHT; + else + abort(); /* XXX */ + Consec_Moves = 0; + } } CKTpzTrapped = new_trap; #ifdef PZDEBUG DEBUG(1) { - if (set[0] && set[1] && set[2]) - fprintf(stderr, "given %.15g %.15g / %.15g %.15g / %.15g %.15g\n", - set[0]->s.real, set[0]->s.imag, set[1]->s.real, set[1]->s.imag, - set[2]->s.real, set[2]->s.imag); - fprintf(stderr, "suggestion(%d/%d/%d | %d): %s\n", - NFlat, NZeros, Max_Zeros, CKTpzTrapped, snames[suggestion]); + if (set[0] && set[1] && set[2]) + fprintf(stderr, "given %.15g %.15g / %.15g %.15g / %.15g %.15g\n", + set[0]->s.real, set[0]->s.imag, set[1]->s.real, set[1]->s.imag, + set[2]->s.real, set[2]->s.imag); + fprintf(stderr, "suggestion(%d/%d/%d | %d): %s\n", + NFlat, NZeros, Max_Zeros, CKTpzTrapped, snames[suggestion]); } #endif return suggestion; @@ -505,13 +525,13 @@ int CKTpzStrat(PZtrial **set) /* CKTpzRunTrial: eval the function at a given 's', fold in deflation */ -int -CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set) +int +CKTpzRunTrial(CKTcircuit* ckt, PZtrial** new_trialp, PZtrial** set) { - PZAN *job = (PZAN *) ckt->CKTcurJob; + PZAN* job = (PZAN*)ckt->CKTcurJob; - PZtrial *match, *new_trial; - PZtrial *p, *prev; + PZtrial* match, * new_trial; + PZtrial* p, * prev; SPcomplex def_frac, diff_frac; double reltol, abstol; int def_mag, diff_mag, error = 0; @@ -522,10 +542,10 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set) new_trial = *new_trialp; if (new_trial->s.imag < 0.0) - new_trial->s.imag *= -1.0; + new_trial->s.imag *= -1.0; /* Insert the trial into the list of Trials, while calculating - the deflation factor from previous zeros */ + the deflation factor from previous zeros */ pretest = 0; shifted = 0; @@ -533,222 +553,231 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set) do { - def_mag = 0; - def_frac.real = 1.0; - def_frac.imag = 0.0; - was_shifted = shifted; - shifted = 0; + def_mag = 0; + def_frac.real = 1.0; + def_frac.imag = 0.0; + was_shifted = shifted; + shifted = 0; - prev = NULL; - match = NULL; + prev = NULL; + match = NULL; - for (p = Trials; p != NULL; p = p->next) { + for (p = Trials; p != NULL; p = p->next) { - C_SUBEQ(diff_frac,p->s,new_trial->s); + C_SUBEQ(diff_frac, p->s, new_trial->s); - if (diff_frac.real < 0.0 - || (diff_frac.real == 0.0 && diff_frac.imag < 0.0)) { - prev = p; - } + if (diff_frac.real < 0.0 + || (diff_frac.real == 0.0 && diff_frac.imag < 0.0)) { + prev = p; + } - if (p->flags & ISAROOT) { - abstol = 1e-5; - reltol = 1e-6; - } else { - abstol = 1e-20; - reltol = 1e-12; - } + if (p->flags & ISAROOT) { + abstol = 1e-5; + reltol = 1e-6; + } + else { + abstol = 1e-20; + reltol = 1e-12; + } - if (diff_frac.imag == 0.0 && - fabs(diff_frac.real) / (fabs(p->s.real) + abstol/reltol) - < reltol) { + if (diff_frac.imag == 0.0 && + fabs(diff_frac.real) / (fabs(p->s.real) + abstol / reltol) + < reltol) { #ifdef PZDEBUG - DEBUG(1) { - fprintf(stderr, - "diff_frac.real = %10g, p->s = %10g, nt = %10g\n", - diff_frac.real, p->s.real, new_trial->s.real); - fprintf(stderr, "ab=%g,rel=%g\n", abstol, reltol); - } + DEBUG(1) { + fprintf(stderr, + "diff_frac.real = %10g, p->s = %10g, nt = %10g\n", + diff_frac.real, p->s.real, new_trial->s.real); + fprintf(stderr, "ab=%g,rel=%g\n", abstol, reltol); + } #endif - if (was_shifted || p->count >= 3 - || !alter(new_trial, set[1], abstol, reltol)) { - /* assume either a root or minima */ - p->count = 0; - pretest = 1; - break; - } else - p->count += 1; /* try to shift */ + if (was_shifted || p->count >= 3 + || !alter(new_trial, set[1], abstol, reltol)) { + /* assume either a root or minima */ + p->count = 0; + pretest = 1; + break; + } + else + p->count += 1; /* try to shift */ - shifted = 1; /* Re-calculate deflation */ - break; + shifted = 1; /* Re-calculate deflation */ + break; - } else { - if (!CKTpzTrapped) - p->count = 0; - if (p->flags & ISAROOT) { - diff_mag = 0; - C_NORM(diff_frac,diff_mag); - if (diff_frac.imag != 0.0) { - C_MAG2(diff_frac); - diff_mag *= 2; - } - C_NORM(diff_frac,diff_mag); + } + else { + if (!CKTpzTrapped) + p->count = 0; + if (p->flags & ISAROOT) { + diff_mag = 0; + C_NORM(diff_frac, diff_mag); + if (diff_frac.imag != 0.0) { + C_MAG2(diff_frac); + diff_mag *= 2; + } + C_NORM(diff_frac, diff_mag); - for (i = p->multiplicity; i > 0; i--) { - C_MUL(def_frac,diff_frac); - def_mag += diff_mag; - C_NORM(def_frac,def_mag); - } - } else if (!match) - match = p; - } - } + for (i = p->multiplicity; i > 0; i--) { + C_MUL(def_frac, diff_frac); + def_mag += diff_mag; + C_NORM(def_frac, def_mag); + } + } + else if (!match) + match = p; + } + } } while (shifted); if (pretest) { #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, "Pre-test taken\n"); + DEBUG(1) fprintf(stderr, "Pre-test taken\n"); - /* XXX Should catch the double-zero right off - * if K is 0.0 - * instead of forcing a re-converge */ - DEBUG(1) { - fprintf(stderr, "NIpzK == %g, mag = %d\n", NIpzK, NIpzK_mag); - fprintf(stderr, "over at %.30g %.30g (new %.30g %.30g, %x)\n", - p->s.real, p->s.imag, new_trial->s.real, new_trial->s.imag, - p->flags); - } + /* XXX Should catch the double-zero right off + * if K is 0.0 + * instead of forcing a re-converge */ + DEBUG(1) { + fprintf(stderr, "NIpzK == %g, mag = %d\n", NIpzK, NIpzK_mag); + fprintf(stderr, "over at %.30g %.30g (new %.30g %.30g, %x)\n", + p->s.real, p->s.imag, new_trial->s.real, new_trial->s.imag, + p->flags); + } #endif - if (!(p->flags & ISAROOT) && CKTpzTrapped == 3 - && NIpzK != 0.0 && NIpzK_mag > -10) { + if (!(p->flags & ISAROOT) && CKTpzTrapped == 3 + && NIpzK != 0.0 && NIpzK_mag > -10) { #ifdef notdef -// if (p->flags & ISAROOT) { -// /* Ugh! muller doesn't work right */ -// new_trial->flags = ISAMINIMA; -// new_trial->s.imag = scalb(NIpzK, (int) (NIpzK_mag / 2)); -// pretest = 0; -// } else { + // if (p->flags & ISAROOT) { + // /* Ugh! muller doesn't work right */ + // new_trial->flags = ISAMINIMA; + // new_trial->s.imag = scalb(NIpzK, (int) (NIpzK_mag / 2)); + // pretest = 0; + // } else { #endif - p->flags |= ISAMINIMA; - tfree(new_trial); - *new_trialp = p; - repeat = 1; - } else if (p->flags & ISAROOT) { + p->flags |= ISAMINIMA; + tfree(new_trial); + *new_trialp = p; + repeat = 1; + } + else if (p->flags & ISAROOT) { #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, "Repeat at %.30g %.30g\n", - p->s.real, p->s.imag); + DEBUG(1) fprintf(stderr, "Repeat at %.30g %.30g\n", + p->s.real, p->s.imag); #endif - *new_trialp = p; - p->flags |= ISAREPEAT; - p->multiplicity += 1; - repeat = 1; - } else { - /* Regular zero, as precise as we can get it */ - error = E_SINGULAR; - } + * new_trialp = p; + p->flags |= ISAREPEAT; + p->multiplicity += 1; + repeat = 1; + } + else { + /* Regular zero, as precise as we can get it */ + error = E_SINGULAR; + } } if (!repeat) { - if (!pretest) { - /* Run the trial */ - ckt->CKTniState |= NIPZSHOULDREORDER; /* XXX */ - if (!(ckt->CKTniState & NIPZSHOULDREORDER)) { - CKTpzLoad(ckt, &new_trial->s); + if (!pretest) { + /* Run the trial */ + ckt->CKTniState |= NIPZSHOULDREORDER; /* XXX */ + if (!(ckt->CKTniState & NIPZSHOULDREORDER)) { + CKTpzLoad(ckt, &new_trial->s); #ifdef PZDEBUG - DEBUG(3) { - printf("Original:\n"); - SMPprint(ckt->CKTmatrix, stdout); - } + DEBUG(3) { + printf("Original:\n"); + SMPprint(ckt->CKTmatrix, stdout); + } #endif - error = SMPcLUfac(ckt->CKTmatrix, ckt->CKTpivotAbsTol); - if (error == E_SINGULAR) { + error = SMPcLUfac(ckt->CKTmatrix, ckt->CKTpivotAbsTol); + if (error == E_SINGULAR) { #ifdef PZDEBUG - DEBUG(1) printf("Needs reordering\n"); + DEBUG(1) printf("Needs reordering\n"); #endif - ckt->CKTniState |= NIPZSHOULDREORDER; - } else if (error != OK) - return error; - } - if (ckt->CKTniState & NIPZSHOULDREORDER) { - CKTpzLoad(ckt, &new_trial->s); - error = SMPcReorder(ckt->CKTmatrix, 1.0e-30, - 0.0 /* 0.1 Piv. Rel. */, - &(job->PZnumswaps)); - } + ckt->CKTniState |= NIPZSHOULDREORDER; + } + else if (error != OK) + return error; + } + if (ckt->CKTniState & NIPZSHOULDREORDER) { + CKTpzLoad(ckt, &new_trial->s); + error = SMPcReorder(ckt->CKTmatrix, 1.0e-30, + 0.0 /* 0.1 Piv. Rel. */, + &(job->PZnumswaps)); + } - if (error != E_SINGULAR) { - ckt->CKTniState &= ~NIPZSHOULDREORDER; + if (error != E_SINGULAR) { + ckt->CKTniState &= ~NIPZSHOULDREORDER; #ifdef PZDEBUG - DEBUG(3) { - printf("Factored:\n"); - SMPprint(ckt->CKTmatrix, NULL); - } + DEBUG(3) { + printf("Factored:\n"); + SMPprint(ckt->CKTmatrix, NULL); + } #endif - error = SMPcDProd(ckt->CKTmatrix, &new_trial->f_raw, - &new_trial->mag_raw); - } - } + error = SMPcDProd(ckt->CKTmatrix, &new_trial->f_raw, + &new_trial->mag_raw); + } + } - if (error == E_SINGULAR || (new_trial->f_raw.real == 0.0 - && new_trial->f_raw.imag == 0.0)) { - new_trial->f_raw.real = 0.0; - new_trial->f_raw.imag = 0.0; - new_trial->mag_raw = 0; - new_trial->f_def.real = 0.0; - new_trial->f_def.imag = 0.0; - new_trial->mag_def = 0; - new_trial->flags = ISAROOT; - /*printf("SMP Det: Singular\n");*/ - } else if (error != OK) - return error; - else { + if (error == E_SINGULAR || (new_trial->f_raw.real == 0.0 + && new_trial->f_raw.imag == 0.0)) { + new_trial->f_raw.real = 0.0; + new_trial->f_raw.imag = 0.0; + new_trial->mag_raw = 0; + new_trial->f_def.real = 0.0; + new_trial->f_def.imag = 0.0; + new_trial->mag_def = 0; + new_trial->flags = ISAROOT; + /*printf("SMP Det: Singular\n");*/ + } + else if (error != OK) + return error; + else { + + /* PZnumswaps is either 0 or 1 */ + new_trial->f_raw.real *= job->PZnumswaps; + new_trial->f_raw.imag *= job->PZnumswaps; - /* PZnumswaps is either 0 or 1 */ - new_trial->f_raw.real *= job->PZnumswaps; - new_trial->f_raw.imag *= job->PZnumswaps; - #ifdef PZDEBUG - printf("SMP Det: (%g,%g)^%d\n", new_trial->f_raw.real, - new_trial->f_raw.imag, new_trial->mag_raw); + printf("SMP Det: (%g,%g)^%d\n", new_trial->f_raw.real, + new_trial->f_raw.imag, new_trial->mag_raw); #endif - new_trial->f_def.real = new_trial->f_raw.real; - new_trial->f_def.imag = new_trial->f_raw.imag; - new_trial->mag_def = new_trial->mag_raw; + new_trial->f_def.real = new_trial->f_raw.real; + new_trial->f_def.imag = new_trial->f_raw.imag; + new_trial->mag_def = new_trial->mag_raw; - C_DIV(new_trial->f_def,def_frac); - new_trial->mag_def -= def_mag; - C_NORM(new_trial->f_def,new_trial->mag_def); - } + C_DIV(new_trial->f_def, def_frac); + new_trial->mag_def -= def_mag; + C_NORM(new_trial->f_def, new_trial->mag_def); + } - /* Link into the rest of the list */ - if (prev) { - new_trial->next = prev->next; - if (prev->next) - prev->next->prev = new_trial; - prev->next = new_trial; - } else { - if (Trials) - Trials->prev = new_trial; - else - ZeroTrial = new_trial; - new_trial->next = Trials; - Trials = new_trial; - } - new_trial->prev = prev; + /* Link into the rest of the list */ + if (prev) { + new_trial->next = prev->next; + if (prev->next) + prev->next->prev = new_trial; + prev->next = new_trial; + } + else { + if (Trials) + Trials->prev = new_trial; + else + ZeroTrial = new_trial; + new_trial->next = Trials; + Trials = new_trial; + } + new_trial->prev = prev; - NTrials += 1; + NTrials += 1; - if (!(new_trial->flags & ISAROOT)) { - if (match) - check_flat(match, new_trial); - else - NFlat = 1; - } + if (!(new_trial->flags & ISAROOT)) { + if (match) + check_flat(match, new_trial); + else + NFlat = 1; + } } #ifdef PZDEBUG @@ -760,94 +789,95 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set) /* Process a zero; inc. zero count, deflate other trials */ -int -CKTpzVerify(PZtrial **set, PZtrial *new_trial) +int +CKTpzVerify(PZtrial** set, PZtrial* new_trial) { - PZtrial *next; + PZtrial* next; int diff_mag; SPcomplex diff_frac; double tdiff; - PZtrial *t, *prev; + PZtrial* t, * prev; NG_IGNORE(set); NZeros += 1; if (new_trial->s.imag != 0.0) - NZeros += 1; + NZeros += 1; NFlat = 0; if (new_trial->multiplicity == 0) { - new_trial->flags |= ISAROOT; - new_trial->multiplicity = 1; + new_trial->flags |= ISAROOT; + new_trial->multiplicity = 1; } prev = NULL; for (t = Trials; t; t = next) { - next = t->next; + next = t->next; - if (t->flags & ISAROOT) { - prev = t; - /* Don't need to bother */ - continue; - } + if (t->flags & ISAROOT) { + prev = t; + /* Don't need to bother */ + continue; + } - C_SUBEQ(diff_frac,new_trial->s,t->s); - if (new_trial->s.imag != 0.0) - C_MAG2(diff_frac); + C_SUBEQ(diff_frac, new_trial->s, t->s); + if (new_trial->s.imag != 0.0) + C_MAG2(diff_frac); - tdiff = diff_frac.real; - /* Note that Verify is called for each time the root is found, so - * multiplicity is not significant - */ - if (diff_frac.real != 0.0) { - diff_mag = 0; - C_NORM(diff_frac,diff_mag); - diff_mag *= -1; - C_DIV(t->f_def,diff_frac); - C_NORM(t->f_def,diff_mag); - t->mag_def += diff_mag; - } + tdiff = diff_frac.real; + /* Note that Verify is called for each time the root is found, so + * multiplicity is not significant + */ + if (diff_frac.real != 0.0) { + diff_mag = 0; + C_NORM(diff_frac, diff_mag); + diff_mag *= -1; + C_DIV(t->f_def, diff_frac); + C_NORM(t->f_def, diff_mag); + t->mag_def += diff_mag; + } - if (t->s.imag != 0.0 - || fabs(tdiff) / (fabs(new_trial->s.real) + 200) < 0.005) { - if (prev) - prev->next = t->next; - if (t->next) - t->next->prev = prev; - NTrials -= 1; + if (t->s.imag != 0.0 + || fabs(tdiff) / (fabs(new_trial->s.real) + 200) < 0.005) { + if (prev) + prev->next = t->next; + if (t->next) + t->next->prev = prev; + NTrials -= 1; #ifdef PZDEBUG - show_trial(t, '-'); + show_trial(t, '-'); #endif - if (t == ZeroTrial) { - if (t->next) - ZeroTrial = t->next; - else if (t->prev) - ZeroTrial = t->prev; - else - ZeroTrial = NULL; - } - if (t == Trials) { - Trials = t->next; - } - tfree(t); - } else { + if (t == ZeroTrial) { + if (t->next) + ZeroTrial = t->next; + else if (t->prev) + ZeroTrial = t->prev; + else + ZeroTrial = NULL; + } + if (t == Trials) { + Trials = t->next; + } + tfree(t); + } + else { - if (prev) - check_flat(prev, t); - else - NFlat = 1; + if (prev) + check_flat(prev, t); + else + NFlat = 1; - if (t->flags & ISAMINIMA) - t->flags &= ~ISAMINIMA; + if (t->flags & ISAMINIMA) + t->flags &= ~ISAMINIMA; - prev = t; + prev = t; #ifdef PZDEBUG - show_trial(t, '+'); + show_trial(t, '+'); #endif - } + } } @@ -860,21 +890,21 @@ CKTpzVerify(PZtrial **set, PZtrial *new_trial) * value to guess at if the search falls of the end of the list */ -static PZtrial * -pzseek(PZtrial *t, int dir) +static PZtrial* +pzseek(PZtrial* t, int dir) { Guess_Param = dir; if (t == NULL) - return NULL; + return NULL; if (dir == 0 && !(t->flags & ISAROOT) && !(t->flags & ISAMINIMA)) - return t; + return t; do { - if (dir >= 0) - t = t->next; - else - t = t->prev; + if (dir >= 0) + t = t->next; + else + t = t->prev; } while (t && ((t->flags & ISAROOT) || (t->flags & ISAMINIMA))); return t; @@ -883,139 +913,151 @@ pzseek(PZtrial *t, int dir) static void clear_trials(int mode) { - PZtrial *t, *next, *prev; + PZtrial* t, * next, * prev; prev = NULL; for (t = Trials; t; t = next) { - next = t->next; - if (mode || !(t->flags & ISAROOT)) { - tfree(t); - } else { - if (prev) - prev->next = t; - else - Trials = t; - t->prev = prev; - prev = t; - } + next = t->next; + if (mode || !(t->flags & ISAROOT)) { + tfree(t); + } + else { + if (prev) + prev->next = t; + else + Trials = t; + t->prev = prev; + prev = t; + } } if (prev) - prev->next = NULL; + prev->next = NULL; else - Trials = NULL; + Trials = NULL; } void -CKTpzUpdateSet(PZtrial **set, PZtrial *new) +CKTpzUpdateSet(PZtrial** set, PZtrial* new) { int this_move; this_move = 0; if (new->s.imag != 0.0) { - set[2] = set[1]; - set[1] = set[0]; - set[0] = new; - } else if (!set[1]) - set[1] = new; + set[2] = set[1]; + set[1] = set[0]; + set[0] = new; + } + else if (!set[1]) + set[1] = new; else if (!set[2] && new->s.real > set[1]->s.real) { - set[2] = new; - } else if (!set[0]) { - set[0] = new; - } else if (new->flags & ISAMINIMA) { - set[1] = new; - } else if (new->s.real < set[0]->s.real) { - set[2] = set[1]; - set[1] = set[0]; - set[0] = new; - this_move = FAR_LEFT; - } else if (new->s.real < set[1]->s.real) { - if (!CKTpzTrapped || new->mag_def < set[1]->mag_def - || (new->mag_def == set[1]->mag_def - && fabs(new->f_def.real) < fabs(set[1]->f_def.real))) { - /* Really should check signs, not just compare fabs( ) */ - set[2] = set[1]; /* XXX = set[2]->prev :: possible opt */ - set[1] = new; - this_move = MID_LEFT; - } else { - set[0] = new; - this_move = NEAR_LEFT; - } - } else if (new->s.real < set[2]->s.real) { - if (!CKTpzTrapped || new->mag_def < set[1]->mag_def - || (new->mag_def == set[1]->mag_def - && fabs(new->f_def.real) < fabs(set[1]->f_def.real))) { - /* Really should check signs, not just compare fabs( ) */ - set[0] = set[1]; - set[1] = new; - this_move = MID_RIGHT; - } else { - set[2] = new; - this_move = NEAR_RIGHT; - } - } else { - set[0] = set[1]; - set[1] = set[2]; - set[2] = new; - this_move = FAR_RIGHT; + set[2] = new; + } + else if (!set[0]) { + set[0] = new; + } + else if (new->flags & ISAMINIMA) { + set[1] = new; + } + else if (new->s.real < set[0]->s.real) { + set[2] = set[1]; + set[1] = set[0]; + set[0] = new; + this_move = FAR_LEFT; + } + else if (new->s.real < set[1]->s.real) { + if (!CKTpzTrapped || new->mag_def < set[1]->mag_def + || (new->mag_def == set[1]->mag_def + && fabs(new->f_def.real) < fabs(set[1]->f_def.real))) { + /* Really should check signs, not just compare fabs( ) */ + set[2] = set[1]; /* XXX = set[2]->prev :: possible opt */ + set[1] = new; + this_move = MID_LEFT; + } + else { + set[0] = new; + this_move = NEAR_LEFT; + } + } + else if (new->s.real < set[2]->s.real) { + if (!CKTpzTrapped || new->mag_def < set[1]->mag_def + || (new->mag_def == set[1]->mag_def + && fabs(new->f_def.real) < fabs(set[1]->f_def.real))) { + /* Really should check signs, not just compare fabs( ) */ + set[0] = set[1]; + set[1] = new; + this_move = MID_RIGHT; + } + else { + set[2] = new; + this_move = NEAR_RIGHT; + } + } + else { + set[0] = set[1]; + set[1] = set[2]; + set[2] = new; + this_move = FAR_RIGHT; } if (CKTpzTrapped && this_move == Last_Move) - Consec_Moves += 1; + Consec_Moves += 1; else - Consec_Moves = 0; + Consec_Moves = 0; Last_Move = this_move; } void -zaddeq(double *a, int *amag, double x, int xmag, double y, int ymag) +zaddeq(double* a, int* amag, double x, int xmag, double y, int ymag) { /* Balance magnitudes . . . */ if (xmag > ymag) { - *amag = xmag; - if (xmag > 50 + ymag) - y = 0.0; - else - for (xmag -= ymag; xmag > 0; xmag--) - y /= 2.0; - } else { - *amag = ymag; - if (ymag > 50 + xmag) - x = 0.0; - else - for (ymag -= xmag; ymag > 0; ymag--) - x /= 2.0; + *amag = xmag; + if (xmag > 50 + ymag) + y = 0.0; + else + for (xmag -= ymag; xmag > 0; xmag--) + y /= 2.0; + } + else { + *amag = ymag; + if (ymag > 50 + xmag) + x = 0.0; + else + for (ymag -= xmag; ymag > 0; ymag--) + x /= 2.0; } *a = x + y; if (*a == 0.0) - *amag = 0; + *amag = 0; else { - while (fabs(*a) > 1.0) { - *a /= 2.0; - *amag += 1; - } - while (fabs(*a) < 0.5) { - *a *= 2.0; - *amag -= 1; - } + while (fabs(*a) > 1.0) { + *a /= 2.0; + *amag += 1; + } + while (fabs(*a) < 0.5) { + *a *= 2.0; + *amag -= 1; + } } } #ifdef PZDEBUG static void -show_trial(PZtrial *new_trial, char x) +show_trial(PZtrial* new_trial, char x) { DEBUG(1) { - if(new_trial) { + if (new_trial) { fprintf(stderr, "%c (%3d/%3d) %.15g %.15g :: %.30g %.30g %d\n", x, - NIter, new_trial->seq_num, new_trial->s.real, new_trial->s.imag, - new_trial->f_def.real, new_trial->f_def.imag, new_trial->mag_def); + NIter, new_trial->seq_num, new_trial->s.real, new_trial->s.imag, + new_trial->f_def.real, new_trial->f_def.imag, new_trial->mag_def); if (new_trial->flags & ISANABERRATION) fprintf(stderr, "*** numerical aberration ***\n"); - } else { + } + else { fprintf(stderr, "%c (%3d/---) new_trial = nil\n", x, NIter); } } @@ -1023,7 +1065,7 @@ show_trial(PZtrial *new_trial, char x) #endif static void -check_flat(PZtrial *a, PZtrial *b) +check_flat(PZtrial* a, PZtrial* b) { int diff_mag; SPcomplex diff_frac; @@ -1031,162 +1073,169 @@ check_flat(PZtrial *a, PZtrial *b) diff_mag = a->mag_def - b->mag_def; if (abs(diff_mag) <= 1) { - if (diff_mag == 1) - mult = 2.0; - else if (diff_mag == -1) - mult = 0.5; - else - mult = 1.0; - C_SUBEQ(diff_frac, mult * a->f_def, b->f_def); - C_MAG2(diff_frac); - if (diff_frac.real < 1.0e-20) - NFlat += 1; + if (diff_mag == 1) + mult = 2.0; + else if (diff_mag == -1) + mult = 0.5; + else + mult = 1.0; + C_SUBEQ(diff_frac, mult * a->f_def, b->f_def); + C_MAG2(diff_frac); + if (diff_frac.real < 1.0e-20) + NFlat += 1; } - /* XXX else NFlat = ?????? */ + /* XXX else NFlat = ?????? */ } /* XXX ZZZ */ int -CKTpzStep(int strat, PZtrial **set) +CKTpzStep(int strat, PZtrial** set) { switch (strat) { case INIT: - if (!set[1]) { - set[1] = pzseek(ZeroTrial, 0); - } else if (!set[2]) - set[2] = pzseek(set[1], 1); - else if (!set[0]) - set[0] = pzseek(set[1], -1); - break; + if (!set[1]) { + set[1] = pzseek(ZeroTrial, 0); + } + else if (!set[2]) + set[2] = pzseek(set[1], 1); + else if (!set[0]) + set[0] = pzseek(set[1], -1); + break; case SKIP_LEFT: - set[0] = pzseek(set[0], -1); - break; + set[0] = pzseek(set[0], -1); + break; case SKIP_RIGHT: - set[2] = pzseek(set[2], 1); - break; + set[2] = pzseek(set[2], 1); + break; case SHIFT_LEFT: - set[2] = set[1]; - set[1] = set[0]; - set[0] = pzseek(set[0], -1); - break; + set[2] = set[1]; + set[1] = set[0]; + set[0] = pzseek(set[0], -1); + break; case SHIFT_RIGHT: - set[0] = set[1]; - set[1] = set[2]; - set[2] = pzseek(set[2], 1); - break; + set[0] = set[1]; + set[1] = set[2]; + set[2] = pzseek(set[2], 1); + break; } if (!set[0] || !set[1] || !set[2]) - return 0; + return 0; else - return 1; + return 1; } void -CKTpzReset(PZtrial **set) +CKTpzReset(PZtrial** set) { CKTpzTrapped = 0; Consec_Moves = 0; set[1] = pzseek(ZeroTrial, 0); if (set[1] != NULL) { - set[0] = pzseek(set[1], -1); - set[2] = pzseek(set[1], 1); - } else { - set[0] = NULL; - set[2] = NULL; + set[0] = pzseek(set[1], -1); + set[2] = pzseek(set[1], 1); + } + else { + set[0] = NULL; + set[2] = NULL; } } static int -alter(PZtrial *new, PZtrial *nearto, double abstol, double reltol) +alter(PZtrial* new, PZtrial* nearto, double abstol, double reltol) { double p1, p2; #ifdef PZDEBUG DEBUG(1) { - fprintf(stderr, "ALTER from: %.30g %.30g\n", - new->s.real, new->s.imag); - if (nearto->prev) - fprintf(stderr, "nt->prev %g\n", nearto->prev->s.real); - if (nearto->next) - fprintf(stderr, "nt->next %g\n", nearto->next->s.real); + fprintf(stderr, "ALTER from: %.30g %.30g\n", + new->s.real, new->s.imag); + if (nearto->prev) + fprintf(stderr, "nt->prev %g\n", nearto->prev->s.real); + if (nearto->next) + fprintf(stderr, "nt->next %g\n", nearto->next->s.real); } #endif if (CKTpzTrapped != 2) { #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, "not 2\n"); + DEBUG(1) fprintf(stderr, "not 2\n"); #endif - p1 = nearto->s.real; - if (nearto->flags & ISAROOT) - p1 -= 1e-6 * nearto->s.real + 1e-5; - if (nearto->prev) { - p1 += nearto->prev->s.real; + p1 = nearto->s.real; + if (nearto->flags & ISAROOT) + p1 -= 1e-6 * nearto->s.real + 1e-5; + if (nearto->prev) { + p1 += nearto->prev->s.real; #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, "p1 %g\n", p1); + DEBUG(1) fprintf(stderr, "p1 %g\n", p1); #endif - } else - p1 -= 10.0 * (fabs(p1) + 1.0); + } + else + p1 -= 10.0 * (fabs(p1) + 1.0); - p1 /= 2.0; - } else - p1 = nearto->s.real; + p1 /= 2.0; + } + else + p1 = nearto->s.real; if (CKTpzTrapped != 1) { #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, "not 1\n"); + DEBUG(1) fprintf(stderr, "not 1\n"); #endif - p2 = nearto->s.real; - if (nearto->flags & ISAROOT) - p2 += 1e-6 * nearto->s.real + 1e-5; /* XXX Would rather use pow(2)*/ - if (nearto->next) { - p2 += nearto->next->s.real; + p2 = nearto->s.real; + if (nearto->flags & ISAROOT) + p2 += 1e-6 * nearto->s.real + 1e-5; /* XXX Would rather use pow(2)*/ + if (nearto->next) { + p2 += nearto->next->s.real; #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, "p2 %g\n", p2); + DEBUG(1) fprintf(stderr, "p2 %g\n", p2); #endif - } else - p2 += 10.0 * (fabs(p2)+ 1.0); + } + else + p2 += 10.0 * (fabs(p2) + 1.0); - p2 /= 2.0; - } else - p2 = nearto->s.real; + p2 /= 2.0; + } + else + p2 = nearto->s.real; if ((nearto->prev && - fabs(p1 - nearto->prev->s.real) / - fabs(nearto->prev->s.real) + abstol/reltol < reltol) - || - (nearto->next && - fabs(p2 - nearto->next->s.real) / - fabs(nearto->next->s.real) + abstol/reltol < reltol)) { + fabs(p1 - nearto->prev->s.real) / + fabs(nearto->prev->s.real) + abstol / reltol < reltol) + || + (nearto->next && + fabs(p2 - nearto->next->s.real) / + fabs(nearto->next->s.real) + abstol / reltol < reltol)) { #ifdef PZDEBUG - DEBUG(1) - fprintf(stderr, "Bailed out\n"); + DEBUG(1) + fprintf(stderr, "Bailed out\n"); #endif - return 0; - } + return 0; + } if (CKTpzTrapped != 2 && nearto->s.real - p1 > p2 - nearto->s.real) { #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, "take p1\n"); + DEBUG(1) fprintf(stderr, "take p1\n"); #endif - new->s.real = p1; - } else { + new->s.real = p1; + } + else { #ifdef PZDEBUG - DEBUG(1) fprintf(stderr, "take p2\n"); + DEBUG(1) fprintf(stderr, "take p2\n"); #endif - new->s.real = p2; + new->s.real = p2; } #ifdef PZDEBUG DEBUG(1) fprintf(stderr, "ALTER to : %.30g %.30g\n", - new->s.real, new->s.imag); + new->s.real, new->s.imag); #endif return 1; From 4368790c5d8d324e8dedb3cc37e188ccd5fdd3df Mon Sep 17 00:00:00 2001 From: dwarning Date: Mon, 11 Sep 2023 14:50:27 +0200 Subject: [PATCH 12/13] remove compiler warning wrt. prototypes --- src/frontend/logicexp.c | 1 + src/spicelib/analysis/cktop.c | 2 -- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/frontend/logicexp.c b/src/frontend/logicexp.c index 32626be0e..13042bde9 100644 --- a/src/frontend/logicexp.c +++ b/src/frontend/logicexp.c @@ -614,6 +614,7 @@ static void ptable_print(PTABLE pt) /* Start of logicexp parser */ static char *get_inst_name(void); +static char *get_temp_name(void); static void aerror(char *s); static BOOL amatch(int t); static BOOL bexpr(void); diff --git a/src/spicelib/analysis/cktop.c b/src/spicelib/analysis/cktop.c index 6bebb681a..1b52c296b 100644 --- a/src/spicelib/analysis/cktop.c +++ b/src/spicelib/analysis/cktop.c @@ -16,8 +16,6 @@ Modified: 2005 Paolo Nenzi - Restructured #include "ngspice/enh.h" #endif -extern bool ft_ngdebug; - static int dynamic_gmin(CKTcircuit *, long int, long int, int); static int spice3_gmin(CKTcircuit *, long int, long int, int); static int new_gmin(CKTcircuit*, long int, long int, int); From acb7f2cd97aef6e8d9aa292293ffeb638d4d922a Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Mon, 11 Sep 2023 17:42:38 +0200 Subject: [PATCH 13/13] Make .ac error messages more verbose, prevent some crash, enable default values. --- src/spicelib/parser/inp2dot.c | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/src/spicelib/parser/inp2dot.c b/src/spicelib/parser/inp2dot.c index e1e545a2e..d6ffe39fe 100644 --- a/src/spicelib/parser/inp2dot.c +++ b/src/spicelib/parser/inp2dot.c @@ -183,6 +183,8 @@ dot_ac(char *line, CKTcircuit *ckt, INPtables *tab, struct card *current, int which; /* which analysis we are performing */ char *steptype; /* ac analysis, type of stepping function */ bool pdef = FALSE; /* issue a warning if default parameters are used */ + char* mline = line; /* for debug printout */ + double startval, stopval; NG_IGNORE(gnode); @@ -195,7 +197,7 @@ dot_ac(char *line, CKTcircuit *ckt, INPtables *tab, struct card *current, IFC(newAnalysis, (ckt, which, "AC Analysis", &foo, task)); INPgetTok(&line, &steptype, 1); /* get DEC, OCT, or LIN */ if (!*steptype || (!ciprefix("dec", steptype) && !ciprefix("oct", steptype) && !ciprefix("lin", steptype))) { - current->error = "Missing DEC, OCT, or LIN\n"; + LITERR("Missing DEC, OCT, or LIN.\n"); return (0); } ptemp.iValue = 1; @@ -203,23 +205,35 @@ dot_ac(char *line, CKTcircuit *ckt, INPtables *tab, struct card *current, tfree(steptype); parm = INPgetValue(ckt, &line, IF_INTEGER, tab); /* number of points */ - if (parm->iValue == 0) + if (parm->iValue < 1) { pdef = TRUE; + parm->iValue = 10; + } GCA(INPapName, (ckt, which, foo, "numsteps", parm)); - if(!isdigit(*line)) + if((*line != '.' && !isdigit(*line)) || (*line == '.' && !isdigit(*(line+1)))) pdef = TRUE; parm = INPgetValue(ckt, &line, IF_REAL, tab); /* fstart */ + startval = parm->rValue; + if (startval <= 0) { + pdef = TRUE; + startval = parm->rValue = 1.; + } GCA(INPapName, (ckt, which, foo, "start", parm)); - parm = INPgetValue(ckt, &line, IF_REAL, tab); /* fstop */ - if (parm->rValue == 0) + if((*line != '.' && !isdigit(*line)) || (*line == '.' && !isdigit(*(line+1)))) pdef = TRUE; + parm = INPgetValue(ckt, &line, IF_REAL, tab); /* fstop */ + stopval = parm->rValue; + if (stopval < startval) { + pdef = TRUE; + parm->rValue = 1000. * startval; + } GCA(INPapName, (ckt, which, foo, "stop", parm)); if (pdef) { fprintf(stderr, "Warning, ngspice assumes default parameter(s) for ac simulation\n"); - fprintf(stderr, " Check your ac or .ac line\n\n"); + fprintf(stderr, " Check your input line '.ac %s'\n\n", mline); } return (0); }