/* gterms.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 mterms_(kappa, b, h__, lafla, poly, pntint, nind, n, 
	neqcon, lpoly, lptint)
doublereal *kappa, *b, *h__, *lafla, *poly, *pntint;
integer *nind, *n, *neqcon, *lpoly, *lptint;
{
    /* System generated locals */
    integer h_dim1, h_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer i__, j, k, nindm1, ndmnp1, dd, jj, kk, km;
    extern /* Subroutine */ int rzrvec_();
    static integer np1;
    extern /* Subroutine */ int rzrmat_();
    static integer ndd;


/*  *********************************************************************** */
/*  THIS SUBROUTINE COMPUTES THE TERMS OF THE QUADRATIC */
/*  INTERPOLATING POLYNOMIAL */
/*                     T         T */
/*      M(X)= KAPPA + G X + 0.5*X H X */

/*  USING FINITE DIFFERENCES, FOUND IN 'FG' */

/*  PARAMETERS: */

/*   N      (INPUT)  DIMENTION OF THE PROBLEM */

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

/*   PNTINT (INPUT)  LIST OF  'NIND' DATA POINTS. THE I-TH POINT OCCUPIES */
/*                   POSITIONS ( I - 1 ) * N + 1 TO I * N. */
/*   POLY   (INPUT)  THE ARRAY CONTAINING COEFFICIENTS OF NEWTON FUNDAMENTAL */
/*                   POLYNOMIALS (AS COMPUTED BY 'NBUILD') */
/*   LAFLA  (INPUT)  THE ARRAY OF FINITE DIFFERENCES */

/*   NEQCON (INPUT)  NUMBER OF LINEARLY INDEP. EQUALITY CONSTRAINTS */

/*   KAPPA  (OUTPUT) THE CONSTANT TERM OF THE INTERPOLATION MODEL. */

/*   G      (OUTPUT) VECTOR OF THE LINEAR TERMS OF THE  INTERPOLATION MODEL. */

/*   H      (OUTPUT) MATRIX OF QUADRATIC TERMS OF THE  INTERPOLATION MODEL. */

/*  ************************************************************************** */


/*  LOCAL VARIABLES */

    /* Parameter adjustments */
    --lafla;
    h_dim1 = *n;
    h_offset = 1 + h_dim1 * 1;
    h__ -= h_offset;
    --b;
    --poly;
    --pntint;

    /* Function Body */
    np1 = *n + 1;
    dd = np1 * (*n + 2) / 2;
    ndd = dd - np1;
    nindm1 = *nind - 1;
    ndmnp1 = nindm1 - *n;

/*  SET THE COEFFICIENT OF THE INTERPOLATION TO ZERO */

    rzrmat_(&h__[h_offset], n, n);
    rzrvec_(&b[1], n);

/*  COMPUTE THE CONSTANT TERM */


/*  INITIALIZE USING CONSTANT BLOCK */

    *kappa = lafla[1] * poly[1];

/*  UPDATE USING THE FINITE DIFFERENCES CORRESPONDING TO LINEAR BLOCK */

/* Computing MIN */
    i__1 = *n - *neqcon;
    k = min(i__1,nindm1);
    i__1 = k;
    for (j = 1; j <= i__1; ++j) {
	*kappa += lafla[j + 1] * poly[np1 * (j - 1) + 2];
/* L10: */
    }
    if (*nind > np1) {

/*  UPDATE USING THE FINITE DIFF. CORRESPONDING TO QUADRATIC BLOCK */

	k = min(ndd,ndmnp1);
	i__1 = k;
	for (j = 1; j <= i__1; ++j) {
	    *kappa += lafla[j + np1] * poly[np1 * *n + 2 + dd * (j - 1)];
/* L20: */
	}
    }
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {

/*  COMPUTE DEGREE ONE TERMS */

	if (*nind > 1) {

/*  UPDATE USING LINEAR BLOCK */

/* Computing MIN */
	    i__2 = *n - *neqcon;
	    k = min(i__2,nindm1);
	    i__2 = k;
	    for (j = 1; j <= i__2; ++j) {
		b[i__] += lafla[j + 1] * poly[i__ + 2 + np1 * (j - 1)];
/* L50: */
	    }
	}
	if (*nind > np1) {

/*  UPDATE USING QUADRATIC BLOCK */

	    k = min(ndd,ndmnp1);
	    i__2 = k;
	    for (j = 1; j <= i__2; ++j) {
		b[i__] += lafla[j + np1] * poly[i__ + np1 * *n + 2 + dd * (j 
			- 1)];
/* L60: */
	    }
	}
/* L40: */
    }

/*  COMPUTE QUADRATIC TERMS USING QUADRATIC BLOCK */

    if (*nind > np1) {

/*  POINTER TO FIRST QUADRATIC  BLOCK IN  'POLY' */

	jj = np1 * np1 + 2;

/*  LOOP OVER THE NUMBER OF SECOND DEGREE POLYNOMIALS */

	km = min(ndd,ndmnp1);
	i__1 = km;
	for (j = 1; j <= i__1; ++j) {
	    k = 1;
	    i__2 = *n;
	    for (i__ = 1; i__ <= i__2; ++i__) {

/*  UPDATE DIAGONAL ELEMENT */

		h__[i__ + i__ * h_dim1] += lafla[j + np1] * 2. * poly[jj + k 
			- 1 + dd * (j - 1)];
		++k;

/*  UPDATE OFF-DIAGONAL ELEMENTS */

		i__3 = *n;
		for (kk = i__ + 1; kk <= i__3; ++kk) {
		    h__[i__ + kk * h_dim1] += lafla[j + np1] * poly[jj + k - 
			    1 + dd * (j - 1)];
		    h__[kk + i__ * h_dim1] = h__[i__ + kk * h_dim1];
		    ++k;
/* L91: */
		}
/* L90: */
	    }
/* L80: */
	}
    }
    return 0;
} /* mterms_ */