pelhomot.cc

Gambit 是一个游戏库理论软件

	dx[i + 2] = *da * dx[i + 2];
	dx[i + 3] = *da * dx[i + 3];
	dx[i + 4] = *da * dx[i + 4];
/* L50: */
    }
    return 0;
} /* dscal_ */

/* end dscal.c */

/**********************************************************************/
/******************* implementation code from f.c *********************/
/**********************************************************************/

/* Subroutine */ int f_(doublereal *x, 
			doublereal *v)
{

/* EVALUATE  F(X)  AND RETURN IN THE VECTOR  V . */

    /* Parameter adjustments */
    --v;
    --x;

    /* Function Body */
    return 0;
} /* f_ */

/* end f.c */

/**********************************************************************/
/***************** implementation code from fjac.c ********************/
/**********************************************************************/

/* Subroutine */ int fjac_(doublereal *x, 
			   doublereal *v, 
			   integer    *k)
{

/* RETURN IN  V  THE KTH COLUMN OF THE JACOBIAN MATRIX OF */
/* F(X) EVALUATED AT  X . */

    /* Parameter adjustments */
    --v;
    --x;

    /* Function Body */
    return 0;
} /* fjac_ */

/* end fjac.c */

/**********************************************************************/
/***************** implementation code from idamax.c ******************/
/**********************************************************************/

integer idamax_(integer    *n, 
		doublereal *dx, 
		integer    *incx)
{
    /* System generated locals */
    integer ret_val, i__1;
    doublereal d__1;

    /* Local variables */
    static doublereal dmax_;
    static integer i, ix;


/*     finds the index of element having max. absolute value. */
/*     jack dongarra, linpack, 3/11/78. */
/*     modified 3/93 to return if incx .le. 0. */


    /* Parameter adjustments */
    --dx;

    /* Function Body */
    ret_val = 0;
    if (*n < 1 || *incx <= 0) {
	return ret_val;
    }
    ret_val = 1;
    if (*n == 1) {
	return ret_val;
    }
    if (*incx == 1) {
	goto L20;
    }

/*        code for increment not equal to 1 */

    ix = 1;
    dmax_ = dblabs(dx[1]);
    ix += *incx;
    i__1 = *n;
    for (i = 2; i <= i__1; ++i) {
	if ((d__1 = dx[ix], dblabs(d__1)) <= dmax_) {
	    goto L5;
	}
	ret_val = i;
	dmax_ = (d__1 = dx[ix], dblabs(d__1));
L5:
	ix += *incx;
/* L10: */
    }
    return ret_val;

/*        code for increment equal to 1 */

L20:
    dmax_ = dblabs(dx[1]);
    i__1 = *n;
    for (i = 2; i <= i__1; ++i) {
	if ((d__1 = dx[i], dblabs(d__1)) <= dmax_) {
	    goto L30;
	}
	ret_val = i;
	dmax_ = (d__1 = dx[i], dblabs(d__1));
L30:
	;
    }
    return ret_val;
} /* idamax_ */

/* end idamax.c */

/**********************************************************************/
/****************** implementation code from root.c *******************/
/**********************************************************************/

double d_sign(double *arg1, double *arg2){
   if (*arg2>=0.0) return dblabs(*arg1);
   else return -1.0*dblabs(*arg1);
}

/* Table of constant values */

/* static integer c__4 = 4; NOW REDUNDANT */
static doublereal c_b17 = 1.;

