cktpzstr.c

spice中支持多层次元件模型仿真的可单独运行的插件源码

/**********
Copyright 1990 Regents of the University of California.  All rights reserved.
**********/

/*
 *	A variant on the "zeroin" method.  This is a bit convoluted.
 */

#include "spice.h"
#include <stdio.h>
#include "pzdefs.h"
#include "util.h"
#include "complex.h"
#include "cktdefs.h"
#include "devdefs.h"
#include "sperror.h"
#include "suffix.h"

#ifdef PZDEBUG
#  ifndef notdef
#    define DEBUG(N)	if (Debug >= (unsigned) (N))
static unsigned int	Debug = 1;
#  else
#    define DEBUG(N)	if (0)
#  endif
#endif

static void	clear_trials( ), check_flat( );
void		CKTpzUpdateSet( ), zaddeq( ), CKTpzReset( );
#ifdef PZDEBUG
static void show_trial( );
#endif

#define NITER_LIM	200

#define	SHIFT_LEFT	2
#define	SHIFT_RIGHT	3
#define	SKIP_LEFT	4
#define	SKIP_RIGHT	5
#define	INIT		6

#define	GUESS		7
#define	SPLIT_LEFT	8
#define	SPLIT_RIGHT	9

#define	MULLER		10
#define	SYM		11
#define	SYM2		12
#define	COMPLEX_INIT	13
#define	COMPLEX_GUESS	14
#define	QUIT		15

#define	NEAR_LEFT	4
#define	MID_LEFT	5
#define	FAR_LEFT	6
#define	NEAR_RIGHT	7
#define	FAR_RIGHT	8
#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"
	};
#endif

#define sgn(X)	((X) < 0 ? -1 : (X) == 0 ? 0 : 1)

#define	ISAROOT		2
#define	ISAREPEAT	4
#define	ISANABERRATION	8
#define	ISAMINIMA	16

extern double	NIpzK;
extern int	NIpzK_mag;

int	CKTpzTrapped;

static int	NZeros, NFlat, Max_Zeros;
static PZtrial	*ZeroTrial, *Trials;
static int	Seq_Num;
static double	Guess_Param;
static double	High_Guess, Low_Guess;
static int	Last_Move, Consec_Moves;
static int	NIter, NTrials;
static int	Aberr_Num;

int PZeval( );
static PZtrial *pzseek( );
static int alter( );

CKTpzFindZeros(ckt, rootinfo, rootcount)
    CKTcircuit  *ckt;
    PZtrial	**rootinfo;
    int		*rootcount;
{
    PZtrial *new_trial;
    PZtrial *neighborhood[3];
    int    done = 0;
    int    strat;
    int    error;
    char   ebuf[513];

    NIpzK = 0.0;
    NIpzK_mag = 0;
    High_Guess = -1.0;
    Low_Guess = 1.0;
    ZeroTrial = 0;
    Trials = 0;
    NZeros = 0;
    NFlat = 0;
    Max_Zeros = SMPmatSize(ckt->CKTmatrix);
    NIter = 0;
    error = OK;
    CKTpzTrapped = 0;
    Aberr_Num = 0;
    NTrials = 0;
    ckt->CKTniState |= NIPZSHOULDREORDER; /* Initial for LU fill-ins */

    Seq_Num = 1;

    CKTpzReset(neighborhood);

    do {

	while ((strat = CKTpzStrat(neighborhood)) < GUESS && !CKTpzTrapped)
	    if (!CKTpzStep(strat, neighborhood)) {
		strat = GUESS;
#ifdef PZDEBUG
		DEBUG(1) fprintf(stderr, "\t\tGuess\n");
#endif
		break;
	    }

	NIter += 1;
	
	/* Evaluate current strategy */
	error = PZeval(strat, neighborhood, &new_trial);
	if (error != OK)
	    return error;

	error = CKTpzRunTrial(ckt, &new_trial, neighborhood);
	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;
	    free(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))( )) {
	    sprintf(ebuf,
		"Pole-Zero analysis interrupted; %d trials, %d roots\n",
		Seq_Num, NZeros); 
	    (*(SPfrontEnd->IFerror))(ERR_WARNING, ebuf, 0);
	    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 */

#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);
#endif

    if (NZeros >= Seq_Num - 1) {
	/* Short */
	clear_trials(ISAROOT);
	*rootinfo = NULL;
	*rootcount = 0;
	ERROR(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) {
	sprintf(ebuf,
    "Pole-zero converging to numerical aberrations; giving up after %d trials",
	    Seq_Num);
	(*(SPfrontEnd->IFerror))(ERR_WARNING, ebuf, 0);
    }

    if (NIter >= NITER_LIM) {
	sprintf(ebuf,
	    "Pole-zero iteration limit reached; giving up after %d trials",
	    Seq_Num);
	(*(SPfrontEnd->IFerror))(ERR_WARNING, ebuf, 0);
    }

    return error;
}

/* PZeval: evaluate an estimation function (given by 'strat') for the next
    guess (returned in a PZtrial) */

/* XXX ZZZ */
int
PZeval(strat, set, new_trial_p)
    int	strat;
    PZtrial	*set[ ];
    PZtrial	**new_trial_p;
{
    int		error;
    double	a, b, est, k, n;
    PZtrial	*new_trial;

    new_trial = NEW(PZtrial);
    new_trial->multiplicity = 0;
    new_trial->count = 0;
    new_trial->seq_num = Seq_Num++;

    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 (Guess_Param > High_Guess)
	    High_Guess = Guess_Param;
	if (Guess_Param < Low_Guess)
	    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);

	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);
