From fae139765e4c93db92fa60ac1bcedee77a1c39c0 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 30 Jul 2008 16:48:25 +0000 Subject: [PATCH] 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. --- Cantera/src/numerics/funcs.cpp | 195 ++++++++++++++++++++++++--------- Cantera/src/numerics/funcs.h | 33 +++++- Cantera/src/numerics/lapack.h | 31 +----- Cantera/src/numerics/polyfit.h | 86 ++++++++++++++- 4 files changed, 258 insertions(+), 87 deletions(-) diff --git a/Cantera/src/numerics/funcs.cpp b/Cantera/src/numerics/funcs.cpp index 0aaa08160..81ce55810 100755 --- a/Cantera/src/numerics/funcs.cpp +++ b/Cantera/src/numerics/funcs.cpp @@ -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; + } } diff --git a/Cantera/src/numerics/funcs.h b/Cantera/src/numerics/funcs.h index d7b98e161..8774f195e 100644 --- a/Cantera/src/numerics/funcs.h +++ b/Cantera/src/numerics/funcs.h @@ -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 diff --git a/Cantera/src/numerics/lapack.h b/Cantera/src/numerics/lapack.h index f4ad15a4c..1ec74b307 100755 --- a/Cantera/src/numerics/lapack.h +++ b/Cantera/src/numerics/lapack.h @@ -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); } diff --git a/Cantera/src/numerics/polyfit.h b/Cantera/src/numerics/polyfit.h index 916f87929..b8436ae34 100755 --- a/Cantera/src/numerics/polyfit.h +++ b/Cantera/src/numerics/polyfit.h @@ -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 -//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