/* Subroutine */ int root_(doublereal *t, 
			   doublereal *ft, 
			   doublereal *b, 
			   doublereal *c, 
			   doublereal *relerr, 
			   doublereal *abserr, 
			   integer    *iflag)
{
    /* System generated locals */
    doublereal d__1, d__2;

    /* Builtin functions */
    double d_sign(double *arg1, double *arg2);

    /* Local variables */
    static doublereal acmb, acbs, a, p, q, u;
    extern doublereal d1mach_(integer *i);
    static integer kount;
    static doublereal ae, fa, fb, fc;
    static integer ic;
    static doublereal re, fx, cmb, tol;


/*  ROOT COMPUTES A ROOT OF THE NONLINEAR EQUATION F(X)=0 */
/*  WHERE F(X) IS A CONTINOUS REAL FUNCTION OF A SINGLE REAL */
/*  VARIABLE X.  THE METHOD USED IS A COMBINATION OF BISECTION */
/*  AND THE SECANT RULE. */

/*  NORMAL INPUT CONSISTS OF A CONTINUOS FUNCTION F AND AN */
/*  INTERVAL (B,C) SUCH THAT F(B)*F(C).LE.0.0.  EACH ITERATION */
/*  FINDS NEW VALUES OF B AND C SUCH THAT THE INTERVAL(B,C) IS */
/*  SHRUNK AND F(B)*F(C).LE.0.0.  THE STOPPING CRITERION IS */

/*          DABS(B-C).LE.2.0*(RELERR*DABS(B)+ABSERR) */

/*  WHERE RELERR=RELATIVE ERROR AND ABSERR=ABSOLUTE ERROR ARE */
/*  INPUT QUANTITIES.  SET THE FLAG, IFLAG, POSITIVE TO INITIALIZE */
/*  THE COMPUTATION.  AS B,C AND IFLAG ARE USED FOR BOTH INPUT AND */
/*  OUTPUT, THEY MUST BE VARIABLES IN THE CALLING PROGRAM. */

/*  IF 0 IS A POSSIBLE ROOT, ONE SHOULD NOT CHOOSE ABSERR=0.0. */

/*  THE OUTPUT VALUE OF B IS THE BETTER APPROXIMATION TO A ROOT */
/*  AS B AND C ARE ALWAYS REDEFINED SO THAT DABS(F(B)).LE.DABS(F(C)). */

/*  TO SOLVE THE EQUATION, ROOT MUST EVALUATE F(X) REPEATEDLY. THIS */
/*  IS DONE IN THE CALLING PROGRAM.  WHEN AN EVALUATION OF F IS */
/*  NEEDED AT T, ROOT RETURNS WITH IFLAG NEGATIVE.  EVALUATE FT=F(T) */
/*  AND CALL ROOT AGAIN.  DO NOT ALTER IFLAG. */

/*  WHEN THE COMPUTATION IS COMPLETE, ROOT RETURNS TO THE CALLING */
/*  PROGRAM WITH IFLAG POSITIVE= */

/*     IFLAG=1  IF F(B)*F(C).LT.0 AND THE STOPPING CRITERION IS MET. */

/*          =2  IF A VALUE B IS FOUND SUCH THAT THE COMPUTED VALUE */
/*              F(B) IS EXACTLY ZERO.  THE INTERVAL (B,C) MAY NOT */
/*              SATISFY THE STOPPING CRITERION. */

/*          =3  IF DABS(F(B)) EXCEEDS THE INPUT VALUES DABS(F(B)), */
/*              DABS(F(C)).  IN THIS CASE IT IS LIKELY THAT B IS CLOSE */
/*              TO A POLE OF F. */

/*          =4  IF NO ODD ORDER ROOT WAS FOUND IN THE INTERVAL.  A */
/*              LOCAL MINIMUM MAY HAVE BEEN OBTAINED. */

/*          =5  IF TOO MANY FUNCTION EVALUATIONS WERE MADE. */
/*              (AS PROGRAMMED, 500 ARE ALLOWED.) */

/*  THIS CODE IS A MODIFICATION OF THE CODE ZEROIN WHICH IS COMPLETELY */
/*  EXPLAINED AND DOCUMENTED IN THE TEXT  NUMERICAL COMPUTING:  AN */
/*  INTRODUCTION,  BY L. F. SHAMPINE AND R. C. ALLEN. */

/*  CALLS  D1MACH . */


    if (*iflag >= 0) {
	goto L100;
    }
    *iflag = dblabs(*iflag);
    switch ((int)*iflag) {
	case 1:  goto L200;
	case 2:  goto L300;
	case 3:  goto L400;
    }
L100:
    u = d1mach_(&c__4);
    re = max(*relerr,u);
    ae = max(*abserr,0.);
    ic = 0;
    acbs = (d__1 = *b - *c, dblabs(d__1));
    a = *c;
    *t = a;
    *iflag = -1;
    return 0;
L200:
    fa = *ft;
    *t = *b;
    *iflag = -2;
    return 0;
L300:
    fb = *ft;
    fc = fa;
    kount = 2;
/* Computing MAX */
    d__1 = dblabs(fb), d__2 = dblabs(fc);
    fx = max(d__1,d__2);
L1:
    if (dblabs(fc) >= dblabs(fb)) {
	goto L2;
    }

/*  INTERCHANGE B AND C SO THAT ABS(F(B)).LE.ABS(F(C)). */

    a = *b;
    fa = fb;
    *b = *c;
    fb = fc;
    *c = a;
    fc = fa;
L2:
    cmb = (*c - *b) * (float).5;
    acmb = dblabs(cmb);
    tol = re * dblabs(*b) + ae;

/*  TEST STOPPING CRITERION AND FUNCTION COUNT. */

    if (acmb <= tol) {
	goto L8;
    }
    if (kount >= 500) {
	goto L12;
    }

/*  CALCULATE NEW ITERATE EXPLICITLY AS B+P/Q */
/*  WHERE WE ARRANGE P.GE.0.  THE IMPLICIT */
/*  FORM IS USED TO PREVENT OVERFLOW. */

    p = (*b - a) * fb;
    q = fa - fb;
    if (p >= (float)0.) {
	goto L3;
    }
    p = -p;
    q = -q;

/*  UPDATE A, CHECK IF REDUCTION IN THE SIZE OF BRACKETING */
/*  INTERVAL IS SATISFACTORY.  IF NOT BISECT UNTIL IT IS. */

L3:
    a = *b;
    fa = fb;
    ++ic;
    if (ic < 4) {
	goto L4;
    }
    if (acmb * (float)8. >= acbs) {
	goto L6;
    }
    ic = 0;
    acbs = acmb;

/*  TEST FOR TOO SMALL A CHANGE. */

L4:
    if (p > dblabs(q) * tol) {
	goto L5;
    }

/*  INCREMENT BY TOLERANCE */

/* sign gives sign(arg2)*|arg1| */
    *b += d_sign(&tol, &cmb);
    goto L7;

/*  ROOT OUGHT TO BE BETWEEN B AND (C+B)/2 */

L5:
    if (p >= cmb * q) {
	goto L6;
    }

/*  USE SECANT RULE. */

    *b += p / q;
    goto L7;

/*  USE BISECTION. */

L6:
    *b = (*c + *b) * (float).5;

/*  HAVE COMPLETED COMPUTATION FOR NEW ITERATE B. */

L7:
    *t = *b;
    *iflag = -3;
    return 0;
L400:
    fb = *ft;
    if (fb == (float)0.) {
	goto L9;
    }
    ++kount;
    if (d_sign(&c_b17, &fb) != d_sign(&c_b17, &fc)) {
	goto L1;
    }
    *c = a;
    fc = fa;
    goto L1;

/* FINISHED.  SET IFLAG. */

L8:
    if (d_sign(&c_b17, &fb) == d_sign(&c_b17, &fc)) {
	goto L11;
    }
    if (dblabs(fb) > fx) {
	goto L10;
    }
    *iflag = 1;
    return 0;
L9:
    *iflag = 2;
    return 0;
L10:
    *iflag = 3;
    return 0;
L11:
    *iflag = 4;
    return 0;
L12:
    *iflag = 5;
    return 0;
} /* root_ */

/* end root.c */

/**********************************************************************/
/***************** implementation code from rootnf.c ******************/
/**********************************************************************/

/* Table of constant values */

/*
static int c__4 = 4;
static int c__1 = 1;
*/

/* Subroutine */ int rootnf_(int    *n, 
			     int    *nfe, 
			     int    *iflag, 
			     double *relerr,
			     double *abserr,
			     double *y,
			     double *yp, 
			     double *yold, 
			     double *ypold, 
			     double *a,
			     double *qr, 
			     double *alpha, 
			     double *tz, 
			     int    *pivot, 
			     double *w, 
			     double *wp, 
			     double *par, 
			     int    *ipar)
{
    /* System generated locals */
    int qr_dim1, qr_offset, i__1;
    double d__1;

    /* Local variables */
    static double dels, aerr, rerr;
    static int judy;
    extern /* Subroutine */ int root_(doublereal *t, 
				      doublereal *ft, 
				      doublereal *b, 
				      doublereal *c, 
				      doublereal *relerr, 
				      doublereal *abserr, 
				      integer    *iflag);
    static double sout;
    extern double dnrm2_(integer    *n, 
			 doublereal *dx, 
			 integer    *incx);
    static double u;
    static int lcode;
    extern double d1mach_(integer *);
    static double qsout, sa, sb;
    static int jw;
    /* extern */ /* Subroutine */ /* int tangnf_(); IN pelutils.h */
    static int np1;


/* ROOTNF  FINDS THE POINT  YBAR = (1, XBAR)  ON THE ZERO CURVE OF THE */
/* HOMOTOPY MAP.  IT STARTS WITH TWO POINTS YOLD=(LAMBDAOLD,XOLD) AND */
/* Y=(LAMBDA,X) SUCH THAT  LAMBDAOLD < 1 <= LAMBDA , AND ALTERNATES */
/* BETWEEN HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION UNTIL */
/* CONVERGENCE. */

/* ON INPUT: */

/* N = DIMENSION OF X AND THE HOMOTOPY MAP. */

/* NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS. */

/* IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE. */

/* RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES.  THE ITERATION IS 
*/
/*    CONSIDERED TO HAVE CONVERGED WHEN A POINT Y=(LAMBDA,X) IS FOUND */
/*    SUCH THAT */

/*    |Y(1) - 1| <= RELERR + ABSERR              AND */

/*    ||Z|| <= RELERR*||X|| + ABSERR  ,          WHERE */

/*    (?,Z) IS THE NEWTON STEP TO Y=(LAMBDA,X). */

/* Y(1:N+1) = POINT (LAMBDA(S), X(S)) ON ZERO CURVE OF HOMOTOPY MAP. */

/* YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY MAP 
*/
/*    AT  Y . */

/* YOLD(1:N+1) = A POINT DIFFERENT FROM  Y  ON THE ZERO CURVE. */

/* YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY */
/*    MAP AT  YOLD . */

/* A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP. */

/* QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1), W(1:N+1), */
/*    WP(1:N+1)  ARE WORK ARRAYS USED FOR THE QR FACTORIZATION (IN THE */
/*    NEWTON STEP CALCULATION) AND THE INTERPOLATION. */

/* PAR(1:*)