Propagated the FTN_TRAILING_UNDERSCORE concept to the dpolft_ and
dpcoef_ fortran functions referenced in funcs.cpp . The previous treatment didn't honor the value of FTN_TRAILING_UNDERSCORE. So, if a different source for the SLATEC routines than the one in ext/math was used, this may have caused some conflict. Also, worked up the doxygen documentation for routines in funcs.cpp. This Fix was suggested by Dave Glaze.
This commit is contained in:
Harry Moffat
parent
2ba2271125
commit
fae139765e
4 changed files with 258 additions and 87 deletions
|
|
@ -1,5 +1,16 @@
|
|||
|
||||
// miscellaneous functions
|
||||
/**
|
||||
* @file funcs.cpp file containing miscellaneous
|
||||
* numerical functions.
|
||||
*/
|
||||
/*
|
||||
* $Author$
|
||||
* $Date$
|
||||
* $Revision$
|
||||
*
|
||||
* Copyright 2001-2003 California Institute of Technology
|
||||
* See file License.txt for licensing information
|
||||
*
|
||||
*/
|
||||
|
||||
#ifdef WIN32
|
||||
#pragma warning(disable:4786)
|
||||
|
|
@ -14,66 +25,142 @@ using namespace std;
|
|||
#include "ct_defs.h"
|
||||
#include "ctexceptions.h"
|
||||
#include "stringUtils.h"
|
||||
#include "funcs.h"
|
||||
#include "polyfit.h"
|
||||
|
||||
#ifndef FTN_TRAILING_UNDERSCORE
|
||||
#define _DPOLFT_ dpolft
|
||||
#define _DPCOEF_ dpcoef
|
||||
#else
|
||||
#define _DPOLFT_ dpolft_
|
||||
#define _DPCOEF_ dpcoef_
|
||||
#endif
|
||||
|
||||
extern "C" {
|
||||
int _DPOLFT_(integer* n, doublereal* x, doublereal* y, doublereal* w,
|
||||
integer* maxdeg, integer* ndeg, doublereal* eps, doublereal* r,
|
||||
integer* ierr, doublereal* a);
|
||||
|
||||
int dpolft_(integer* n, doublereal* x, doublereal* y, doublereal* w,
|
||||
integer* maxdeg, integer* ndeg, doublereal* eps, doublereal* r,
|
||||
integer* ierr, doublereal* a);
|
||||
|
||||
int dpcoef_(integer* l, doublereal* c, doublereal* tc, doublereal* a);
|
||||
int _DPCOEF_(integer* l, doublereal* c, doublereal* tc, doublereal* a);
|
||||
}
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
/**
|
||||
* Linearly interpolate a function defined on a discrete grid.
|
||||
* vector xpts contains a monotonic sequence of grid points, and
|
||||
* vector fpts contains function values defined at these points.
|
||||
* The value returned is the linear interpolate at point x.
|
||||
* If x is outside the range of xpts, the value of fpts at the
|
||||
* nearest end is returned.
|
||||
*/
|
||||
// Linearly interpolate a function defined on a discrete grid.
|
||||
/*
|
||||
* Vector xpts contains a monotonic sequence of grid points, and
|
||||
* vector fpts contains function values defined at these points.
|
||||
* The value returned is the linear interpolate at point x.
|
||||
* If x is outside the range of xpts, the value of fpts at the
|
||||
* nearest end is returned.
|
||||
*
|
||||
* @param x value of the x coordinate
|
||||
* @param xpts value of the grid points
|
||||
* @param fpts value of the interpolant at the grid points
|
||||
*
|
||||
* @return Returned value is the value of of the interpolated
|
||||
* function at x.
|
||||
*/
|
||||
doublereal linearInterp(doublereal x, const vector_fp& xpts,
|
||||
const vector_fp& fpts) {
|
||||
if (x <= xpts[0])
|
||||
return fpts[0];
|
||||
if (x >= xpts.back())
|
||||
return fpts.back();
|
||||
vector_fp::const_iterator loc =
|
||||
lower_bound(xpts.begin(), xpts.end(), x);
|
||||
int iloc = int(loc - xpts.begin()) - 1;
|
||||
doublereal ff = fpts[iloc] +
|
||||
(x - xpts[iloc])*(fpts[iloc + 1]
|
||||
- fpts[iloc])/(xpts[iloc + 1] - xpts[iloc]);
|
||||
return ff;
|
||||
}
|
||||
|
||||
doublereal linearInterp(doublereal x, const vector_fp& xpts,
|
||||
const vector_fp& fpts) {
|
||||
if (x <= xpts[0])
|
||||
return fpts[0];
|
||||
if (x >= xpts.back())
|
||||
return fpts.back();
|
||||
vector_fp::const_iterator loc =
|
||||
lower_bound(xpts.begin(), xpts.end(), x);
|
||||
int iloc = int(loc - xpts.begin()) - 1;
|
||||
doublereal ff = fpts[iloc] +
|
||||
(x - xpts[iloc])*(fpts[iloc + 1]
|
||||
- fpts[iloc])/(xpts[iloc + 1] - xpts[iloc]);
|
||||
return ff;
|
||||
}
|
||||
|
||||
|
||||
|
||||
doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w,
|
||||
int maxdeg, int& ndeg, doublereal eps, doublereal* r) {
|
||||
integer nn = n;
|
||||
integer mdeg = maxdeg;
|
||||
integer ndg = ndeg;
|
||||
doublereal epss = eps;
|
||||
integer ierr;
|
||||
int worksize = 3*n + 3*maxdeg + 3;
|
||||
vector_fp awork(worksize,0.0);
|
||||
vector_fp coeffs(n+1, 0.0);
|
||||
doublereal zer = 0.0;
|
||||
|
||||
dpolft_(&nn, x, y, w, &mdeg, &ndg, &epss, &coeffs[0],
|
||||
&ierr, &awork[0]);
|
||||
if (ierr != 1) throw CanteraError("polyfit",
|
||||
"DPOLFT returned error code IERR = " + int2str(ierr) +
|
||||
"while attempting to fit " + int2str(n) + " data points "
|
||||
+ "to a polynomial of degree " + int2str(maxdeg));
|
||||
ndeg = ndg;
|
||||
dpcoef_(&ndg, &zer, r, &awork[0]);
|
||||
return epss;
|
||||
}
|
||||
//! Fits a polynomial function to a set of data points
|
||||
/*!
|
||||
* Given a collection of points X(I) and a set of values Y(I) which
|
||||
* correspond to some function or measurement at each of the X(I),
|
||||
* subroutine DPOLFT computes the weighted least-squares polynomial
|
||||
* fits of all degrees up to some degree either specified by the user
|
||||
* or determined by the routine. The fits thus obtained are in
|
||||
* orthogonal polynomial form. Subroutine DP1VLU may then be
|
||||
* called to evaluate the fitted polynomials and any of their
|
||||
* derivatives at any point. The subroutine DPCOEF may be used to
|
||||
* express the polynomial fits as powers of (X-C) for any specified
|
||||
* point C.
|
||||
*
|
||||
* @param n The number of data points.
|
||||
*
|
||||
* @param x A set of grid points on which the data is specified.
|
||||
* The array of values of the independent variable. These
|
||||
* values may appear in any order and need not all be
|
||||
* distinct. There are n of them.
|
||||
*
|
||||
* @param y array of corresponding function values. There are n of them
|
||||
*
|
||||
* @param w array of positive values to be used as weights. If
|
||||
* W[0] is negative, DPOLFT will set all the weights
|
||||
* to 1.0, which means unweighted least squares error
|
||||
* will be minimized. To minimize relative error, the
|
||||
* user should set the weights to: W(I) = 1.0/Y(I)**2,
|
||||
* I = 1,...,N .
|
||||
*
|
||||
* @param maxdeg maximum degree to be allowed for polynomial fit.
|
||||
* MAXDEG may be any non-negative integer less than N.
|
||||
* Note -- MAXDEG cannot be equal to N-1 when a
|
||||
* statistical test is to be used for degree selection,
|
||||
* i.e., when input value of EPS is negative.
|
||||
*
|
||||
* @param ndeg output degree of the fit computed.
|
||||
*
|
||||
* @param eps Specifies the criterion to be used in determining
|
||||
* the degree of fit to be computed.
|
||||
* (1) If EPS is input negative, DPOLFT chooses the
|
||||
* degree based on a statistical F test of
|
||||
* significance. One of three possible
|
||||
* significance levels will be used: .01, .05 or
|
||||
* .10. If EPS=-1.0 , the routine will
|
||||
* automatically select one of these levels based
|
||||
* on the number of data points and the maximum
|
||||
* degree to be considered. If EPS is input as
|
||||
* -.01, -.05, or -.10, a significance level of
|
||||
* .01, .05, or .10, respectively, will be used.
|
||||
* (2) If EPS is set to 0., DPOLFT computes the
|
||||
* polynomials of degrees 0 through MAXDEG .
|
||||
* (3) If EPS is input positive, EPS is the RMS
|
||||
* error tolerance which must be satisfied by the
|
||||
* fitted polynomial. DPOLFT will increase the
|
||||
* degree of fit until this criterion is met or
|
||||
* until the maximum degree is reached.
|
||||
*
|
||||
* @param r Output vector containing the first LL+1 Taylor coefficients
|
||||
* where LL=ABS(ndeg).
|
||||
* P(X) = r[0] + r[1]*(X-C) + ... + r[ndeg] * (X-C)**ndeg
|
||||
* ( here C = 0.0)
|
||||
*
|
||||
* @return returns the RMS error of the polynomial of degree ndeg .
|
||||
*/
|
||||
doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w,
|
||||
int maxdeg, int& ndeg, doublereal eps, doublereal* r) {
|
||||
integer nn = n;
|
||||
integer mdeg = maxdeg;
|
||||
integer ndg = ndeg;
|
||||
doublereal epss = eps;
|
||||
integer ierr;
|
||||
int worksize = 3*n + 3*maxdeg + 3;
|
||||
vector_fp awork(worksize,0.0);
|
||||
vector_fp coeffs(n+1, 0.0);
|
||||
doublereal zer = 0.0;
|
||||
|
||||
_DPOLFT_(&nn, x, y, w, &mdeg, &ndg, &epss, &coeffs[0],
|
||||
&ierr, &awork[0]);
|
||||
if (ierr != 1) throw CanteraError("polyfit",
|
||||
"DPOLFT returned error code IERR = " + int2str(ierr) +
|
||||
"while attempting to fit " + int2str(n) + " data points "
|
||||
+ "to a polynomial of degree " + int2str(maxdeg));
|
||||
ndeg = ndg;
|
||||
_DPCOEF_(&ndg, &zer, r, &awork[0]);
|
||||
return epss;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,11 +1,40 @@
|
|||
/**
|
||||
* @file funcs.h Header for a file containing miscellaneous
|
||||
* numerical functions.
|
||||
*/
|
||||
/*
|
||||
* $Author$
|
||||
* $Date$
|
||||
* $Revision$
|
||||
*
|
||||
* Copyright 2001-2003 California Institute of Technology
|
||||
* See file License.txt for licensing information
|
||||
*
|
||||
*/
|
||||
#ifndef CT_FUNCS_H
|
||||
#define CT_FUNCS_H
|
||||
|
||||
#include "ct_defs.h"
|
||||
|
||||
namespace Cantera {
|
||||
doublereal linearInterp(doublereal x, const vector_fp& xpts,
|
||||
const vector_fp& fpts);
|
||||
|
||||
//! Linearly interpolate a function defined on a discrete grid.
|
||||
/*!
|
||||
* Vector xpts contains a monotonic sequence of grid points, and
|
||||
* vector fpts contains function values defined at these points.
|
||||
* The value returned is the linear interpolate at point x.
|
||||
* If x is outside the range of xpts, the value of fpts at the
|
||||
* nearest end is returned.
|
||||
*
|
||||
* @param x value of the x coordinate
|
||||
* @param xpts value of the grid points
|
||||
* @param fpts value of the interpolant at the grid points
|
||||
*
|
||||
* @return Returned value is the value of of the interpolated
|
||||
* function at x.
|
||||
*/
|
||||
doublereal linearInterp(doublereal x, const vector_fp& xpts,
|
||||
const vector_fp& fpts);
|
||||
}
|
||||
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -27,35 +27,16 @@
|
|||
extern "C"
|
||||
{
|
||||
|
||||
// /* Subroutine */ int dgelss_(integer *m, integer *n, integer *nrhs,
|
||||
// doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *
|
||||
// s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork,
|
||||
// integer *info);
|
||||
void dcopy_(const integer *n, const doublereal *dx, const integer *incx, doublereal *dy, const integer *incy);
|
||||
|
||||
///* Subroutine */ int dgetrs_(char *trans, integer *n, integer *nrhs,
|
||||
// doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *
|
||||
// ldb, integer *info, int len);
|
||||
doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy);
|
||||
|
||||
///* Subroutine */ int dgetrf_(integer *m, integer *n, doublereal *a, integer *
|
||||
// lda, integer *ipiv, integer *info);
|
||||
void daxpy_(integer* n, doublereal* a, doublereal* x, integer* incx,
|
||||
doublereal* y, integer* incy);
|
||||
|
||||
///* Subroutine */ int dgetri_(integer *n, doublereal *a, integer *
|
||||
// lda, integer *ipiv, doublereal *work, integer *lwork, integer *info);
|
||||
void dscal_(integer *n, doublereal *da, doublereal *dx, integer *incx);
|
||||
|
||||
// /* Subroutine */ int dgemv_(char *trans, integer *m, integer *n, doublereal *
|
||||
// alpha, doublereal *a, integer *lda, doublereal *x, integer *incx,
|
||||
// doublereal *beta, doublereal *y, integer *incy, int len);
|
||||
|
||||
void dcopy_(const integer *n, const doublereal *dx, const integer *incx, doublereal *dy, const integer *incy);
|
||||
//doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy);
|
||||
|
||||
doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy);
|
||||
|
||||
void daxpy_(integer* n, doublereal* a, doublereal* x, integer* incx,
|
||||
doublereal* y, integer* incy);
|
||||
void dscal_(integer *n, doublereal *da, doublereal *dx, integer *incx);
|
||||
integer idamax_(integer* n, doublereal* a, integer* incx);
|
||||
|
||||
integer idamax_(integer* n, doublereal* a, integer* incx);
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,19 +1,93 @@
|
|||
|
||||
/* C interface for Fortran DPOLFT subroutine */
|
||||
/**
|
||||
* @file polyfit.h C interface for Fortran DPOLFT subroutine
|
||||
*/
|
||||
/*
|
||||
* $Author$
|
||||
* $Date$
|
||||
* $Revision$
|
||||
*
|
||||
* Copyright 2001-2003 California Institute of Technology
|
||||
* See file License.txt for licensing information
|
||||
*
|
||||
*/
|
||||
|
||||
#ifndef CT_POLYFIT_H
|
||||
#define CT_POLYFIT_H
|
||||
|
||||
#include <vector>
|
||||
//using namespace std;
|
||||
|
||||
#include "ct_defs.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w,
|
||||
int maxdeg, int& ndeg, doublereal eps, doublereal* r);
|
||||
|
||||
//! Fits a polynomial function to a set of data points
|
||||
/*!
|
||||
* Given a collection of points X(I) and a set of values Y(I) which
|
||||
* correspond to some function or measurement at each of the X(I),
|
||||
* subroutine DPOLFT computes the weighted least-squares polynomial
|
||||
* fits of all degrees up to some degree either specified by the user
|
||||
* or determined by the routine. The fits thus obtained are in
|
||||
* orthogonal polynomial form. Subroutine DP1VLU may then be
|
||||
* called to evaluate the fitted polynomials and any of their
|
||||
* derivatives at any point. The subroutine DPCOEF may be used to
|
||||
* express the polynomial fits as powers of (X-C) for any specified
|
||||
* point C.
|
||||
*
|
||||
* @param n The number of data points.
|
||||
*
|
||||
* @param x A set of grid points on which the data is specified.
|
||||
* The array of values of the independent variable. These
|
||||
* values may appear in any order and need not all be
|
||||
* distinct. There are n of them.
|
||||
*
|
||||
* @param y array of corresponding function values. There are n of them
|
||||
*
|
||||
* @param w array of positive values to be used as weights. If
|
||||
* W[0] is negative, DPOLFT will set all the weights
|
||||
* to 1.0, which means unweighted least squares error
|
||||
* will be minimized. To minimize relative error, the
|
||||
* user should set the weights to: W(I) = 1.0/Y(I)**2,
|
||||
* I = 1,...,N .
|
||||
*
|
||||
* @param maxdeg maximum degree to be allowed for polynomial fit.
|
||||
* MAXDEG may be any non-negative integer less than N.
|
||||
* Note -- MAXDEG cannot be equal to N-1 when a
|
||||
* statistical test is to be used for degree selection,
|
||||
* i.e., when input value of EPS is negative.
|
||||
*
|
||||
* @param ndeg output degree of the fit computed.
|
||||
*
|
||||
* @param eps Specifies the criterion to be used in determining
|
||||
* the degree of fit to be computed.
|
||||
* (1) If EPS is input negative, DPOLFT chooses the
|
||||
* degree based on a statistical F test of
|
||||
* significance. One of three possible
|
||||
* significance levels will be used: .01, .05 or
|
||||
* .10. If EPS=-1.0 , the routine will
|
||||
* automatically select one of these levels based
|
||||
* on the number of data points and the maximum
|
||||
* degree to be considered. If EPS is input as
|
||||
* -.01, -.05, or -.10, a significance level of
|
||||
* .01, .05, or .10, respectively, will be used.
|
||||
* (2) If EPS is set to 0., DPOLFT computes the
|
||||
* polynomials of degrees 0 through MAXDEG .
|
||||
* (3) If EPS is input positive, EPS is the RMS
|
||||
* error tolerance which must be satisfied by the
|
||||
* fitted polynomial. DPOLFT will increase the
|
||||
* degree of fit until this criterion is met or
|
||||
* until the maximum degree is reached.
|
||||
*
|
||||
* @param r Output vector containing the first ndeg+1 Taylor coefficients
|
||||
*
|
||||
* P(X) = r[0] + r[1]*(X-C) + ... + r[ndeg] * (X-C)**ndeg
|
||||
* ( here C = 0.0)
|
||||
*
|
||||
* @return Returned value is the value of the rms of the interpolated
|
||||
* function at x.
|
||||
*/
|
||||
doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w,
|
||||
int maxdeg, int& ndeg, doublereal eps, doublereal* r);
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue