pelhomot.cc

Gambit 是一个游戏库理论软件

	alpha[k] = alphak;
	beta = (float)1. / (sigma - qrkk * alphak);
	qr[k + k * qr_dim1] = qrkk - alphak;
	i__2 = np1;
	for (j = kp1; j <= i__2; ++j) {
/* L80: */
	    i__3 = *n - k + 1;
	    y[j] = beta * ddot_(&i__3, &qr[k + k * qr_dim1], &c__1, &qr[k + j 
		    * qr_dim1], &c__1);
	}
	i__3 = np1;
	for (j = kp1; j <= i__3; ++j) {
	    i__2 = *n;
	    for (i = k; i <= i__2; ++i) {
		qr[i + j * qr_dim1] -= qr[i + k * qr_dim1] * y[j];
/* L90: */
	    }
/* Computing 2nd power */
	    d__1 = qr[k + j * qr_dim1];
	    sum[j] -= d__1 * d__1;
/* L100: */
	}
L500:
	;
    }
    alpha[*n] = qr[*n + *n * qr_dim1];
    return 0;
} /* dcpose_ */

/* end dcpose.c */

/**********************************************************************/
/****************** implementation code from ddot.c *******************/
/**********************************************************************/

doublereal ddot_(integer    *n,
		 doublereal *dx,
		 integer    *incx,
		 doublereal *dy,
		 integer    *incy)
{
    /* System generated locals */
    integer i__1;
    doublereal ret_val;

    /* Local variables */
    static integer i, m;
    static doublereal dtemp;
    static integer ix, iy, mp1;


/*     forms the dot product of two vectors. */
/*     uses unrolled loops for increments equal to one. */
/*     jack dongarra, linpack, 3/11/78. */


    /* Parameter adjustments */
    --dy;
    --dx;

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

/*        code for unequal increments or equal increments */
/*          not equal to 1 */

    ix = 1;
    iy = 1;
    if (*incx < 0) {
	ix = (-(*n) + 1) * *incx + 1;
    }
    if (*incy < 0) {
	iy = (-(*n) + 1) * *incy + 1;
    }
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	dtemp += dx[ix] * dy[iy];
	ix += *incx;
	iy += *incy;
/* L10: */
    }
    ret_val = dtemp;
    return ret_val;

/*        code for both increments equal to 1 */


/*        clean-up loop */

L20:
    m = *n % 5;
    if (m == 0) {
	goto L40;
    }
    i__1 = m;
    for (i = 1; i <= i__1; ++i) {
	dtemp += dx[i] * dy[i];
/* L30: */
    }
    if (*n < 5) {
	goto L60;
    }
L40:
    mp1 = m + 1;
    i__1 = *n;
    for (i = mp1; i <= i__1; i += 5) {
	dtemp = dtemp + dx[i] * dy[i] + dx[i + 1] * dy[i + 1] + dx[i + 2] * 
		dy[i + 2] + dx[i + 3] * dy[i + 3] + dx[i + 4] * dy[i + 4];
/* L50: */
    }
L60:
    ret_val = dtemp;
    return ret_val;
} /* ddot_ */

/* end ddot.c */

/**********************************************************************/
/******************* implementation code from divp.c *********************/
/**********************************************************************/

/* Table of constant values */

static integer c__2 = 2;

