luke
/
ngspice
mirror of https://git.code.sf.net/p/ngspice/ngspice
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Formatting cktpzstr.c

This commit is contained in:
Holger Vogt 2023-09-11 14:28:49 +02:00
parent dfeb0bdb4e
commit b7993bb530
1 changed files with 762 additions and 713 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

283
src/spicelib/analysis/cktpzstr.c
View File

@ -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,7 +67,7 @@ static void show_trial(PZtrial *new_trial, char x);
#define MID_RIGHT 9
#ifdef PZDEBUG
static char *snames[ ] = {
static char* snames[] = {
"none",
"none",
"shift left",
@ -85,7 +85,7 @@ static char *snames[ ] = {
"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;
@ -167,23 +167,27 @@ CKTpzFindZeros(CKTcircuit *ckt, PZtrial **rootinfo, int *rootcount)
if (CKTpzVerify(neighborhood, new_trial)) {
NIter = 0;
CKTpzReset(neighborhood);
} else
}
else
/* XXX Verify fails ?!? */
CKTpzUpdateSet(neighborhood, new_trial);
} else if (new_trial->flags & ISANABERRATION) {
}
else if (new_trial->flags & ISANABERRATION) {
CKTpzReset(neighborhood);
Aberr_Num += 1;
tfree(new_trial);
} else if (new_trial->flags & ISAMINIMA) {
}
else if (new_trial->flags & ISAMINIMA) {
neighborhood[0] = NULL;
neighborhood[1] = new_trial;
neighborhood[2] = NULL;
} else {
}
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);
SPfrontEnd->IFerrorf(ERR_WARNING, "Pole-Zero analysis interrupted; %d trials, %d roots\n", Seq_Num, NZeros);
error = E_PAUSE;
break;
}
@ -207,18 +211,19 @@ CKTpzFindZeros(CKTcircuit *ckt, PZtrial **rootinfo, int *rootcount)
*rootinfo = NULL;
*rootcount = 0;
MERROR(E_SHORT, "The input signal is shorted on the way to the output");
} else
}
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;
@ -248,7 +253,8 @@ PZeval(int strat, PZtrial **set, PZtrial **new_trial_p)
Guess_Param = High_Guess * 10.0;
else
Guess_Param = 1.0;
} else {
}
else {
if (Low_Guess < 0.0)
Guess_Param = Low_Guess * 10.0;
else
@ -281,7 +287,8 @@ PZeval(int strat, PZtrial **set, PZtrial **new_trial_p)
#endif
new_trial->s.real = (set[0]->s.real + set[1]->s.real) / 2.0;
}
} else if (CKTpzTrapped == 2) {
}
else if (CKTpzTrapped == 2) {
if (new_trial->s.real < set[1]->s.real
|| new_trial->s.real > set[2]->s.real) {
#ifdef PZDEBUG
@ -292,7 +299,8 @@ PZeval(int strat, PZtrial **set, PZtrial **new_trial_p)
#endif
new_trial->s.real = (set[1]->s.real + set[2]->s.real) / 2.0;
}
} else if (CKTpzTrapped == 3) {
}
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)
@ -348,7 +356,8 @@ PZeval(int strat, PZtrial **set, PZtrial **new_trial_p)
NIpzK_mag += 1;
}
new_trial->s.imag = NIpzK;
} else
}
else
new_trial->s.imag = 10000.0;
/*
@ -364,7 +373,8 @@ PZeval(int strat, PZtrial **set, PZtrial **new_trial_p)
if (!set[2]) {
new_trial->s.real = set[0]->s.real;
new_trial->s.imag = 1.0e8;
} else {
}
else {
new_trial->s.real = set[0]->s.real;
new_trial->s.imag = 1.0e12;
}
@ -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;
@ -410,23 +420,28 @@ int CKTpzStrat(PZtrial **set)
if (set[1] && (set[1]->flags & ISAMINIMA)) {
suggestion = COMPLEX_INIT;
} else if (set[0] && set[0]->s.imag != 0.0) {
}
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]) {
}
else if (!set[0] || !set[1] || !set[2]) {
suggestion = INIT;
} else {
}
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)) {
}
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 {
}
else {
zaddeq(&a, &a_mag, set[1]->f_def.real, set[1]->mag_def,
-set[0]->f_def.real, set[0]->mag_def);
@ -445,7 +460,8 @@ int CKTpzStrat(PZtrial **set)
suggestion = SKIP_RIGHT;
else
suggestion = SKIP_LEFT;
} else if (sgn(a) != -sgn(b)) {
}
else if (sgn(a) != -sgn(b)) {
if (a == 0.0)
suggestion = SKIP_LEFT;
else if (b == 0.0)
@ -454,23 +470,27 @@ int CKTpzStrat(PZtrial **set)
suggestion = SHIFT_LEFT;
else
suggestion = SHIFT_RIGHT;
} else if (sgn(a) == -sgn(set[1]->f_def.real)) {
}
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)
}
else if (k1 > k2)
suggestion = SKIP_RIGHT;
else
suggestion = SKIP_LEFT;
} else {
}
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
}
else if (a_mag > b_mag || (a_mag == b_mag
&& fabs(a) > fabs(b)))
suggestion = SPLIT_LEFT;
else
@ -484,7 +504,7 @@ int CKTpzStrat(PZtrial **set)
else if (Last_Move == MID_RIGHT || Last_Move == NEAR_LEFT)
suggestion = SPLIT_RIGHT;
else
abort( ); /* XXX */
abort(); /* XXX */
Consec_Moves = 0;
}
}
@ -506,12 +526,12 @@ int CKTpzStrat(PZtrial **set)
/* CKTpzRunTrial: eval the function at a given 's', fold in deflation */
int
CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set)
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;
@ -544,7 +564,7 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set)
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)) {
@ -554,13 +574,14 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set)
if (p->flags & ISAROOT) {
abstol = 1e-5;
reltol = 1e-6;
} else {
}
else {
abstol = 1e-20;
reltol = 1e-12;
}
if (diff_frac.imag == 0.0 &&
fabs(diff_frac.real) / (fabs(p->s.real) + abstol/reltol)
fabs(diff_frac.real) / (fabs(p->s.real) + abstol / reltol)
< reltol) {
#ifdef PZDEBUG
@ -577,30 +598,33 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set)
p->count = 0;
pretest = 1;
break;
} else
}
else
p->count += 1; /* try to shift */
shifted = 1; /* Re-calculate deflation */
break;
} else {
}
else {
if (!CKTpzTrapped)
p->count = 0;
if (p->flags & ISAROOT) {
diff_mag = 0;
C_NORM(diff_frac,diff_mag);
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);
C_NORM(diff_frac, diff_mag);
for (i = p->multiplicity; i > 0; i--) {
C_MUL(def_frac,diff_frac);
C_MUL(def_frac, diff_frac);
def_mag += diff_mag;
C_NORM(def_frac,def_mag);
C_NORM(def_frac, def_mag);
}
} else if (!match)
}
else if (!match)
match = p;
}
}
@ -625,27 +649,29 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set)
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) {
}
else if (p->flags & ISAROOT) {
#ifdef PZDEBUG
DEBUG(1) fprintf(stderr, "Repeat at %.30g %.30g\n",
p->s.real, p->s.imag);
#endif
*new_trialp = p;
* new_trialp = p;
p->flags |= ISAREPEAT;
p->multiplicity += 1;
repeat = 1;
} else {
}
else {
/* Regular zero, as precise as we can get it */
error = E_SINGULAR;
}
@ -670,7 +696,8 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set)
DEBUG(1) printf("Needs reordering\n");
#endif
ckt->CKTniState |= NIPZSHOULDREORDER;
} else if (error != OK)
}
else if (error != OK)
return error;
}
if (ckt->CKTniState & NIPZSHOULDREORDER) {
@ -703,7 +730,8 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set)
new_trial->mag_def = 0;
new_trial->flags = ISAROOT;
/*printf("SMP Det: Singular\n");*/
} else if (error != OK)
}
else if (error != OK)
return error;
else {
@ -720,9 +748,9 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set)
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);
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_NORM(new_trial->f_def, new_trial->mag_def);
}
/* Link into the rest of the list */
@ -731,7 +759,8 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set)
if (prev->next)
prev->next->prev = new_trial;
prev->next = new_trial;
} else {
}
else {
if (Trials)
Trials->prev = new_trial;
else
@ -761,14 +790,14 @@ CKTpzRunTrial(CKTcircuit *ckt, PZtrial **new_trialp, PZtrial **set)
/* Process a zero; inc. zero count, deflate other trials */
int
CKTpzVerify(PZtrial **set, PZtrial *new_trial)
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);
@ -794,7 +823,7 @@ CKTpzVerify(PZtrial **set, PZtrial *new_trial)
continue;
}
C_SUBEQ(diff_frac,new_trial->s,t->s);
C_SUBEQ(diff_frac, new_trial->s, t->s);
if (new_trial->s.imag != 0.0)
C_MAG2(diff_frac);
@ -804,10 +833,10 @@ CKTpzVerify(PZtrial **set, PZtrial *new_trial)
*/
if (diff_frac.real != 0.0) {
diff_mag = 0;
C_NORM(diff_frac,diff_mag);
C_NORM(diff_frac, diff_mag);
diff_mag *= -1;
C_DIV(t->f_def,diff_frac);
C_NORM(t->f_def,diff_mag);
C_DIV(t->f_def, diff_frac);
C_NORM(t->f_def, diff_mag);
t->mag_def += diff_mag;
}
@ -833,7 +862,8 @@ CKTpzVerify(PZtrial **set, PZtrial *new_trial)
Trials = t->next;
}
tfree(t);
} else {
}
else {
if (prev)
check_flat(prev, t);
@ -860,8 +890,8 @@ 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)
@ -883,7 +913,7 @@ pzseek(PZtrial *t, int dir)
static void
clear_trials(int mode)
{
PZtrial *t, *next, *prev;
PZtrial* t, * next, * prev;
prev = NULL;
@ -891,7 +921,8 @@ clear_trials(int mode)
next = t->next;
if (mode || !(t->flags & ISAROOT)) {
tfree(t);
} else {
}
else {
if (prev)
prev->next = t;
else
@ -908,7 +939,7 @@ clear_trials(int mode)
}
void
CKTpzUpdateSet(PZtrial **set, PZtrial *new)
CKTpzUpdateSet(PZtrial** set, PZtrial* new)
{
int this_move;
@ -918,20 +949,25 @@ CKTpzUpdateSet(PZtrial **set, PZtrial *new)
set[2] = set[1];
set[1] = set[0];
set[0] = new;
} else if (!set[1])
}
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]) {
}
else if (!set[0]) {
set[0] = new;
} else if (new->flags & ISAMINIMA) {
}
else if (new->flags & ISAMINIMA) {
set[1] = new;
} else if (new->s.real < set[0]->s.real) {
}
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) {
}
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))) {
@ -939,11 +975,13 @@ CKTpzUpdateSet(PZtrial **set, PZtrial *new)
set[2] = set[1]; /* XXX = set[2]->prev :: possible opt */
set[1] = new;
this_move = MID_LEFT;
} else {
}
else {
set[0] = new;
this_move = NEAR_LEFT;
}
} else if (new->s.real < set[2]->s.real) {
}
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))) {
@ -951,11 +989,13 @@ CKTpzUpdateSet(PZtrial **set, PZtrial *new)
set[0] = set[1];
set[1] = new;
this_move = MID_RIGHT;
} else {
}
else {
set[2] = new;
this_move = NEAR_RIGHT;
}
} else {
}
else {
set[0] = set[1];
set[1] = set[2];
set[2] = new;
@ -970,7 +1010,7 @@ CKTpzUpdateSet(PZtrial **set, PZtrial *new)
}
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) {
@ -980,7 +1020,8 @@ zaddeq(double *a, int *amag, double x, int xmag, double y, int ymag)
else
for (xmag -= ymag; xmag > 0; xmag--)
y /= 2.0;
} else {
}
else {
*amag = ymag;
if (ymag > 50 + xmag)
x = 0.0;
@ -1006,16 +1047,17 @@ zaddeq(double *a, int *amag, double x, int xmag, double y, int ymag)
#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);
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;
@ -1047,13 +1089,14 @@ check_flat(PZtrial *a, PZtrial *b)
/* 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])
}
else if (!set[2])
set[2] = pzseek(set[1], 1);
else if (!set[0])
set[0] = pzseek(set[1], -1);
@ -1087,7 +1130,7 @@ CKTpzStep(int strat, PZtrial **set)
}
void
CKTpzReset(PZtrial **set)
CKTpzReset(PZtrial** set)
{
CKTpzTrapped = 0;
Consec_Moves = 0;
@ -1096,14 +1139,15 @@ CKTpzReset(PZtrial **set)
if (set[1] != NULL) {
set[0] = pzseek(set[1], -1);
set[2] = pzseek(set[1], 1);
} else {
}
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;
@ -1130,11 +1174,13 @@ alter(PZtrial *new, PZtrial *nearto, double abstol, double reltol)
#ifdef PZDEBUG
DEBUG(1) fprintf(stderr, "p1 %g\n", p1);
#endif
} else
}
else
p1 -= 10.0 * (fabs(p1) + 1.0);
p1 /= 2.0;
} else
}
else
p1 = nearto->s.real;
if (CKTpzTrapped != 1) {
@ -1149,20 +1195,22 @@ alter(PZtrial *new, PZtrial *nearto, double abstol, double reltol)
#ifdef PZDEBUG
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
}
else
p2 = nearto->s.real;
if ((nearto->prev &&
fabs(p1 - nearto->prev->s.real) /
fabs(nearto->prev->s.real) + abstol/reltol < reltol)
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(nearto->next->s.real) + abstol / reltol < reltol)) {
#ifdef PZDEBUG
DEBUG(1)
@ -1177,7 +1225,8 @@ alter(PZtrial *new, PZtrial *nearto, double abstol, double reltol)
DEBUG(1) fprintf(stderr, "take p1\n");
#endif
new->s.real = p1;
} else {
}
else {
#ifdef PZDEBUG
DEBUG(1) fprintf(stderr, "take p2\n");
#endif