cktpzstr.c

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

		    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;
	    }
	}

    } while (shifted);

    if (pretest) {

#ifdef PZDEBUG
	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);
	}
#endif
	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 {
#endif
	    p->flags |= ISAMINIMA;
	    free(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);
#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;
	}

    }

    if (!repeat) {
	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);
		}
#endif
		error = SMPcLUfac(ckt->CKTmatrix, ckt->CKTpivotAbsTol);
		if (error == E_SINGULAR) {
#ifdef PZDEBUG
		    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. */,
		    &(((PZAN *) ckt->CKTcurJob)->PZnumswaps));
	    }

	    if (error != E_SINGULAR) {
		ckt->CKTniState &= ~NIPZSHOULDREORDER;
#ifdef PZDEBUG
		DEBUG(3) {
		    printf("Factored:\n");
		    SMPprint(ckt->CKTmatrix, stdout);
		}
#endif
		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 {

	    /* PZnumswaps is either 0 or 1 */
	    new_trial->f_raw.real *= ((PZAN *) ckt->CKTcurJob)->PZnumswaps;
	    new_trial->f_raw.imag *= ((PZAN *) ckt->CKTcurJob)->PZnumswaps;
	    /*
	    printf("SMP Det: (%g,%g)^%d\n", new_trial->f_raw.real,
		    new_trial->f_raw.imag, 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);
	}

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

	if (!(new_trial->flags & ISAROOT)) {
	    if (match)
		check_flat(match, new_trial);
	    else
		NFlat = 1;
	}
    }

#ifdef PZDEBUG
    show_trial(new_trial, '*');
#endif

    return OK;
}

/* Process a zero; inc. zero count, deflate other trials */

CKTpzVerify(set, new_trial)
    PZtrial	**set;
    PZtrial	*new_trial;
{
    PZtrial	*next;
    int		diff_mag;
    SPcomplex	diff_frac;
    double	tdiff;

    PZtrial *t, *prev;

    NZeros += 1;
    if (new_trial->s.imag != 0.0)
	NZeros += 1;
    NFlat = 0;

    if (new_trial->multiplicity == 0) {
	new_trial->flags |= ISAROOT;
	new_trial->multiplicity = 1;
    }

    prev = NULL;

    for (t = Trials; t; t = next) {

	next = t->next;

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

	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;
#ifdef PZDEBUG
	    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;
	    }
	    free(t);
	} else {

	    if (prev)
		check_flat(prev, t);
	    else
		NFlat = 1;

	    if (t->flags & ISAMINIMA)
		t->flags &= ~ISAMINIMA;

	    prev = t;
#ifdef PZDEBUG
	    show_trial(t, '+');
#endif
	}

    }

    return 1;	/* always ok */
}

/* pzseek: search the trial list (given a starting point) for the first
 *	non-zero entry; direction: -1 for prev, 1 for next, 0 for next
 *	-or- first.  Also, sets "Guess_Param" at the next reasonable
 *	value to guess at if the search falls of the end of the list
 */

static PZtrial *
pzseek(t, dir)
    PZtrial	*t;
    int		dir;
{
    Guess_Param = dir;
    if (t == NULL)
	return NULL;

    if (dir == 0 && !(t->flags & ISAROOT) && !(t->flags & ISAMINIMA))
	return t;

    do {
	if (dir >= 0)
	    t = t->next;
	else
	    t = t->prev;
    } while (t && ((t->flags & ISAROOT) || (t->flags & ISAMINIMA)));

    return t;
}

static void
clear_trials(mode)
    int mode;
{
    PZtrial *t, *next, *prev;

    prev = NULL;

    for (t = Trials; t; t = next) {
	next = t->next;
	if (mode || !(t->flags & ISAROOT)) {
	    free(t);
	} else {
	    if (prev)
		prev->next = t;
	    else
		Trials = t;
	    t->prev = prev;
	    prev = t;
	}
    }

    if (prev)
	prev->next = NULL;
    else
	Trials = NULL;
}

void
CKTpzUpdateSet(set, new)
    PZtrial	*set[3], *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;
    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;
    }

    if (CKTpzTrapped && this_move == Last_Move)
	Consec_Moves += 1;
    else
	Consec_Moves = 0;
    Last_Move = this_move;
}

void
zaddeq(a, amag, x, xmag, y, ymag)
    double	*a, x, y;
    int		*amag, xmag, 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;
    }

    *a = x + y;
    if (*a == 0.0)
	*amag = 0;
    else {
	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(new_trial, x)
    PZtrial *new_trial;
    char	x;
{
    DEBUG(1) 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);
    DEBUG(1)
	if (new_trial->flags & ISANABERRATION) {
	    fprintf(stderr, "*** numerical aberration ***\n");
	}
}
#endif

static void
check_flat(a, b)
    PZtrial	*a, *b;
{
    int		diff_mag;
    SPcomplex	diff_frac;
    double	mult;

    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;
    }
	/* XXX else NFlat = ?????? */
}

/* XXX ZZZ */
int
CKTpzStep(strat, set)
    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;

    case SKIP_LEFT:
	set[0] = pzseek(set[0], -1);
	break;

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

    case SHIFT_RIGHT:
	set[0] = set[1];
	set[1] = set[2];
	set[2] = pzseek(set[2], 1);
	break;

    }
    if (!set[0] || !set[1] || !set[2])
	return 0;
    else
	return 1;
}

void
CKTpzReset(set)
    PZtrial	*set[3];
{
    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;
    }
}

static int
alter(new, nearto, abstol, reltol)
    PZtrial	*new, *nearto;
    double	abstol, reltol;
{
    double	p1, p2;

#ifdef PZDEBUG
    DEBUG(1) fprintf(stderr, "ALTER from: %.30g %.30g\n",
	new->s.real, new->s.imag);
    DEBUG(1) fprintf(stderr, "nt->next %g\n", nearto->prev->s.real);
    DEBUG(1) fprintf(stderr, "nt->next %g\n", nearto->next->s.real);
#endif

    if (CKTpzTrapped != 2) {
#ifdef PZDEBUG
	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;
#ifdef PZDEBUG
	    DEBUG(1) fprintf(stderr, "p1 %g\n", p1);
#endif
	} else
	    p1 -= 10.0 * (FABS(p1) + 1.0);

	p1 /= 2.0;
    } else
	p1 = nearto->s.real;

    if (CKTpzTrapped != 1) {
#ifdef PZDEBUG
	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;
#ifdef PZDEBUG
	    DEBUG(1) fprintf(stderr, "p2 %g\n", p2);
#endif
	} else
	    p2 += 10.0 * (FABS(p2)+ 1.0);

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

#ifdef PZDEBUG
	DEBUG(1)
		fprintf(stderr, "Bailed out\n");
#endif

	return 0;
	}

    if (CKTpzTrapped != 2 && nearto->s.real - p1 > p2 - nearto->s.real) {
#ifdef PZDEBUG
	DEBUG(1) fprintf(stderr, "take p1\n");
#endif
	new->s.real = p1;
    } else {
#ifdef PZDEBUG
	DEBUG(1) fprintf(stderr, "take p2\n");
#endif
	new->s.real = p2;
    }

#ifdef PZDEBUG
    DEBUG(1) fprintf(stderr, "ALTER to  : %.30g %.30g\n",
	new->s.real, new->s.imag);
#endif
    return 1;

}