/* ptexch.f -- translated by f2c (version 20000817).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"
#include "blaswrap.h"

/* Common Block Declarations */

struct {
    integer iout, iprint;
    doublereal mcheps, cnstol;
} dfocm_;

#define dfocm_1 dfocm_

/* Copyright (C) 2000, International Business Machines */
/* Corporation and others.  All Rights Reserved. */
/* Subroutine */ int ptexch_(poly, xnew, pntint, ipoly, nind, n, lpoly, 
	lptint, pivthr, pivval, vmax, fail)
doublereal *poly, *xnew, *pntint;
integer *ipoly, *nind, *n, *lpoly, *lptint;
doublereal *pivthr, *pivval, *vmax;
logical *fail;
{
    /* System generated locals */
    integer i__1;
    doublereal d__1;

    /* Local variables */
    static integer base, i__, block;
    extern /* Subroutine */ int evalx_();
    static integer dd;
    extern /* Subroutine */ int ptrepl_();
    static integer np1, np2;
    static doublereal val;


/*  ********************************************************************** */
/*  THIS SUBROUTINE INCLUDES A POINT XNEW IN THE INTERPOLATION SET BY */
/*  UPDATING THE SET OF NEWTON FUNDAMENTAL POLYNOMIALS ACCORDINGLY. */
/*  FIRST A POINT WHICH IS  REPLACED BY XNEW IS DETERMINED, THEN */
/*  SUBROUTINE "PTREPL" IS CALLED, WHICH PERFORMS THE UPDATES. */
/*  IPOLY IS THE INDEX OF THE POINT THAT IS REPLACED BY 'XGNEW'. */
/*  THE SET OF NEWTON POLYNOMIALS IS UPDATED ACCORDINGLY. */

/*  IF THERE ARE QUADRATIC POLYNOMIALS IN THE INTERPOLATION BASIS THEN */
/*  WE EVALUATE EACH POLYNOMIAL IN THE QUADRATIC BLOCK AT 'XNEW', */
/*  OTHERWISE WE EVALUATE POLYNOMIALS OF THE LINEAR BLOCK AT 'XNEW'. */
/*  WE CHOOSE THE POLYNOMIAL  THAT HAS THE LARGEST VALUE AT XNEW. */
/*  IF THIS VALUE IS LARGER THEN THE PIVOT THRESHOLD THEN WE SET */
/*  TO THE INDEX OF THAT POLYNOMIAL. */
/*  OTHERWISE  WE DECLARE THAT WE FAILED TO INCLUDE XGNEW IN THE */
/*  INTERPOLATION SET AND RETURN. */

/*  PARAMETERS */
/*    POLY   (INPUT/OUTPUT) THE SET OF NEWTON POLYNOMIALS */

/*    XNEW   (INPUT) POINT THAT SHOULD BE ADDED */

/*    PNTINT (INPUT) SET OF INTERPOLATION POINTS BEFORE XNEW IS ADDED */

/*    NIND   (INPUT) NUMBER OF INTERPOLATION POINTS */

/*    N      (INPUT) PROBLEM DIMENSION */

/*    PIVTHR (INPUT) PIVOT THRESHOLD VALUE */

/*    PIVVAL (INPUT) ARRAY OF PIVOT VALUES */

/*    FAIL   (OUTPUT) INDICATES IF THE POINT INTERCHANGE FAILED */

/*    IPOLY  (OUTPUT) THE INDEX OF THE POINT THAT WAS REPLACED, IF FAIL=.FALSE. */

/*    VMAX   (OUTPUT) THE VALUE OF THE PIVOT, IF FAIL=.FALSE. */
/*  ******************************************************************* */


/*  PARAMETER VARIABLES */


/*  COMMON VARIABLES */


/*  LOCAL VARIABLES */


/*  SUBROUTINES AND FUNCTIONS CALLED: */

/*    APPLICATION     :       EVALX, PTREPL */
/*    FORTRAN SUPPLIED:       ABS */

    /* Parameter adjustments */
    --pivval;
    --xnew;
    --poly;
    --pntint;

    /* Function Body */
    *fail = FALSE_;
    base = *ipoly;
    np1 = *n + 1;
    np2 = *n + 2;
    dd = np1 * np2 / 2;

/*  CHECK IF INTERPOLATION SET CONTAINS ELEMENTS FROM QUADRATIC BLOCK OR NOT */

    if (*nind <= np1) {
	block = 1;
    } else {
	block = 2;
    }

/*  TO IDENTIFY IPOLY WE  EVALUATE ALL NEWTON POLYNOMIALS IN THE LAST */
/*  BLOCK AT THIS POINT AND CHOOSE THE ONE THAT GIVES THE LARGE PIVOT. */
/*  IF THIS PIVOT IS LARGER THEN THE PIVOT THRESHOLD THEN WE SET IPOLY TO */
/*  THE INDEX OF THIS POLYNOMIAL. */

    *vmax = 0.;
    *ipoly = 0;
    if (block == 2) {

/*  IF WE START WITH THE QUADRATIC BLOCK */

	i__1 = *nind;
	for (i__ = np2; i__ <= i__1; ++i__) {
	    if (i__ != base) {
		evalx_(&val, &xnew[1], &poly[1], &i__, n, lpoly);
		if (abs(val) > abs(*vmax)) {
		    *vmax = val;
		    *ipoly = i__;
		}
	    }
/* L10: */
	}
	if (*ipoly == 0 || (d__1 = *vmax * pivval[*ipoly], abs(d__1)) < *
		pivthr) {
	    *fail = TRUE_;
	    return 0;
	}
    } else {

/*  IF WE START FROM THE LINEAR BLOCK */

	i__1 = *nind;
	for (i__ = 2; i__ <= i__1; ++i__) {
	    evalx_(&val, &xnew[1], &poly[1], &i__, n, lpoly);
	    if (abs(val) > abs(*vmax)) {
		*vmax = val;
		*ipoly = i__;
	    }
/* L20: */
	}
	if (*ipoly == 0 || (d__1 = *vmax * pivval[*ipoly], abs(d__1)) < *
		pivthr) {
	    *fail = TRUE_;
	    return 0;
	}
    }

/*  UPDATE THE NEWTON POLYNOMIALS, SO THAT POLYNOMIAL WITH THE */
/*  INDEX IPOLY CORRESPONDS TO XGNEW IN THE INTERPOLATION SET */

    ptrepl_(&xnew[1], ipoly, vmax, &pntint[1], &poly[1], nind, n, lpoly, 
	    lptint);
    return 0;
} /* ptexch_ */