pelhomot.cc

Gambit 是一个游戏库理论软件

/* Homotopies.c */

/*
  This is the implementation code from the various files that used 
to reside in the Cont subdirectory of Pelican0.80/source.
*/

#include "pelhomot.h"

/**********************************************************************/
/*************** implementation code from Hom_params.c ****************/
/**********************************************************************/

double Hom_tol=.000001;  /* tolerence for path tracking */

void print_Hom_params(FILE *outfile){
#ifdef HOM_PRINT
 fprintf(outfile," Hominuation Parameters:\n");
 fprintf(outfile,"      Hom_tol=%g\n",Hom_tol)
#endif
;
}

/* end Hom_params.c */

/************************************************************************/
/***************************** code from fixpnf.c ***********************/
/************************************************************************/

#define TRUE_ (1)
#define FALSE_ (0)
double get_abs_homog();
double coord_r, coord_i;

/* Table of constant values */

static integer c__4 = 4;
static integer c__1 = 1;

/* Subroutine */ int fixpnf_0_(int     n__, 
			       int    *n, 
			       double *y, 
			       int    *iflag, 
			       double *arcre, 
			       double *arcae, 
			       double *ansre, 
			       double *ansae, 
			       int    *trace, 
			       double *a, 
			       int    *nfe, 
			       double *arclen, 
			       double *yp, 
			       double *yold, 
			       double *ypold, 
			       double *qr, 
			       double *alpha, 
			       double *tz, 
			       int    *pivot, 
			       double *w, 
			       double *wp, 
			       double *z0, 
			       double *z1, 
			       double *sspar, 
			       double *par, 
			       int    *ipar,

			       int    tweak) /* added to adjust step size */
{
    /* System generated locals */
    int qr_dim1, qr_offset, i__1, i__2;
    double d__1 /* ,abx UNUSED */;

    /* Builtin functions */
    /*    double sqrt(double); CANT DECLARE BUILTINS UNDER C++ */
    integer s_wsfe(), do_fio(), e_wsfe(); /* in Hom_params.c these are int's */

    /* Local variables */
    static int nfec;
    static double hold;
    static int iter;
    extern double dnrm2_(integer    *n, 
			 doublereal *dx, 
			 integer    *incx);
    static double h, s;
    static long int crash;
    static int limit;
    extern double d1mach_(integer *);
    static long int start;
    static int nc, iflagc, jw;
    static double abserr, relerr;
    extern /* Subroutine */ int stepnf_(integer    *n, 
					integer    *nfe, 
					integer    *iflag, 
					logical    *start, 
					logical    *crash, 
					doublereal *hold, 
					doublereal *h, 
					doublereal *relerr, 
					doublereal *abserr, 
					doublereal *s, 
					doublereal *y, 
					doublereal *yp, 
					doublereal *yold, 
					doublereal *ypold, 
					doublereal *a, 
					doublereal *qr, 
					doublereal *alpha, 
					doublereal *tz, 
					integer    *pivot, 
					doublereal *w, 
					doublereal *wp, 
					doublereal *z0, 
					doublereal *z1,
					doublereal *sspar, 
					doublereal *par, 
					integer    *ipar);

    static double curtol;

    extern /* 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);
    static int np1;
    static long int polsys;



/* SUBROUTINE  FIXPNF  FINDS A FIXED POINT OR ZERO OF THE */
/* N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE */
/* OF A GENERAL HOMOTOPY MAP RHO(A,LAMBDA,X).  FOR THE FIXED */
/* POINT PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL */
/* INTO ITSELF.  THE EQUATION  X = F(X)  IS SOLVED BY */
/* FOLLOWING THE ZERO CURVE OF THE HOMOTOPY MAP */

/*  LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A)  , */

/* STARTING FROM LAMBDA = 0, X = A.  THE CURVE IS PARAMETERIZED */
/* BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY */
/* DIFFERENTIAL EQUATION  D(HOMOTOPY MAP)/DS = 0  FOR */
/* Y(S) = (LAMBDA(S), X(S)) USING A HERMITE CUBIC PREDICTOR AND A */
/* CORRECTOR WHICH RETURNS TO THE ZERO CURVE ALONG THE FLOW NORMAL */
/* TO THE DAVIDENKO FLOW (WHICH CONSISTS OF THE INTEGRAL CURVES OF */
/* D(HOMOTOPY MAP)/DS ). */

/* FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP */
/* SUCH THAT FOR SOME R > 0,  X*F(X) >= 0  WHENEVER NORM(X) = R. */
/* THE EQUATION  F(X) = 0  IS SOLVED BY FOLLOWING THE ZERO CURVE */
/* OF THE HOMOTOPY MAP */

/*   LAMBDA*F(X) + (1 - LAMBDA)*(X - A) */

/* EMANATING FROM LAMBDA = 0, X = A. */

/*  A  MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS. */

/* FOR THE CURVE TRACKING PROBLEM RHO(A,LAMBDA,X) IS ASSUMED TO */
/* BE A C2 MAP FROM E**M X [0,1) X E**N INTO E**N, WHICH FOR */
/* ALMOST ALL PARAMETER VECTORS A IN SOME NONEMPTY OPEN SUBSET */
/* OF E**M SATISFIES */

/*  RANK [D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX] = N */

/* FOR ALL POINTS (LAMBDA,X) SUCH THAT RHO(A,LAMBDA,X)=0.  IT IS */
/* FURTHER ASSUMED THAT */

/*           RANK [ D RHO(A,0,X0)/DX ] = N  . */

/* WITH A FIXED, THE ZERO CURVE OF RHO(A,LAMBDA,X) EMANATING */
/* FROM  LAMBDA = 0, X = X0  IS TRACKED UNTIL  LAMBDA = 1  BY */
/* SOLVING THE ORDINARY DIFFERENTIAL EQUATION */
/* D RHO(A,LAMBDA(S),X(S))/DS = 0  FOR  Y(S) = (LAMBDA(S), X(S)), */
/* WHERE S IS ARC LENGTH ALONG THE ZERO CURVE.  ALSO THE HOMOTOPY */
/* MAP RHO(A,LAMBDA,X) IS ASSUMED TO BE CONSTRUCTED SUCH THAT */

/*              D LAMBDA(0)/DS > 0  . */


/* FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER MUST SUPPLY */
/* A SUBROUTINE  F(X,V)  WHICH EVALUATES F(X) AT X AND RETURNS THE */
/* VECTOR F(X) IN V, AND A SUBROUTINE  FJAC(X,V,K)  WHICH RETURNS IN V */
/* THE KTH COLUMN OF THE JACOBIAN MATRIX OF F(X) EVALUATED AT X.  FOR */
/* THE CURVE TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE */
/*  RHO(A,LAMBDA,X,V,PAR,IPAR)  WHICH EVALUATES THE HOMOTOPY MAP RHO AT */
/* (A,LAMBDA,X) AND RETURNS THE VECTOR RHO(A,LAMBDA,X) IN V, AND A */
/* SUBROUTINE  RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR)  WHICH RETURNS IN V THE KTH
 */
/* COLUMN OF THE N X (N+1) JACOBIAN MATRIX [D RHO/D LAMBDA, D RHO/DX] */
/* EVALUATED AT (A,LAMBDA,X).  FIXPNF  DIRECTLY OR INDIRECTLY USES */
/* THE SUBROUTINES  STEPNF , TANGNF , ROOTNF , ROOT , F (OR  RHO ), */
/* FJAC (OR  RHOJAC ), D1MACH , AND THE BLAS FUNCTIONS  DDOT  AND */
/* DNRM2 .  ONLY  D1MACH  CONTAINS MACHINE DEPENDENT CONSTANTS. */
/* NO OTHER MODIFICATIONS BY THE USER ARE REQUIRED. */


/* ON INPUT: */

/* N  IS THE DIMENSION OF X, F(X), AND RHO(A,LAMBDA,X). */

/* Y  IS AN ARRRAY OF LENGTH  N + 1.  (Y(2),...,Y(N+1)) = A  IS THE */
/*    STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND */
/*    ZERO FINDING PROBLEMS.  (Y(2),...,Y(N+1)) = X0  FOR THE CURVE */
/*    TRACKING PROBLEM. */

/* IFLAG  CAN BE -2, -1, 0, 2, OR 3.  IFLAG  SHOULD BE 0 ON THE */
/*    FIRST CALL TO  FIXPNF  FOR THE PROBLEM  X=F(X), -1 FOR THE */
/*    PROBLEM  F(X)=0, AND -2 FOR THE PROBLEM  RHO(A,LAMBDA,X)=0. */
/*    IN CERTAIN SITUATIONS  IFLAG  IS SET TO 2 OR 3 BY  FIXPNF, */
/*    AND  FIXPNF  CAN BE CALLED AGAIN WITHOUT CHANGING  IFLAG. */

/* ARCRE , ARCAE  ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY, */
/*    ALLOWED THE NORMAL FLOW ITERATION ALONG THE ZERO CURVE.  IF */
/*    ARC?E .LE. 0.0  ON INPUT IT IS RESET TO  .5*SQRT(ANS?E) . */
/*    NORMALLY  ARC?E SHOULD BE CONSIDERABLY LARGER THAN  ANS?E . */

/* ANSRE , ANSAE  ARE THE RELATIVE AND ABSOLUTE ERROR VALUES USED FOR */
/*    THE ANSWER AT LAMBDA = 1.  THE ACCEPTED ANSWER  Y = (LAMBDA, X) */
/*    SATISFIES */

/*       |Y(1) - 1|  .LE.  ANSRE + ANSAE           .AND. */

/*       ||Z||  .LE.  ANSRE*||X|| + ANSAE          WHERE */

/*    (.,Z) IS THE NEWTON STEP TO Y. */

/* TRACE  IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR */
/*    INTERMEDIATE OUTPUT.  IF  TRACE .GT. 0  THE POINTS COMPUTED ON */
/*    THE ZERO CURVE ARE WRITTEN TO I/O UNIT  TRACE . */

/* A(1:*)  CONTAINS THE PARAMETER VECTOR  A .  FOR THE FIXED POINT */
/*    AND ZERO FINDING PROBLEMS, A  NEED NOT BE INITIALIZED BY THE */
/*    USER, AND IS ASSUMED TO HAVE LENGTH  N.  FOR THE CURVE */
/*    TRACKING PROBLEM, A  MUST BE INITIALIZED BY THE USER. */

/* YP(1:N+1)  IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO */
/*    THE ZERO CURVE AT THE CURRENT POINT  Y . */

/* YOLD(1:N+1)  IS A WORK ARRAY CONTAINING THE PREVIOUS POINT FOUND */
/*    ON THE ZERO CURVE. */

/* YPOLD(1:N+1)  IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO */
/*    THE ZERO CURVE AT  YOLD . */

/* 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) , Z0(1:N+1) , Z1(1:N+1)  ARE ALL WORK ARRAYS USED BY */
/*    STEPNF  TO CALCULATE THE TANGENT VECTORS AND NEWTON STEPS. */

/* SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P)  IS */
/*    A VECTOR OF PARAMETERS USED FOR THE OPTIMAL STEP SIZE ESTIMATION. */
/*    IF  SSPAR(J) .LE. 0.0  ON INPUT, IT IS RESET TO A DEFAULT VALUE */
/*    BY  FIXPNF .  OTHERWISE THE INPUT VALUE OF  SSPAR(J)  IS USED. */
/*    SEE THE COMMENTS BELOW AND IN  STEPNF  FOR MORE INFORMATION ABOUT */
/*    THESE CONSTANTS. */

/* PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, */
/*    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES */
/*    RHO, RHOJAC. */


/* ON OUTPUT: */

/* N , TRACE , A  ARE UNCHANGED. */

/* Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = X, AND Y IS AN APPROXIMATE */
/*    ZERO OF THE HOMOTOPY MAP.  NORMALLY LAMBDA = 1 AND X IS A */
/*    FIXED POINT(ZERO) OF F(X).  IN ABNORMAL SITUATIONS LAMBDA */
/*    MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO). */

/* IFLAG = */
/*  -2   CAUSES  FIXPNF  TO INITIALIZE EVERYTHING FOR THE PROBLEM */
/*       RHO(A,LAMBDA,X) = 0 (USE ON FIRST CALL). */

/*  -1   CAUSES  FIXPNF  TO INITIALIZE EVERYTHING FOR THE PROBLEM */
/*       F(X) = 0 (USE ON FIRST CALL). */

/*   0   CAUSES  FIXPNF  TO INITIALIZE EVERYTHING FOR THE PROBLEM */
/*       X = F(X) (USE ON FIRST CALL). */

/*   1   NORMAL RETURN. */

/*   2   SPECIFIED ERROR TOLERANCE CANNOT BE MET.  SOME OR ALL OF */
/*       ARCRE , ARCAE , ANSRE , ANSAE  HAVE BEEN INCREASED TO */
/*       SUITABLE VALUES.  TO CONTINUE, JUST CALL  FIXPNF  AGAIN */
/*       WITHOUT CHANGING ANY PARAMETERS. */

/*   3   STEPNF  HAS BEEN CALLED 1000 TIMES.  TO CONTINUE, CALL */
/*       FIXPNF  AGAIN WITHOUT CHANGING ANY PARAMETERS. */

/*   4   JACOBIAN MATRIX DOES NOT HAVE FULL RANK.  THE ALGORITHM */
/*       HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE */
/*       FOLLOWED ANY FURTHER). */

/*   5   THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE */
/*       HOMOTOPY MAP AND IS NOT MAKING PROGRESS.  THE ERROR TOLERANCES */
/*       ARC?E  AND  ANS?E  WERE TOO LENIENT.  THE PROBLEM SHOULD BE */
/*       RESTARTED BY CALLING  FIXPNF  WITH SMALLER ERROR TOLERANCES */
/*       AND  IFLAG = 0 (-1, -2). */

/*   6   THE NORMAL FLOW NEWTON ITERATION IN  STEPNF  OR  ROOTNF */
/*       FAILED TO CONVERGE.  THE ERROR TOLERANCES  ANS?E  MAY BE TOO */
/*       STRINGENT. */

/*   7   ILLEGAL INPUT PARAMETERS, A FATAL ERROR. */

/* ARCRE , ARCAE , ANSRE , ANSAE  ARE UNCHANGED AFTER A NORMAL RETURN */
/*    (IFLAG = 1).  THEY ARE INCREASED TO APPROPRIATE VALUES ON THE */
/*    RETURN  IFLAG = 2 . */

/* NFE  IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF */
/*    JACOBIAN EVALUATIONS). */

/* ARCLEN  IS THE LENGTH OF THE PATH FOLLOWED. */




/* ***** ARRAY DECLARATIONS. ***** */


/* ***** END OF DIMENSIONAL INFORMATION. ***** */


/* LIMITD  IS AN UPPER BOUND ON THE NUMBER OF STEPS.  IT MAY BE */
/* CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT: */

/* SWITCH FROM THE TOLERANCE  ARC?E  TO THE (FINER) TOLERANCE  ANS?E  IF 
*/
/* THE CURVATURE OF ANY COMPONENT OF  Y  EXCEEDS  CURSW. */



/* :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  :  : 
*/
/* SET LOGICAL SWITCH TO REFLECT ENTRY POINT. */
    /* Parameter adjustments */
    --y;
    --a;
    --yp;
    --yold;
    --ypold;
    qr_dim1 = *n;
    qr_offset = qr_dim1 + 1;
    qr -= qr_offset;
    --alpha;
    --tz;
    --pivot;
    --w;
    --wp;
    --z0;
    --z1;
    --sspar;
    --par;
    --ipar;

    /* Function Body */
    switch(n__) {
	case 1: goto L_polynf;
	}

    polsys = FALSE_;
    goto L11;

L_polynf:
    polsys = TRUE_;
L11:

    if (*n <= 0 || *ansre <= (float)0. || *ansae < (float)0.) {
	*iflag = 7;
    }
    if (*iflag >= -2 && *iflag <= 0) {
	goto L20;
    }
    if (*iflag == 2) {
	goto L120;
    }
    if (*iflag == 3) {
	goto L90;
    }
/* ONLY VALID INPUT FOR  IFLAG  IS -2, -1, 0, 2, 3. */
    *iflag = 7;
    return 0;

/* *****  INITIALIZATION BLOCK.  ***** */

L20:
    *arclen = (float)0.;
    if (*arcre <= (float)0.) {
	*arcre = sqrt(*ansre) * (float).5;
    }
    if (*arcae <= (float)0.) {
	*arcae = sqrt(*ansae) * (float).5;
    }
    nc = *n;
    nfec = 0;
    iflagc = *iflag;
    np1 = *n + 1;
/* SET INITIAL CONDITIONS FOR FIRST CALL TO  STEPNF . */
    start = TRUE_;
    crash = FALSE_;
    hold = (float)1.;
    h = (float).1;
    s = (float)0.;
    ypold[1] = (float)1.;
    yp[1] = (float)1.;
    y[1] = (float)0.;
    i__1 = np1;
    for (jw = 2; jw <= i__1; ++jw) {
	ypold[jw] = (float)0.;
	yp[jw] = (float)0.;
/* L40: */
    }
/* SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS. */
/* LET Z[K] DENOTE THE NEWTON ITERATES ALONG THE FLOW NORMAL TO THE */
/* DAVIDENKO FLOW AND Y THEIR LIMIT. */
/* IDEAL CONTRACTION FACTOR:  ||Z[2] - Z[1]|| / ||Z[1] - Z[0]|| */
    if (sspar[1] <= (float)0.) {
	sspar[1] = (float).5;
    }
/* IDEAL RESIDUAL FACTOR:  ||RHO(A, Z[1])|| / ||RHO(A, Z[0])|| */
    if (sspar[2] <= (float)0.) {
	sspar[2] = (float).01;
    }
/* IDEAL DISTANCE FACTOR:  ||Z[1] - Y|| / ||Z[0] - Y|| */
    if (sspar[3] <= (float)0.) {
	sspar[3] = (float).5;
    }
/* MINIMUM STEP SIZE  HMIN . */
    if (sspar[4] <= (float)0.) {
	sspar[4] = (sqrt(*n + (float)1.) + (float)4.) * d1mach_(&c__4);
    }
/* MAXIMUM STEP SIZE  HMAX . */
    if (sspar[5] <= (float)0.) {
	sspar[5] = (float)1.;
    }
/* MINIMUM STEP SIZE REDUCTION FACTOR  BMIN . */
    if (sspar[6] <= (float)0.) {
	sspar[6] = (float).1;
    }
/* MAXIMUM STEP SIZE EXPANSION FACTOR  BMAX . */
    if (sspar[7] <= (float)0.) {
	sspar[7] = (float)3.;
    }
/* ASSUMED OPERATING ORDER  P . */
    if (sspar[8] <= (float)0.) {
	sspar[8] = (float)2.;
    }

    /* Adjustments of step size */
    if (tweak == 1)
      sspar[5] = (float).1; 
    else if (tweak == 2)
      sspar[5] = (float).01;
    else if (tweak == 3)
      sspar[5] = (float).001; 


/* LOAD  A  FOR THE FIXED POINT AND ZERO FINDING PROBLEMS. */
    if (iflagc >= -1) {
	i__1 = np1;
	for (jw = 2; jw <= i__1; ++jw) {
	    a[jw - 1] = y[jw];
/* L60: */
	}
    }
L90:
    limit = 1000;