/* Subroutine */ int divp_(doublereal *xxxx, 
			   doublereal *yyyy, 
			   doublereal *zzzz, 
			   integer    *ierr)
{
    static doublereal xnum, denom;
    extern doublereal d1mach_(integer *i);


/* THIS SUBROUTINE PERFORMS DIVISION  OF COMPLEX NUMBERS: */
/* ZZZZ = XXXX/YYYY */

/* ON INPUT: */

/* XXXX  IS AN ARRAY OF LENGTH TWO REPRESENTING THE FIRST COMPLEX */
/*       NUMBER, WHERE XXXX(1) = REAL PART OF XXXX AND XXXX(2) = */
/*       IMAGINARY PART OF XXXX. */

/* YYYY  IS AN ARRAY OF LENGTH TWO REPRESENTING THE SECOND COMPLEX */
/*       NUMBER, WHERE YYYY(1) = REAL PART OF YYYY AND YYYY(2) = */
/*       IMAGINARY PART OF YYYY. */

/* ON OUTPUT: */

/* ZZZZ  IS AN ARRAY OF LENGTH TWO REPRESENTING THE RESULT OF */
/*       THE DIVISION, ZZZZ = XXXX/YYYY, WHERE ZZZZ(1) = */
/*       REAL PART OF ZZZZ AND ZZZZ(2) = IMAGINARY PART OF ZZZZ. */

/* IERR = */
/*  1   IF DIVISION WOULD HAVE CAUSED OVERFLOW.  IN THIS CASE, THE */
/*      APPROPRIATE PARTS OF ZZZZ ARE SET EQUAL TO THE LARGEST */
/*      FLOATING POINT NUMBER, AS GIVEN BY FUNCTION  D1MACH . */

/*  0   IF DIVISION DOES NOT CAUSE OVERFLOW. */

/* DECLARATION OF INPUT */

/* DECLARATION OF OUTPUT */

/* DECLARATION OF VARIABLES */

    /* Parameter adjustments */
    --zzzz;
    --yyyy;
    --xxxx;

    /* Function Body */
    *ierr = 0;
    denom = yyyy[1] * yyyy[1] + yyyy[2] * yyyy[2];
    xnum = xxxx[1] * yyyy[1] + xxxx[2] * yyyy[2];
    if (dblabs(denom) >= (float)1. || (dblabs(denom) < (float)1. && dblabs(xnum) / 
	    d1mach_(&c__2) < dblabs(denom)) ) {
	zzzz[1] = xnum / denom;
    } else {
	zzzz[1] = d1mach_(&c__2);
	*ierr = 1;
    }
    xnum = xxxx[2] * yyyy[1] - xxxx[1] * yyyy[2];
    if (dblabs(denom) >= (float)1. || (dblabs(denom) < (float)1. && dblabs(xnum) / 
	    d1mach_(&c__2) < dblabs(denom))) {
	zzzz[2] = xnum / denom;
    } else {
	zzzz[2] = d1mach_(&c__2);
	*ierr = 1;
    }
    return 0;
} /* divp_ */

/* end divp.c */

/**********************************************************************/
/****************** implementation code from dnrm2.c ******************/
/**********************************************************************/