#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) {
#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);
#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) {
#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);
#endif
		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!");
#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;
		}
	    }
	}

	break;

    case COMPLEX_INIT:
	/* Not automatic */
#ifdef PZDEBUG
	DEBUG(1) fprintf(stderr, "\tPZ minima at: %-30g %d\n",
	    NIpzK, NIpzK_mag);
#endif

	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".
	 */

	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.
	 */

	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;

    case MULLER:
	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;

    case SPLIT_RIGHT:
	new_trial->s.real = (set[2]->s.real + 2 * set[1]->s.real) / 3.0;
	error = OK;
	break;

    default:
	ERROR(E_PANIC, "Step type unkown");
	break;
    }

    *new_trial_p = new_trial;
    return error;
}

/* CKTpzStrat: given three points, determine a good direction or method for
    guessing the next zero */

/* XXX ZZZ what is a strategy for complex hunting? */
CKTpzStrat(set)
    PZtrial	**set;
{
    int	suggestion;
    double a, b, q;
    double d1, d2;
    double diff;
    int diff_mag;
    int	a_mag, b_mag, q_mag;
    double x, y;
    int	xmag;
    double k1, k2;
    int	new_trap;

    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 {

	    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) {

		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]);
    }
#endif
    return suggestion;
}

/* CKTpzRunTrial: eval the function at a given 's', fold in deflation */

CKTpzRunTrial(ckt, new_trialp, set)
    CKTcircuit	*ckt;
    PZtrial	**new_trialp;
    PZtrial	*set[ ];
{
    PZtrial	*t, *match, *base, *new_trial;
    PZtrial	*p, *prev;
    SPcomplex	def_frac, diff_frac;
    double	num, den;
    double	reltol, abstol;
    int		e, def_mag, diff_mag, error;
    int		inserted, i;
    int		pretest, shifted, was_shifted;
    int		repeat;

    new_trial = *new_trialp;

    if (new_trial->s.imag < 0.0)
	new_trial->s.imag *= -1.0;

    /* Insert the trial into the list of Trials, while calculating
	the deflation factor from previous zeros */

    pretest = 0;
    shifted = 0;
    repeat = 0;

    do {

	def_mag = 0;
	def_frac.real = 1.0;
	def_frac.imag = 0.0;
	was_shifted = shifted;
	shifted = 0;

	prev = NULL;
	base = NULL;
	match = NULL;

	for (p = Trials; p != NULL; p = p->next) {

	    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 (p->flags & ISAMINIMA)
		    base = p;
	    }

	    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) {

#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);
		    }
#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 */

		shifted = 1;	/* Re-calculate deflation */
		break;

	    } else {
		if (!CKTpzTrapped)
		    p->count = 0;
		if (p->flags & ISAROOT) {