/* getnp.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"

/* Copyright (C) 2000, International Business Machines */
/* Corporation and others.  All Rights Reserved. */



/* Subroutine */ int getnp_(ipoly, poly, lpoly, n, kappa, g, h__)
integer *ipoly;
doublereal *poly;
integer *lpoly, *n;
doublereal *kappa, *g, *h__;
{
    /* System generated locals */
    integer h_dim1, h_offset, i__1, i__2;

    /* Local variables */
    static integer kbeg, i__, j, k, dd, np1;


/*     collects  the coefficient of the polynomial IPOLY into */
/*     doefficient of a quadratic kappa+g'x+0.5 x'H'x */
/*     INPUT */
/*     POLY   = Newton polynomials */
/*     IPOLY  = The index of the polynomial we are processing */
/*     OUTPUT */
/*     KAPPA = the constant of the polynomail */
/*     G     = the linear coeffitients */
/*     H     = the n by n matrix of quadratic doeffitients */

/*     Local variables */

    /* Parameter adjustments */
    --poly;
    h_dim1 = *n;
    h_offset = 1 + h_dim1 * 1;
    h__ -= h_offset;
    --g;

    /* Function Body */
    np1 = *n + 1;
    dd = np1 * (*n + 2) / 2;
    *kappa = 0.;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	g[i__] = 0.;
	i__2 = *n;
	for (j = 1; j <= i__2; ++j) {
	    h__[i__ + j * h_dim1] = 0.;
/* L20: */
	}
/* L10: */
    }
    if (*ipoly == 1) {
	*kappa = poly[1];
    } else if (*ipoly <= np1) {
	kbeg = (*ipoly - 2) * np1 + 2;
	*kappa = poly[kbeg];
	i__1 = *n;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    g[i__] = poly[kbeg + i__];
/* L30: */
	}
    } else {
	kbeg = *n * np1 + 2 + (*ipoly - np1 - 1) * dd;
	*kappa = poly[kbeg];
	i__1 = *n;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    g[i__] = poly[kbeg + i__];
/* L40: */
	}
	k = kbeg + *n + 1;
	i__1 = *n;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    h__[i__ + i__ * h_dim1] = poly[k] * 2;
	    ++k;
	    i__2 = *n;
	    for (j = i__ + 1; j <= i__2; ++j) {
		h__[i__ + j * h_dim1] = poly[k];
		h__[j + i__ * h_dim1] = poly[k];
		++k;
/* L60: */
	    }
/* L50: */
	}
    }
    return 0;
} /* getnp_ */