doublereal dnrm2_(integer    *n, 
		  doublereal *dx, 
		  integer    *incx)
{
    /* Initialized data */

    static doublereal zero = 0.;
    static doublereal one = 1.;
    static doublereal cutlo = 8.232e-11;
    static doublereal cuthi = 1.304e19;

    /* Format strings */
    static char fmt_30[] = "";
    static char fmt_50[] = "";
    static char fmt_70[] = "";
    static char fmt_110[] = "";

    /* System generated locals */
    integer i__1;
    doublereal ret_val, d__1;

    /* Builtin functions */
    /*    double sqrt(double); CANT DECLARE BUILTINS UNDER C++ */

    /* Local variables */
    static doublereal xmax;
    static integer next, i, j, ix;
    static doublereal hitest, sum;

    /* Assigned format variables */
    char *next_fmt;

    /* Parameter adjustments */
    --dx;

    /* Function Body */

/*     euclidean norm of the n-vector stored in dx() with storage */
/*     increment incx . */
/*     if    n .le. 0 return with result = 0. */
/*     if n .ge. 1 then incx must be .ge. 1 */

/*           c.l.lawson, 1978 jan 08 */
/*     modified to correct failure to update ix, 1/25/92. */
/*     modified 3/93 to return if incx .le. 0. */

/*     four phase method     using two built-in constants that are */
/*     hopefully applicable to all machines. */
/*         cutlo = maximum of  dsqrt(u/eps)  over all known machines. */
/*         cuthi = minimum of  dsqrt(v)      over all known machines. */
/*     where */
/*         eps = smallest no. such that eps + 1. .gt. 1. */
/*         u   = smallest positive no.   (underflow limit) */
/*         v   = largest  no.            (overflow  limit) */

/*     brief outline of algorithm.. */

/*     phase 1    scans zero components. */
/*     move to phase 2 when a component is nonzero and .le. cutlo */
/*     move to phase 3 when a component is .gt. cutlo */
/*     move to phase 4 when a component is .ge. cuthi/m */
/*     where m = n for x() real and m = 2*n for complex. */

/*     values for cutlo and cuthi.. */
/*     from the environmental parameters listed in the imsl converter */
/*     document the limiting values are as follows.. */
/*     cutlo, s.p.   u/eps = 2**(-102) for  honeywell.  close seconds are 
*/
/*                   univac and dec at 2**(-103) */
/*                   thus cutlo = 2**(-51) = 4.44089e-16 */
/*     cuthi, s.p.   v = 2**127 for univac, honeywell, and dec. */
/*                   thus cuthi = 2**(63.5) = 1.30438e19 */
/*     cutlo, d.p.   u/eps = 2**(-67) for honeywell and dec. */
/*                   thus cutlo = 2**(-33.5) = 8.23181d-11 */
/*     cuthi, d.p.   same as s.p.  cuthi = 1.30438d19 */
/*     data cutlo, cuthi / 8.232d-11,  1.304d19 / */
/*     data cutlo, cuthi / 4.441e-16,  1.304e19 / */

    if (*n > 0 && *incx > 0) {
	goto L10;
    }
    ret_val = zero;
    goto L300;

L10:
    next = 0;
    next_fmt = fmt_30;
    sum = zero;
    i = 1;
    ix = 1;
/*                                                 begin main loop */
L20:
    switch ((int)next) {
	case 0: goto L30;
	case 1: goto L50;
	case 2: goto L70;
	case 3: goto L110;
    }
L30:
    if ((d__1 = dx[i], dblabs(d__1)) > cutlo) {
	goto L85;
    }
    next = 1;
    next_fmt = fmt_50;
    xmax = zero;

/*                        phase 1.  sum is zero */

L50:
    if (dx[i] == zero) {
	goto L200;
    }
    if ((d__1 = dx[i], dblabs(d__1)) > cutlo) {
	goto L85;
    }

/*                                prepare for phase 2. */
    next = 2;
    next_fmt = fmt_70;
    goto L105;

/*                                prepare for phase 4. */

L100:
    ix = j;
    next = 3;
    next_fmt = fmt_110;
    sum = sum / dx[i] / dx[i];
L105:
    xmax = (d__1 = dx[i], dblabs(d__1));
    goto L115;

/*                   phase 2.  sum is small. */
/*                             scale to avoid destructive underflow. */

L70:
    if ((d__1 = dx[i], dblabs(d__1)) > cutlo) {
	goto L75;
    }

/*                     common code for phases 2 and 4. */
/*                     in phase 4 sum is large.  scale to avoid overflow. 
*/

L110:
    if ((d__1 = dx[i], dblabs(d__1)) <= xmax) {
	goto L115;
    }
/* Computing 2nd power */
    d__1 = xmax / dx[i];
    sum = one + sum * (d__1 * d__1);
    xmax = (d__1 = dx[i], dblabs(d__1));
    goto L200;

L115:
/* Computing 2nd power */
    d__1 = dx[i] / xmax;
    sum += d__1 * d__1;
    goto L200;


/*                  prepare for phase 3. */

L75:
    sum = sum * xmax * xmax;


/*     for real or d.p. set hitest = cuthi/n */
/*     for complex      set hitest = cuthi/(2*n) */

L85:
    hitest = cuthi / (real) (*n);

/*                   phase 3.  sum is mid-range.  no scaling. */

    i__1 = *n;
    for (j = ix; j <= i__1; ++j) {
	if ((d__1 = dx[i], dblabs(d__1)) >= hitest) {
	    goto L100;
	}
/* Computing 2nd power */
	d__1 = dx[i];
	sum += d__1 * d__1;
	i += *incx;
/* L95: */
    }
    ret_val = sqrt(sum);
    goto L300;

L200:
    ++ix;
    i += *incx;
    if (ix <= *n) {
	goto L20;
    }

/*              end of main loop. */

/*              compute square root and adjust for scaling. */

    ret_val = xmax * sqrt(sum);
L300:
    return ret_val;
} /* dnrm2_ */

/* end dnrm2.c */

/**********************************************************************/
/***************** implementation code from dscal.c *******************/
/**********************************************************************/

/* Subroutine */ int dscal_(integer    *n, 
			    doublereal *da, 
			    doublereal *dx, 
			    integer    *incx)
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    static integer i, m, nincx, mp1;


/*     scales a vector by a constant. */
/*     uses unrolled loops for increment equal to one. */
/*     jack dongarra, linpack, 3/11/78. */
/*     modified 3/93 to return if incx .le. 0. */


    /* Parameter adjustments */
    --dx;

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

/*        code for increment not equal to 1 */

    nincx = *n * *incx;
    i__1 = nincx;
    i__2 = *incx;
    for (i = 1; i__2 < 0 ? i >= i__1 : i <= i__1; i += i__2) {
	dx[i] = *da * dx[i];
/* L10: */
    }
    return 0;

/*        code for increment equal to 1 */


/*        clean-up loop */

L20:
    m = *n % 5;
    if (m == 0) {
	goto L40;
    }
    i__2 = m;
    for (i = 1; i <= i__2; ++i) {
	dx[i] = *da * dx[i];
/* L30: */
    }
    if (*n < 5) {
	return 0;
    }
L40:
    mp1 = m + 1;
    i__2 = *n;
    for (i = mp1; i <= i__2; i += 5) {
	dx[i] = *da * dx[i];
	dx[i + 1] = *da * dx[i + 1];