First commit of this directory
This commit is contained in:
Harry Moffat
parent
507055dd15
commit
6f4a6f6879
27 changed files with 15831 additions and 0 deletions
1
ext/f2c_math/.cvsignore
Normal file
1
ext/f2c_math/.cvsignore
Normal file
|
|
@ -0,0 +1 @@
|
|||
Makefile
|
||||
96
ext/f2c_math/Makefile.in
Normal file
96
ext/f2c_math/Makefile.in
Normal file
|
|
@ -0,0 +1,96 @@
|
|||
#/bin/sh
|
||||
#
|
||||
# $Source$
|
||||
# $Author$
|
||||
# $Revision$
|
||||
# $Date$
|
||||
#
|
||||
|
||||
.SUFFIXES :
|
||||
.SUFFIXES : .c .cpp .d .o
|
||||
|
||||
# the directory where the Cantera libraries are located
|
||||
CANTERA_LIBDIR=@buildlib@
|
||||
|
||||
# the directory where Cantera include files may be found.
|
||||
CANTERA_INCDIR=@ctroot@/build/include/cantera
|
||||
|
||||
# the C++ compiler
|
||||
CXX = @CXX@
|
||||
|
||||
# the C compiler
|
||||
CC = @CC@
|
||||
|
||||
# C++ compile flags
|
||||
CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT)
|
||||
|
||||
# Local include files
|
||||
CXX_INCLUDES=-I../f2c_libs
|
||||
|
||||
# How to compile the dependency file
|
||||
.c.d:
|
||||
g++ -MM $(CXX_FLAGS) $(CXX_INCLUDES) $*.c > $*.d
|
||||
|
||||
.cpp.d:
|
||||
g++ -MM $(CXX_FLAGS) $(CXX_INCLUDES) $*.cpp > $*.d
|
||||
|
||||
# How to compile a C file
|
||||
.c.o:
|
||||
@CC@ -c $< @DEFS@ $(CXX_FLAGS) $(CXX_INCLUDES)
|
||||
|
||||
# How to compile a Cpp file
|
||||
.cpp.o:
|
||||
@CXX@ -c $< @DEFS@ $(CXX_FLAGS) $(CXX_INCLUDES)
|
||||
|
||||
|
||||
# -----------------------------------------------
|
||||
|
||||
LIB = @buildlib@/libctmath.a
|
||||
|
||||
all: $(LIB)
|
||||
|
||||
OBJS = \
|
||||
mach.o \
|
||||
ddaspk.o \
|
||||
dgbefa.o \
|
||||
dgbsl.o \
|
||||
dgefa.o \
|
||||
dgesl.o \
|
||||
dp1vlu.o \
|
||||
dpcoef.o \
|
||||
dpolft.o \
|
||||
fdump.o \
|
||||
j4save.o \
|
||||
pcoef.o \
|
||||
polfit.o \
|
||||
pvalue.o \
|
||||
xercnt.o \
|
||||
xerhlt.o \
|
||||
xermsg.o \
|
||||
xerprn.o \
|
||||
xersve.o \
|
||||
xgetua.o
|
||||
|
||||
SRCS = $(OBJS:.o=.cpp)
|
||||
|
||||
# List of dependency files to be created
|
||||
DEPENDS=$(OBJS:.o=.d)
|
||||
|
||||
# How to make the static library
|
||||
$(LIB): $(OBJS)
|
||||
@ARCHIVE@ $(LIB) $(OBJS)
|
||||
|
||||
# ------------------------------------------------
|
||||
# Utility Targets
|
||||
|
||||
clean:
|
||||
$(RM) $(OBJS) $(LIB) *.d .depends
|
||||
|
||||
# depends target
|
||||
depends:
|
||||
$(RM) *.d .depends
|
||||
@MAKE@ .depends
|
||||
|
||||
.depends: $(DEPENDS)
|
||||
cat *.d > .depends
|
||||
|
||||
646
ext/f2c_math/cblas.h
Normal file
646
ext/f2c_math/cblas.h
Normal file
|
|
@ -0,0 +1,646 @@
|
|||
// -*- C++ -*-
|
||||
|
||||
// ============================================= //
|
||||
// die double-Versionen der BLAS Level 1 und 2 //
|
||||
// ============================================= //
|
||||
|
||||
#ifndef CBLAS1_H
|
||||
// ============================================================================
|
||||
|
||||
// generate a plane rotation
|
||||
void drotg( double *a, double *b, double *c, double *s );
|
||||
|
||||
|
||||
#if 0
|
||||
// generate a modified plane rotation
|
||||
void drotmg( double *d1, double *d2, double *a, double b, double *param );
|
||||
#endif
|
||||
|
||||
// apply a plane rotation
|
||||
void drot( int n, double *x, int incx, double *y, int incy, double c,
|
||||
double s );
|
||||
|
||||
|
||||
#if 0
|
||||
// apply a modified plane rotation
|
||||
void drotm( int n, double *x, int incx, double *y, int incy, double *param );
|
||||
#endif
|
||||
|
||||
|
||||
// x <=> y
|
||||
void dswap( int n, double *x, int incx, double *y, int incy );
|
||||
|
||||
|
||||
// x <= a*x
|
||||
void dscal( int n, double alpha, double *x, int incx );
|
||||
|
||||
|
||||
// y <= x
|
||||
void dcopy( int n, const double *x, int incx, double *y, int incy );
|
||||
|
||||
|
||||
// y <= a*x+y
|
||||
void daxpy( int n, double alpha, const double *x, int incx, double *y,
|
||||
int incy );
|
||||
|
||||
|
||||
// dot <= x^T*y
|
||||
double ddot( int n, const double *x, int incx, const double *y, int incy );
|
||||
|
||||
|
||||
// dnrm2 <= |x|_2
|
||||
double dnrm2( int n, const double *x, int incx );
|
||||
|
||||
|
||||
// asum <= |x|_1
|
||||
double dasum( int n, const double *x, int incx );
|
||||
|
||||
|
||||
// idamax <= first k such that |x_k| = max|x_i|
|
||||
int idamax( int n, const double *x, int incx );
|
||||
|
||||
// ============================================================================
|
||||
#endif // CBLAS1_H
|
||||
|
||||
|
||||
#ifndef CBLAS2_H
|
||||
// ============================================================================
|
||||
|
||||
|
||||
enum MatrixTranspose { NoTranspose=0, Transpose=1, ConjugateTranspose=2 };
|
||||
enum MatrixTriangle { UpperTriangle=0, LowerTriangle=1 };
|
||||
enum MatrixUnitTriangular { UnitTriangular=0, NotUnitTriangular=1 };
|
||||
|
||||
|
||||
// ============================================================================
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y, y <= alpha*A^T*x + beta*y, A-(m,n)
|
||||
void dgemv( MatrixTranspose trans, int m, int n, double alpha,
|
||||
const double *A, int ldA, const double *x, int incx,
|
||||
double beta, double *y, int incy );
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y, y <= alpha*A^T*x + beta*y, A-(m,n)
|
||||
void dgbmv( MatrixTranspose trans, int m, int n, int kl, int ku, double alpha,
|
||||
const double *A, int ldA, const double *x, int incx, double *beta,
|
||||
double *y, int incy );
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y
|
||||
void dsymv( MatrixTriangle uplo, int n, double alpha, const double *A, int ldA,
|
||||
const double *x, int incx, double beta, double *y, int incy );
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y
|
||||
void dsbmv( MatrixTriangle uplo, int n, int k, double alpha, double *A,
|
||||
int ldA, const double *x, int incx, double beta, double *y,
|
||||
int *incy );
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y
|
||||
void dspmv( MatrixTriangle uplo, int n, double alpha, const double *AP,
|
||||
const double *x, int incx, double beta, double *y, int incy );
|
||||
|
||||
|
||||
// x <= A*x, x <= A^T*x
|
||||
void dtrmv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, const double *A, int ldA,
|
||||
double *x, int incx );
|
||||
|
||||
|
||||
// x <= A*x, x <= A^T*x
|
||||
void dtbmv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, int k, const double *A, int ldA,
|
||||
double *x, int incx );
|
||||
|
||||
|
||||
// x <= A*x, x <= A^T*x
|
||||
void dtpmv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, int k, const double *AP,
|
||||
double *x, int incx );
|
||||
|
||||
|
||||
// x <= A^{-1}*x, x <= A^{-T}*x
|
||||
void dtrsv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, const double *A, int ldA,
|
||||
double *x, int incx );
|
||||
|
||||
|
||||
// x <= A^{-1}*x, x <= A^{-T}*x
|
||||
void dtbsv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, int k, const double *A, int ldA,
|
||||
double *x, int incx );
|
||||
|
||||
|
||||
// x <= A^{-1}*x, x <= A^{-T}*x
|
||||
void dtpsv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, int k, const double *AP,
|
||||
double *x, int incx );
|
||||
|
||||
|
||||
// A <= alpha*x*y^T + A, A-(m,n)
|
||||
void dger( int m, int n, double alpha, const double *x, int incx,
|
||||
const double *y, int incy, double *A, int ldA );
|
||||
|
||||
|
||||
// A <= alpha*x*x^T + A
|
||||
void dsyr( MatrixTriangle uplo, int n, double alpha, const double *x,
|
||||
int incx, double *A, int ldA );
|
||||
|
||||
|
||||
// A <= alpha*x*x^T + A
|
||||
void dspr( MatrixTriangle uplo, int n, double alpha, const double *x,
|
||||
int incx, double *AP );
|
||||
|
||||
|
||||
// A <= alpha*x*y^T + alpha*y*x^T + A
|
||||
void dsyr2( MatrixTriangle uplo, int n, double alpha, const double *x,
|
||||
int incx, const double *y, int incy, double *A, int ldA );
|
||||
|
||||
|
||||
// A <= alpha*x*y^T + alpha*y*x^T + A
|
||||
void dspr2( MatrixTriangle uplo, int n, double alpha, const double *x,
|
||||
int incx, const double *y, int incy, double *AP );
|
||||
|
||||
// ============================================================================
|
||||
#endif // CBLAS2_H
|
||||
|
||||
|
||||
#ifndef BLAS1_H
|
||||
#define BLAS1_H
|
||||
// ============================================================================
|
||||
|
||||
// generate a plane rotation
|
||||
extern "C"
|
||||
void drotg_( double *a, double *b, double *c, double *s );
|
||||
|
||||
|
||||
#if 0
|
||||
// generate a modified plane rotation
|
||||
extern "C"
|
||||
void drotmg_( double *d1, double *d2, double *a, double *b, double *param );
|
||||
#endif
|
||||
|
||||
|
||||
// apply a plane rotation
|
||||
extern "C"
|
||||
void drot_( int *n, double *x, int *incx, double *y, int *incy,
|
||||
double *c, double *s );
|
||||
|
||||
|
||||
#if 0
|
||||
// apply a modified plane rotation
|
||||
extern "C"
|
||||
void drotm_( int *n, double *x, int *incx, double *y, int *incy,
|
||||
double *param );
|
||||
#endif
|
||||
|
||||
|
||||
// x <=> y
|
||||
extern "C"
|
||||
void dswap_( const int *n, double *x, const int *incx, double *y,
|
||||
const int *incy );
|
||||
|
||||
// x <= a*x
|
||||
extern "C"
|
||||
void dscal_( const int *n, const double *alpha, double *x, const int *incx );
|
||||
|
||||
|
||||
// y <= x
|
||||
extern "C"
|
||||
void dcopy_( const int *n, const double *x, const int *incx, double *y,
|
||||
const int *incy );
|
||||
|
||||
|
||||
// y <= a*x+y
|
||||
extern "C"
|
||||
void daxpy_( const int *n, const double *alpha, const double *x,
|
||||
const int *incx, double *y, const int *incy );
|
||||
|
||||
|
||||
// dot <= x^T*y
|
||||
extern "C"
|
||||
double ddot_( const int *n, const double *x, const int *incx, const double *y,
|
||||
const int *incy );
|
||||
|
||||
|
||||
// dnrm2 <= |x|_2
|
||||
extern "C"
|
||||
double dnrm2_( const int *n, const double *x, const int *incx );
|
||||
|
||||
|
||||
// asum <= |x|_1
|
||||
extern "C"
|
||||
double dasum_( const int *n, const double *x, const int *incx );
|
||||
|
||||
|
||||
// idamax <= first k such that |x_k| = max|x_i|
|
||||
extern "C"
|
||||
int idamax_( const int *n, const double *x, const int *incx );
|
||||
|
||||
|
||||
// ============================================================================
|
||||
#endif // BLAS1_H
|
||||
|
||||
|
||||
|
||||
#ifndef CBLAS1_H
|
||||
#define CBLAS1_H
|
||||
// ============================================================================
|
||||
|
||||
|
||||
#ifdef __linux__ // muss dnorm2 f"ur linux neu implementieren
|
||||
# include <math.h>
|
||||
#endif
|
||||
|
||||
inline
|
||||
void drotg( double *a, double *b, double *c, double *s ) {
|
||||
drotg_(a,b,c,s);
|
||||
}
|
||||
|
||||
#if 0
|
||||
inline
|
||||
void drotmg( double *d1, double *d2, double *a, double b, double *param ) {
|
||||
drotmg_(d1,d2,a,&b,param);
|
||||
}
|
||||
#endif
|
||||
|
||||
inline
|
||||
void drot( int n, double *x, int incx, double *y, int incy, double c,
|
||||
double s ) {
|
||||
drot_(&n,x,&incx,y,&incy,&c,&s);
|
||||
}
|
||||
|
||||
#if 0
|
||||
inline
|
||||
void drotm( int n, double *x, int incx, double *y, int incy, double *param ) {
|
||||
drotm_(&n,x,&incx,y,&incy,param);
|
||||
}
|
||||
#endif
|
||||
|
||||
inline
|
||||
void dswap( int n, double *x, int incx, double *y, int incy ) {
|
||||
dswap_(&n,x,&incx,y,&incy);
|
||||
}
|
||||
|
||||
inline
|
||||
void dscal( int n, double alpha, double *x, int incx ) {
|
||||
int nn = n;
|
||||
int incxx = incx;
|
||||
double aa = alpha;
|
||||
dscal_(&nn,&aa,x,&incxx);
|
||||
}
|
||||
|
||||
inline
|
||||
void dcopy( int n, const double *x, int incx, double *y, int incy ) {
|
||||
int nn = n;
|
||||
int incxx = incx;
|
||||
int incyy = incy;
|
||||
dcopy_(&nn,x,&incxx,y,&incyy);
|
||||
}
|
||||
|
||||
inline
|
||||
void daxpy( int n, double alpha, const double *x, int incx, double *y,
|
||||
int incy ) {
|
||||
double aa = alpha;
|
||||
int incxx = incx;
|
||||
int incyy = incy;
|
||||
daxpy_(&n,&aa,x,&incxx,y,&incyy);
|
||||
}
|
||||
|
||||
inline
|
||||
double ddot( int n, const double *x, int incx, const double *y, int incy ) {
|
||||
int nn = n;
|
||||
int incxx = incx;
|
||||
int incyy = incy;
|
||||
return ddot_(&nn,x,&incxx,y,&incyy);
|
||||
}
|
||||
|
||||
inline
|
||||
double dnrm2( int n, const double *x, int incx ) {
|
||||
int nn = n;
|
||||
int incxx = incx;
|
||||
#ifdef __linux__ // fehlerhafte Berechnung
|
||||
double d=0.;
|
||||
while ( nn-- )
|
||||
d+=(*x)*(*x), x+=incxx;
|
||||
return sqrt(d);
|
||||
#else // unter nicht-Linux korrekt
|
||||
return dnrm2_(&nn,x,&incxx);
|
||||
#endif
|
||||
}
|
||||
|
||||
inline
|
||||
double dasum( int n, const double *x, int incx ) {
|
||||
return dasum_(&n,x,&incx);
|
||||
}
|
||||
|
||||
inline
|
||||
int idamax( int n, const double *x, int incx ) {
|
||||
return idamax_(&n,x,&incx);
|
||||
}
|
||||
|
||||
|
||||
// ============================================================================
|
||||
#endif // CBLAS1_H
|
||||
|
||||
|
||||
#ifndef BLAS2_H
|
||||
#define BLAS2_H
|
||||
// ============================================================================
|
||||
|
||||
// y <= alpha*A*x + beta*y, y <= alpha*A^T*x + beta*y, A-(m,n)
|
||||
//extern "C"
|
||||
//void dgemv_( const char *trans, const int *m, const int *n,
|
||||
// const double *alpha, const double *A, const int *ldA,
|
||||
// const double *x, const int *incx,
|
||||
// const double *beta, double *y, const int *incy );
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y, y <= alpha*A^T*x + beta*y, A-(m,n)
|
||||
extern "C"
|
||||
void dgbmv_( const char *trans, const int *m, const int *n, const int *kl,
|
||||
const int *ku, const double *alpha, const double *A,
|
||||
const int *ldA, const double *x, const int *incx,
|
||||
const double *beta, double *y, const int *incy );
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y
|
||||
extern "C"
|
||||
void dsymv_( const char *uplo, const int *n, const double *alpha,
|
||||
const double *A, const int *ldA, const double *x, const int *incx,
|
||||
const double *beta, double *y, const int *incy );
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y
|
||||
extern "C"
|
||||
void dsbmv_( const char *uplo, const int *n, const int *k, const double *alpha,
|
||||
const double *A, const int *ldA, const double *x, const int *incx,
|
||||
const double *beta, double *y, const int *incy );
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y
|
||||
extern "C"
|
||||
void dspmv_( const char *uplo, const int *n, const double *alpha,
|
||||
const double *AP, const double *x, const int *incx,
|
||||
const double *beta, double *y, const int *incy );
|
||||
|
||||
|
||||
// x <= A*x, x <= A^T*x
|
||||
extern "C"
|
||||
void dtrmv_( const char *uplo, const char *trans, const char *diag,
|
||||
const int *n, const double *A, const int *ldA,
|
||||
double *x, const int *incx );
|
||||
|
||||
|
||||
// x <= A*x, x <= A^T*x
|
||||
extern "C"
|
||||
void dtbmv_( const char *uplo, const char *trans, const char *diag,
|
||||
const int *n, const int *k, const double *A, const int *ldA,
|
||||
double *x, const int *incx );
|
||||
|
||||
|
||||
// x <= A*x, x <= A^T*x
|
||||
extern "C"
|
||||
void dtpmv_( const char *uplo, const char *trans, const char *diag,
|
||||
const int *n, const double *AP, double *x, const int *incx );
|
||||
|
||||
|
||||
// x <= A^{-1}*x, x <= A^{-T}*x
|
||||
extern "C"
|
||||
void dtrsv_( const char *uplo, const char *trans, const char *diag,
|
||||
const int *n, const double *A, const int *ldA,
|
||||
double *x, const int *incx );
|
||||
|
||||
|
||||
// x <= A^{-1}*x, x <= A^{-T}*x
|
||||
extern "C"
|
||||
void dtbsv_( const char *uplo, const char *trans, const char *diag,
|
||||
const int *n, const int *k, const double *A, const int *ldA,
|
||||
double *x, const int *incx );
|
||||
|
||||
|
||||
// x <= A^{-1}*x, x <= A^{-T}*x
|
||||
extern "C"
|
||||
void dtpsv_( const char *uplo, const char *trans, const char *diag,
|
||||
const int *n, const double *AP, double *x, const int *incx );
|
||||
|
||||
|
||||
// A <= alpha*x*y^T + A, A-(m,n)
|
||||
extern "C"
|
||||
void dger_( const int *m, const int *n, const double *alpha, const double *x,
|
||||
const int *incx, const double *y, const int *incy, double *A,
|
||||
const int *ldA );
|
||||
|
||||
|
||||
// A <= alpha*x*x^T + A
|
||||
extern "C"
|
||||
void dsyr_( const char *uplo, const int *n, const double *alpha,
|
||||
const double *x, const int *incx, double *A, const int *ldA );
|
||||
|
||||
|
||||
// A <= alpha*x*x^T + A
|
||||
extern "C"
|
||||
void dspr_( const char *uplo, const int *n, const double *alpha,
|
||||
const double *x, const int *incx, double *AP );
|
||||
|
||||
|
||||
// A <= alpha*x*y^T + alpha*y*x^T + A
|
||||
extern "C"
|
||||
void dsyr2_( const char *uplo, const int *n, const double *alpha,
|
||||
const double *x, const int *incx, const double *y,
|
||||
const int *incy, double *A, const int *ldA );
|
||||
|
||||
|
||||
// A <= alpha*x*y^T + alpha*y*x^T + A
|
||||
extern "C"
|
||||
void dspr2_( const char *uplo, const int *n, const double *alpha,
|
||||
const double *x, const int *incx, const double *y,
|
||||
const int *incy, double *AP );
|
||||
|
||||
// ============================================================================
|
||||
#endif // BLAS2_H
|
||||
|
||||
|
||||
#ifndef CBLAS2_H
|
||||
#define CBLAS2_H
|
||||
// ============================================================================
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y, y <= alpha*A^T*x + beta*y, A-(m,n)
|
||||
inline
|
||||
void dgemv( MatrixTranspose trans, int m, int n, double alpha,
|
||||
const double *A, int ldA, const double *x, int incx,
|
||||
double beta, double *y, int incy ) {
|
||||
const char *T[3] = { "N", "T", 0 };
|
||||
int mm = m;
|
||||
int nn = n;
|
||||
double aa = alpha;
|
||||
double bb = beta;
|
||||
int ldaa = ldA;
|
||||
int incxx = incx;
|
||||
int incyy = incy;
|
||||
dgemv_(T[(int)trans],&mm,&nn,&aa,A,&ldaa,x,&incxx,&bb,y,&incyy,1);
|
||||
}
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y, y <= alpha*A^T*x + beta*y, A-(m,n)
|
||||
inline
|
||||
void dgbmv( MatrixTranspose trans, int m, int n, int kl, int ku, double alpha,
|
||||
const double *A, int ldA, const double *x, int incx, double beta,
|
||||
double *y, int incy ) {
|
||||
const char *T[3] = { "N", "T" };
|
||||
dgbmv_(T[(int)trans],&m,&n,&kl,&ku,&alpha,A,&ldA,x,&incx,&beta,y,&incy);
|
||||
}
|
||||
|
||||
// y <= alpha*A*x + beta*y
|
||||
inline
|
||||
void dsymv( MatrixTriangle uplo, int n, double alpha, const double *A, int ldA,
|
||||
const double *x, int incx, double beta, double *y, int incy ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
dsymv_(UL[(int)uplo],&n,&alpha,A,&ldA,x,&incx,&beta,y,&incy);
|
||||
}
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y
|
||||
inline
|
||||
void dsbmv( MatrixTriangle uplo, int n, int k, double alpha, double *A,
|
||||
int ldA, const double *x, int incx, double beta, double *y,
|
||||
int incy ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
dsbmv_(UL[(int)uplo],&n,&k,&alpha,A,&ldA,x,&incx,&beta,y,&incy);
|
||||
}
|
||||
|
||||
|
||||
// y <= alpha*A*x + beta*y
|
||||
inline
|
||||
void dspmv( MatrixTriangle uplo, int n, double alpha, const double *AP,
|
||||
const double *x, int incx, double beta, double *y, int incy ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
dspmv_(UL[(int)uplo],&n,&alpha,AP,x,&incx,&beta,y,&incy);
|
||||
}
|
||||
|
||||
|
||||
// x <= A*x, x <= A^T*x
|
||||
inline
|
||||
void dtrmv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, const double *A, int ldA,
|
||||
double *x, int incx ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
const char *T[3] = { "N", "T", 0 };
|
||||
const char *D[2] = { "U", "N" };
|
||||
dtrmv_(UL[(int)uplo],T[(int)trans],D[(int)diag],&n,A,&ldA,x,&incx);
|
||||
}
|
||||
|
||||
|
||||
// x <= A*x, x <= A^T*x
|
||||
inline
|
||||
void dtbmv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, int k, const double *A, int ldA,
|
||||
double *x, int incx ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
const char *T[3] = { "N", "T", 0 };
|
||||
const char *D[2] = { "U", "N" };
|
||||
dtbmv_(UL[(int)uplo],T[(int)trans],D[(int)diag],&n,&k,A,&ldA,x,&incx);
|
||||
}
|
||||
|
||||
|
||||
// x <= A*x, x <= A^T*x
|
||||
inline
|
||||
void dtpmv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, const double *AP,
|
||||
double *x, int incx ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
const char *T[3] = { "N", "T", 0 };
|
||||
const char *D[2] = { "U", "N" };
|
||||
dtpmv_(UL[(int)uplo],T[(int)trans],D[(int)diag],&n,AP,x,&incx);
|
||||
}
|
||||
|
||||
|
||||
// x <= A^{-1}*x, x <= A^{-T}*x
|
||||
inline
|
||||
void dtrsv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, const double *A, int ldA,
|
||||
double *x, int incx ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
const char *T[3] = { "N", "T", 0 };
|
||||
const char *D[2] = { "U", "N" };
|
||||
dtrsv_(UL[(int)uplo],T[(int)trans],D[(int)diag],&n,A,&ldA,x,&incx);
|
||||
}
|
||||
|
||||
|
||||
// x <= A^{-1}*x, x <= A^{-T}*x
|
||||
inline
|
||||
void dtbsv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, int k, const double *A, int ldA,
|
||||
double *x, int incx ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
const char *T[3] = { "N", "T", 0 };
|
||||
const char *D[2] = { "U", "N" };
|
||||
dtbsv_(UL[(int)uplo],T[(int)trans],D[(int)diag],&n,&k,A,&ldA,x,&incx);
|
||||
}
|
||||
|
||||
|
||||
// x <= A^{-1}*x, x <= A^{-T}*x
|
||||
inline
|
||||
void dtpsv( MatrixTriangle uplo, MatrixTranspose trans,
|
||||
MatrixUnitTriangular diag, int n, const double *AP,
|
||||
double *x, int incx ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
const char *T[3] = { "N", "T", 0 };
|
||||
const char *D[2] = { "U", "N" };
|
||||
int nn = n;
|
||||
int incxx = incx;
|
||||
dtpsv_(UL[(int)uplo],T[(int)trans],D[(int)diag],&nn,AP,x,&incxx);
|
||||
}
|
||||
|
||||
|
||||
// A <= alpha*x*y^T + A, A-(m,n)
|
||||
inline
|
||||
void dger( int m, int n, double alpha, const double *x, int incx,
|
||||
const double *y, int incy, double *A, int ldA ) {
|
||||
dger_(&m,&n,&alpha,x,&incx,y,&incy,A,&ldA);
|
||||
}
|
||||
|
||||
|
||||
// A <= alpha*x*x^T + A
|
||||
inline
|
||||
void dsyr( MatrixTriangle uplo, int n, double alpha, const double *x,
|
||||
int incx, double *A, int ldA ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
dsyr_(UL[(int)uplo],&n,&alpha,x,&incx,A,&ldA);
|
||||
}
|
||||
|
||||
|
||||
// A <= alpha*x*x^T + A
|
||||
inline
|
||||
void dspr( MatrixTriangle uplo, int n, double alpha, const double *x,
|
||||
int incx, double *AP ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
dspr_(UL[(int)uplo],&n,&alpha,x,&incx,AP);
|
||||
}
|
||||
|
||||
|
||||
// A <= alpha*x*y^T + alpha*y*x^T + A
|
||||
inline
|
||||
void dsyr2( MatrixTriangle uplo, int n, double alpha, const double *x,
|
||||
int incx, const double *y, int incy, double *A, int ldA ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
dsyr2_(UL[(int)uplo],&n,&alpha,x,&incx,y,&incy,A,&ldA);
|
||||
}
|
||||
|
||||
|
||||
// A <= alpha*x*y^T + alpha*y*x^T + A
|
||||
inline
|
||||
void dspr2( MatrixTriangle uplo, int n, double alpha, const double *x,
|
||||
int incx, const double *y, int incy, double *AP ) {
|
||||
const char *UL[2] = { "U", "L" };
|
||||
dspr2_(UL[(int)uplo],&n,&alpha,x,&incx,y,&incy,AP);
|
||||
}
|
||||
|
||||
|
||||
// ============================================================================
|
||||
|
||||
|
||||
#endif // CBLAS2_H
|
||||
355
ext/f2c_math/daux.c
Normal file
355
ext/f2c_math/daux.c
Normal file
|
|
@ -0,0 +1,355 @@
|
|||
/* daux.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__1 = 1;
|
||||
static integer c__0 = 0;
|
||||
static logical c_false = FALSE_;
|
||||
static integer c__2 = 2;
|
||||
static logical c_true = TRUE_;
|
||||
|
||||
/* DOUBLE PRECISION FUNCTION D1MACH (IDUM) */
|
||||
/* INTEGER IDUM */
|
||||
/* C----------------------------------------------------------------------- */
|
||||
/* C THIS ROUTINE COMPUTES THE UNIT ROUNDOFF OF THE MACHINE IN DOUBLE */
|
||||
/* C PRECISION. THIS IS DEFINED AS THE SMALLEST POSITIVE MACHINE NUMBER */
|
||||
/* C U SUCH THAT 1.0D0 + U .NE. 1.0D0 (IN DOUBLE PRECISION). */
|
||||
/* C----------------------------------------------------------------------- */
|
||||
/* DOUBLE PRECISION U, COMP */
|
||||
/* U = 1.0D0 */
|
||||
/* 10 U = U*0.5D0 */
|
||||
/* COMP = 1.0D0 + U */
|
||||
/* IF (COMP .NE. 1.0D0) GO TO 10 */
|
||||
/* D1MACH = U*2.0D0 */
|
||||
/* RETURN */
|
||||
/* C----------------------- END OF FUNCTION D1MACH ------------------------ */
|
||||
/* END */
|
||||
/* DECK XERRWD */
|
||||
/* Subroutine */ int xerrwd_(char *msg, integer *nmes, integer *nerr, integer
|
||||
*level, integer *ni, integer *i1, integer *i2, integer *nr,
|
||||
doublereal *r1, doublereal *r2, ftnlen msg_len)
|
||||
{
|
||||
/* Format strings */
|
||||
static char fmt_10[] = "(1x,a)";
|
||||
static char fmt_20[] = "(6x,\002In above message, I1 =\002,i10)";
|
||||
static char fmt_30[] = "(6x,\002In above message, I1 =\002,i10,3x,\002I"
|
||||
"2 =\002,i10)";
|
||||
static char fmt_40[] = "(6x,\002In above message, R1 =\002,d21.13)";
|
||||
static char fmt_50[] = "(6x,\002In above, R1 =\002,d21.13,3x,\002R2 "
|
||||
"=\002,d21.13)";
|
||||
|
||||
/* Builtin functions */
|
||||
integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
|
||||
/* Subroutine */ int s_stop(char *, ftnlen);
|
||||
|
||||
/* Local variables */
|
||||
extern integer ixsav_(integer *, integer *, logical *);
|
||||
static integer lunit, mesflg;
|
||||
|
||||
/* Fortran I/O blocks */
|
||||
static cilist io___3 = { 0, 0, 0, fmt_10, 0 };
|
||||
static cilist io___4 = { 0, 0, 0, fmt_20, 0 };
|
||||
static cilist io___5 = { 0, 0, 0, fmt_30, 0 };
|
||||
static cilist io___6 = { 0, 0, 0, fmt_40, 0 };
|
||||
static cilist io___7 = { 0, 0, 0, fmt_50, 0 };
|
||||
|
||||
|
||||
/* ***BEGIN PROLOGUE XERRWD */
|
||||
/* ***SUBSIDIARY */
|
||||
/* ***PURPOSE Write error message with values. */
|
||||
/* ***LIBRARY MATHLIB */
|
||||
/* ***CATEGORY R3C */
|
||||
/* ***TYPE DOUBLE PRECISION (XERRWV-S, XERRWD-D) */
|
||||
/* ***AUTHOR Hindmarsh, Alan C., (LLNL) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* Subroutines XERRWD, XSETF, XSETUN, and the function routine IXSAV, */
|
||||
/* as given here, constitute a simplified version of the SLATEC error */
|
||||
/* handling package. */
|
||||
|
||||
/* All arguments are input arguments. */
|
||||
|
||||
/* MSG = The message (character array). */
|
||||
/* NMES = The length of MSG (number of characters). */
|
||||
/* NERR = The error number (not used). */
|
||||
/* LEVEL = The error level.. */
|
||||
/* 0 or 1 means recoverable (control returns to caller). */
|
||||
/* 2 means fatal (run is aborted--see note below). */
|
||||
/* NI = Number of integers (0, 1, or 2) to be printed with message. */
|
||||
/* I1,I2 = Integers to be printed, depending on NI. */
|
||||
/* NR = Number of reals (0, 1, or 2) to be printed with message. */
|
||||
/* R1,R2 = Reals to be printed, depending on NR. */
|
||||
|
||||
/* Note.. this routine is machine-dependent and specialized for use */
|
||||
/* in limited context, in the following ways.. */
|
||||
/* 1. The argument MSG is assumed to be of type CHARACTER, and */
|
||||
/* the message is printed with a format of (1X,A). */
|
||||
/* 2. The message is assumed to take only one line. */
|
||||
/* Multi-line messages are generated by repeated calls. */
|
||||
/* 3. If LEVEL = 2, control passes to the statement STOP */
|
||||
/* to abort the run. This statement may be machine-dependent. */
|
||||
/* 4. R1 and R2 are assumed to be in double precision and are printed */
|
||||
/* in D21.13 format. */
|
||||
|
||||
/* ***ROUTINES CALLED IXSAV */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 920831 DATE WRITTEN */
|
||||
/* 921118 Replaced MFLGSV/LUNSAV by IXSAV. (ACH) */
|
||||
/* 930329 Modified prologue to SLATEC format. (FNF) */
|
||||
/* 930407 Changed MSG from CHARACTER*1 array to variable. (FNF) */
|
||||
/* 930922 Minor cosmetic change. (FNF) */
|
||||
/* ***END PROLOGUE XERRWD */
|
||||
|
||||
/* *Internal Notes: */
|
||||
|
||||
/* For a different default logical unit number, IXSAV (or a subsidiary */
|
||||
/* routine that it calls) will need to be modified. */
|
||||
/* For a different run-abort command, change the statement following */
|
||||
/* statement 100 at the end. */
|
||||
/* ----------------------------------------------------------------------- */
|
||||
/* Subroutines called by XERRWD.. None */
|
||||
/* Function routine called by XERRWD.. IXSAV */
|
||||
/* ----------------------------------------------------------------------- */
|
||||
/* **End */
|
||||
|
||||
/* Declare arguments. */
|
||||
|
||||
|
||||
/* Declare local variables. */
|
||||
|
||||
|
||||
/* Get logical unit number and message print flag. */
|
||||
|
||||
/* ***FIRST EXECUTABLE STATEMENT XERRWD */
|
||||
lunit = ixsav_(&c__1, &c__0, &c_false);
|
||||
mesflg = ixsav_(&c__2, &c__0, &c_false);
|
||||
if (mesflg == 0) {
|
||||
goto L100;
|
||||
}
|
||||
|
||||
/* Write the message. */
|
||||
|
||||
io___3.ciunit = lunit;
|
||||
s_wsfe(&io___3);
|
||||
do_fio(&c__1, msg, msg_len);
|
||||
e_wsfe();
|
||||
if (*ni == 1) {
|
||||
io___4.ciunit = lunit;
|
||||
s_wsfe(&io___4);
|
||||
do_fio(&c__1, (char *)&(*i1), (ftnlen)sizeof(integer));
|
||||
e_wsfe();
|
||||
}
|
||||
if (*ni == 2) {
|
||||
io___5.ciunit = lunit;
|
||||
s_wsfe(&io___5);
|
||||
do_fio(&c__1, (char *)&(*i1), (ftnlen)sizeof(integer));
|
||||
do_fio(&c__1, (char *)&(*i2), (ftnlen)sizeof(integer));
|
||||
e_wsfe();
|
||||
}
|
||||
if (*nr == 1) {
|
||||
io___6.ciunit = lunit;
|
||||
s_wsfe(&io___6);
|
||||
do_fio(&c__1, (char *)&(*r1), (ftnlen)sizeof(doublereal));
|
||||
e_wsfe();
|
||||
}
|
||||
if (*nr == 2) {
|
||||
io___7.ciunit = lunit;
|
||||
s_wsfe(&io___7);
|
||||
do_fio(&c__1, (char *)&(*r1), (ftnlen)sizeof(doublereal));
|
||||
do_fio(&c__1, (char *)&(*r2), (ftnlen)sizeof(doublereal));
|
||||
e_wsfe();
|
||||
}
|
||||
|
||||
/* Abort the run if LEVEL = 2. */
|
||||
|
||||
L100:
|
||||
if (*level != 2) {
|
||||
return 0;
|
||||
}
|
||||
s_stop("", (ftnlen)0);
|
||||
/* ----------------------- End of Subroutine XERRWD ---------------------- */
|
||||
return 0;
|
||||
} /* xerrwd_ */
|
||||
|
||||
/* DECK XSETF */
|
||||
/* Subroutine */ int xsetf_(integer *mflag)
|
||||
{
|
||||
static integer junk;
|
||||
extern integer ixsav_(integer *, integer *, logical *);
|
||||
|
||||
/* ***BEGIN PROLOGUE XSETF */
|
||||
/* ***PURPOSE Reset the error print control flag. */
|
||||
/* ***LIBRARY MATHLIB */
|
||||
/* ***CATEGORY R3A */
|
||||
/* ***TYPE ALL (XSETF-A) */
|
||||
/* ***KEYWORDS ERROR CONTROL */
|
||||
/* ***AUTHOR Hindmarsh, Alan C., (LLNL) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* XSETF sets the error print control flag to MFLAG: */
|
||||
/* MFLAG=1 means print all messages (the default). */
|
||||
/* MFLAG=0 means no printing. */
|
||||
|
||||
/* ***SEE ALSO XERMSG, XERRWD, XERRWV */
|
||||
/* ***REFERENCES (NONE) */
|
||||
/* ***ROUTINES CALLED IXSAV */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 921118 DATE WRITTEN */
|
||||
/* 930329 Added SLATEC format prologue. (FNF) */
|
||||
/* 930407 Corrected SEE ALSO section. (FNF) */
|
||||
/* 930922 Made user-callable, and other cosmetic changes. (FNF) */
|
||||
/* ***END PROLOGUE XSETF */
|
||||
|
||||
/* Subroutines called by XSETF.. None */
|
||||
/* Function routine called by XSETF.. IXSAV */
|
||||
/* ----------------------------------------------------------------------- */
|
||||
/* **End */
|
||||
|
||||
/* ***FIRST EXECUTABLE STATEMENT XSETF */
|
||||
if (*mflag == 0 || *mflag == 1) {
|
||||
junk = ixsav_(&c__2, mflag, &c_true);
|
||||
}
|
||||
return 0;
|
||||
/* ----------------------- End of Subroutine XSETF ----------------------- */
|
||||
} /* xsetf_ */
|
||||
|
||||
/* DECK XSETUN */
|
||||
/* Subroutine */ int xsetun_(integer *lun)
|
||||
{
|
||||
static integer junk;
|
||||
extern integer ixsav_(integer *, integer *, logical *);
|
||||
|
||||
/* ***BEGIN PROLOGUE XSETUN */
|
||||
/* ***PURPOSE Reset the logical unit number for error messages. */
|
||||
/* ***LIBRARY MATHLIB */
|
||||
/* ***CATEGORY R3B */
|
||||
/* ***TYPE ALL (XSETUN-A) */
|
||||
/* ***KEYWORDS ERROR CONTROL */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* XSETUN sets the logical unit number for error messages to LUN. */
|
||||
|
||||
/* ***AUTHOR Hindmarsh, Alan C., (LLNL) */
|
||||
/* ***SEE ALSO XERMSG, XERRWD, XERRWV */
|
||||
/* ***REFERENCES (NONE) */
|
||||
/* ***ROUTINES CALLED IXSAV */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 921118 DATE WRITTEN */
|
||||
/* 930329 Added SLATEC format prologue. (FNF) */
|
||||
/* 930407 Corrected SEE ALSO section. (FNF) */
|
||||
/* 930922 Made user-callable, and other cosmetic changes. (FNF) */
|
||||
/* ***END PROLOGUE XSETUN */
|
||||
|
||||
/* Subroutines called by XSETUN.. None */
|
||||
/* Function routine called by XSETUN.. IXSAV */
|
||||
/* ----------------------------------------------------------------------- */
|
||||
/* **End */
|
||||
|
||||
/* ***FIRST EXECUTABLE STATEMENT XSETUN */
|
||||
if (*lun > 0) {
|
||||
junk = ixsav_(&c__1, lun, &c_true);
|
||||
}
|
||||
return 0;
|
||||
/* ----------------------- End of Subroutine XSETUN ---------------------- */
|
||||
} /* xsetun_ */
|
||||
|
||||
/* DECK IXSAV */
|
||||
integer ixsav_(integer *ipar, integer *ivalue, logical *iset)
|
||||
{
|
||||
/* Initialized data */
|
||||
|
||||
static integer lunit = -1;
|
||||
static integer lundef = 6;
|
||||
static integer mesflg = 1;
|
||||
|
||||
/* System generated locals */
|
||||
integer ret_val;
|
||||
|
||||
/* ***BEGIN PROLOGUE IXSAV */
|
||||
/* ***SUBSIDIARY */
|
||||
/* ***PURPOSE Save and recall error message control parameters. */
|
||||
/* ***LIBRARY MATHLIB */
|
||||
/* ***CATEGORY R3C */
|
||||
/* ***TYPE ALL (IXSAV-A) */
|
||||
/* ***AUTHOR Hindmarsh, Alan C., (LLNL) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* IXSAV saves and recalls one of two error message parameters: */
|
||||
/* LUNIT, the logical unit number to which messages are printed, and */
|
||||
/* MESFLG, the message print flag. */
|
||||
/* This is a modification of the SLATEC library routine J4SAVE. */
|
||||
|
||||
/* Saved local variables.. */
|
||||
/* LUNIT = Logical unit number for messages. */
|
||||
/* LUNDEF = Default logical unit number, data-loaded to 6 below */
|
||||
/* (may be machine-dependent). */
|
||||
/* MESFLG = Print control flag.. */
|
||||
/* 1 means print all messages (the default). */
|
||||
/* 0 means no printing. */
|
||||
|
||||
/* On input.. */
|
||||
/* IPAR = Parameter indicator (1 for LUNIT, 2 for MESFLG). */
|
||||
/* IVALUE = The value to be set for the parameter, if ISET = .TRUE. */
|
||||
/* ISET = Logical flag to indicate whether to read or write. */
|
||||
/* If ISET = .TRUE., the parameter will be given */
|
||||
/* the value IVALUE. If ISET = .FALSE., the parameter */
|
||||
/* will be unchanged, and IVALUE is a dummy argument. */
|
||||
|
||||
/* On return.. */
|
||||
/* IXSAV = The (old) value of the parameter. */
|
||||
|
||||
/* ***SEE ALSO XERMSG, XERRWD, XERRWV */
|
||||
/* ***ROUTINES CALLED NONE */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 921118 DATE WRITTEN */
|
||||
/* 930329 Modified prologue to SLATEC format. (FNF) */
|
||||
/* 941025 Minor modification re default unit number. (ACH) */
|
||||
/* ***END PROLOGUE IXSAV */
|
||||
|
||||
/* **End */
|
||||
/* ----------------------------------------------------------------------- */
|
||||
/* ----------------------------------------------------------------------- */
|
||||
/* The following Fortran-77 declaration is to cause the values of the */
|
||||
/* listed (local) variables to be saved between calls to this routine. */
|
||||
/* ----------------------------------------------------------------------- */
|
||||
|
||||
/* ***FIRST EXECUTABLE STATEMENT IXSAV */
|
||||
if (*ipar == 1) {
|
||||
if (lunit == -1) {
|
||||
lunit = lundef;
|
||||
}
|
||||
ret_val = lunit;
|
||||
if (*iset) {
|
||||
lunit = *ivalue;
|
||||
}
|
||||
}
|
||||
|
||||
if (*ipar == 2) {
|
||||
ret_val = mesflg;
|
||||
if (*iset) {
|
||||
mesflg = *ivalue;
|
||||
}
|
||||
}
|
||||
|
||||
return ret_val;
|
||||
/* ----------------------- End of Function IXSAV ------------------------- */
|
||||
} /* ixsav_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
8783
ext/f2c_math/ddaspk.c
Normal file
8783
ext/f2c_math/ddaspk.c
Normal file
File diff suppressed because it is too large
Load diff
253
ext/f2c_math/dgbefa.c
Normal file
253
ext/f2c_math/dgbefa.c
Normal file
|
|
@ -0,0 +1,253 @@
|
|||
/* dgbefa.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__1 = 1;
|
||||
|
||||
/* Subroutine */ int dgbfa_(doublereal *abd, integer *lda, integer *n,
|
||||
integer *ml, integer *mu, integer *ipvt, integer *info)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer abd_dim1, abd_offset, i__1, i__2, i__3, i__4;
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, j, k, l, m;
|
||||
static doublereal t;
|
||||
static integer i0, j0, j1, lm, mm, ju, jz, kp1, nm1;
|
||||
extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
|
||||
integer *), daxpy_(integer *, doublereal *, doublereal *, integer
|
||||
*, doublereal *, integer *);
|
||||
extern integer idamax_(integer *, doublereal *, integer *);
|
||||
|
||||
|
||||
/* dgbfa factors a double precision band matrix by elimination. */
|
||||
|
||||
/* dgbfa is usually called by dgbco, but it can be called */
|
||||
/* directly with a saving in time if rcond is not needed. */
|
||||
|
||||
/* on entry */
|
||||
|
||||
/* abd double precision(lda, n) */
|
||||
/* contains the matrix in band storage. the columns */
|
||||
/* of the matrix are stored in the columns of abd and */
|
||||
/* the diagonals of the matrix are stored in rows */
|
||||
/* ml+1 through 2*ml+mu+1 of abd . */
|
||||
/* see the comments below for details. */
|
||||
|
||||
/* lda integer */
|
||||
/* the leading dimension of the array abd . */
|
||||
/* lda must be .ge. 2*ml + mu + 1 . */
|
||||
|
||||
/* n integer */
|
||||
/* the order of the original matrix. */
|
||||
|
||||
/* ml integer */
|
||||
/* number of diagonals below the main diagonal. */
|
||||
/* 0 .le. ml .lt. n . */
|
||||
|
||||
/* mu integer */
|
||||
/* number of diagonals above the main diagonal. */
|
||||
/* 0 .le. mu .lt. n . */
|
||||
/* more efficient if ml .le. mu . */
|
||||
/* on return */
|
||||
|
||||
/* abd an upper triangular matrix in band storage and */
|
||||
/* the multipliers which were used to obtain it. */
|
||||
/* the factorization can be written a = l*u where */
|
||||
/* l is a product of permutation and unit lower */
|
||||
/* triangular matrices and u is upper triangular. */
|
||||
|
||||
/* ipvt integer(n) */
|
||||
/* an integer vector of pivot indices. */
|
||||
|
||||
/* info integer */
|
||||
/* = 0 normal value. */
|
||||
/* = k if u(k,k) .eq. 0.0 . this is not an error */
|
||||
/* condition for this subroutine, but it does */
|
||||
/* indicate that dgbsl will divide by zero if */
|
||||
/* called. use rcond in dgbco for a reliable */
|
||||
/* indication of singularity. */
|
||||
|
||||
/* band storage */
|
||||
|
||||
/* if a is a band matrix, the following program segment */
|
||||
/* will set up the input. */
|
||||
|
||||
/* ml = (band width below the diagonal) */
|
||||
/* mu = (band width above the diagonal) */
|
||||
/* m = ml + mu + 1 */
|
||||
/* do 20 j = 1, n */
|
||||
/* i1 = max0(1, j-mu) */
|
||||
/* i2 = min0(n, j+ml) */
|
||||
/* do 10 i = i1, i2 */
|
||||
/* k = i - j + m */
|
||||
/* abd(k,j) = a(i,j) */
|
||||
/* 10 continue */
|
||||
/* 20 continue */
|
||||
|
||||
/* this uses rows ml+1 through 2*ml+mu+1 of abd . */
|
||||
/* in addition, the first ml rows in abd are used for */
|
||||
/* elements generated during the triangularization. */
|
||||
/* the total number of rows needed in abd is 2*ml+mu+1 . */
|
||||
/* the ml+mu by ml+mu upper left triangle and the */
|
||||
/* ml by ml lower right triangle are not referenced. */
|
||||
|
||||
/* linpack. this version dated 08/14/78 . */
|
||||
/* cleve moler, university of new mexico, argonne national lab. */
|
||||
|
||||
/* subroutines and functions */
|
||||
|
||||
/* blas daxpy,dscal,idamax */
|
||||
/* fortran max0,min0 */
|
||||
|
||||
/* internal variables */
|
||||
|
||||
|
||||
|
||||
/* Parameter adjustments */
|
||||
abd_dim1 = *lda;
|
||||
abd_offset = 1 + abd_dim1;
|
||||
abd -= abd_offset;
|
||||
--ipvt;
|
||||
|
||||
/* Function Body */
|
||||
m = *ml + *mu + 1;
|
||||
*info = 0;
|
||||
|
||||
/* zero initial fill-in columns */
|
||||
|
||||
j0 = *mu + 2;
|
||||
j1 = min(*n,m) - 1;
|
||||
if (j1 < j0) {
|
||||
goto L30;
|
||||
}
|
||||
i__1 = j1;
|
||||
for (jz = j0; jz <= i__1; ++jz) {
|
||||
i0 = m + 1 - jz;
|
||||
i__2 = *ml;
|
||||
for (i__ = i0; i__ <= i__2; ++i__) {
|
||||
abd[i__ + jz * abd_dim1] = 0.;
|
||||
/* L10: */
|
||||
}
|
||||
/* L20: */
|
||||
}
|
||||
L30:
|
||||
jz = j1;
|
||||
ju = 0;
|
||||
|
||||
/* gaussian elimination with partial pivoting */
|
||||
|
||||
nm1 = *n - 1;
|
||||
if (nm1 < 1) {
|
||||
goto L130;
|
||||
}
|
||||
i__1 = nm1;
|
||||
for (k = 1; k <= i__1; ++k) {
|
||||
kp1 = k + 1;
|
||||
|
||||
/* zero next fill-in column */
|
||||
|
||||
++jz;
|
||||
if (jz > *n) {
|
||||
goto L50;
|
||||
}
|
||||
if (*ml < 1) {
|
||||
goto L50;
|
||||
}
|
||||
i__2 = *ml;
|
||||
for (i__ = 1; i__ <= i__2; ++i__) {
|
||||
abd[i__ + jz * abd_dim1] = 0.;
|
||||
/* L40: */
|
||||
}
|
||||
L50:
|
||||
|
||||
/* find l = pivot index */
|
||||
|
||||
/* Computing MIN */
|
||||
i__2 = *ml, i__3 = *n - k;
|
||||
lm = min(i__2,i__3);
|
||||
i__2 = lm + 1;
|
||||
l = idamax_(&i__2, &abd[m + k * abd_dim1], &c__1) + m - 1;
|
||||
ipvt[k] = l + k - m;
|
||||
|
||||
/* zero pivot implies this column already triangularized */
|
||||
|
||||
if (abd[l + k * abd_dim1] == 0.) {
|
||||
goto L100;
|
||||
}
|
||||
|
||||
/* interchange if necessary */
|
||||
|
||||
if (l == m) {
|
||||
goto L60;
|
||||
}
|
||||
t = abd[l + k * abd_dim1];
|
||||
abd[l + k * abd_dim1] = abd[m + k * abd_dim1];
|
||||
abd[m + k * abd_dim1] = t;
|
||||
L60:
|
||||
|
||||
/* compute multipliers */
|
||||
|
||||
t = -1. / abd[m + k * abd_dim1];
|
||||
dscal_(&lm, &t, &abd[m + 1 + k * abd_dim1], &c__1);
|
||||
|
||||
/* row elimination with column indexing */
|
||||
|
||||
/* Computing MIN */
|
||||
/* Computing MAX */
|
||||
i__3 = ju, i__4 = *mu + ipvt[k];
|
||||
i__2 = max(i__3,i__4);
|
||||
ju = min(i__2,*n);
|
||||
mm = m;
|
||||
if (ju < kp1) {
|
||||
goto L90;
|
||||
}
|
||||
i__2 = ju;
|
||||
for (j = kp1; j <= i__2; ++j) {
|
||||
--l;
|
||||
--mm;
|
||||
t = abd[l + j * abd_dim1];
|
||||
if (l == mm) {
|
||||
goto L70;
|
||||
}
|
||||
abd[l + j * abd_dim1] = abd[mm + j * abd_dim1];
|
||||
abd[mm + j * abd_dim1] = t;
|
||||
L70:
|
||||
daxpy_(&lm, &t, &abd[m + 1 + k * abd_dim1], &c__1, &abd[mm + 1 +
|
||||
j * abd_dim1], &c__1);
|
||||
/* L80: */
|
||||
}
|
||||
L90:
|
||||
goto L110;
|
||||
L100:
|
||||
*info = k;
|
||||
L110:
|
||||
/* L120: */
|
||||
;
|
||||
}
|
||||
L130:
|
||||
ipvt[*n] = *n;
|
||||
if (abd[m + *n * abd_dim1] == 0.) {
|
||||
*info = *n;
|
||||
}
|
||||
return 0;
|
||||
} /* dgbfa_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
206
ext/f2c_math/dgbsl.c
Normal file
206
ext/f2c_math/dgbsl.c
Normal file
|
|
@ -0,0 +1,206 @@
|
|||
/* dgbsl.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__1 = 1;
|
||||
|
||||
/* Subroutine */ int dgbsl_(doublereal *abd, integer *lda, integer *n,
|
||||
integer *ml, integer *mu, integer *ipvt, doublereal *b, integer *job)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer abd_dim1, abd_offset, i__1, i__2, i__3;
|
||||
|
||||
/* Local variables */
|
||||
static integer k, l, m;
|
||||
static doublereal t;
|
||||
static integer kb, la, lb, lm, nm1;
|
||||
extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
|
||||
integer *);
|
||||
extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
|
||||
integer *, doublereal *, integer *);
|
||||
|
||||
|
||||
/* dgbsl solves the double precision band system */
|
||||
/* a * x = b or trans(a) * x = b */
|
||||
/* using the factors computed by dgbco or dgbfa. */
|
||||
|
||||
/* on entry */
|
||||
|
||||
/* abd double precision(lda, n) */
|
||||
/* the output from dgbco or dgbfa. */
|
||||
|
||||
/* lda integer */
|
||||
/* the leading dimension of the array abd . */
|
||||
|
||||
/* n integer */
|
||||
/* the order of the original matrix. */
|
||||
|
||||
/* ml integer */
|
||||
/* number of diagonals below the main diagonal. */
|
||||
|
||||
/* mu integer */
|
||||
/* number of diagonals above the main diagonal. */
|
||||
|
||||
/* ipvt integer(n) */
|
||||
/* the pivot vector from dgbco or dgbfa. */
|
||||
|
||||
/* b double precision(n) */
|
||||
/* the right hand side vector. */
|
||||
|
||||
/* job integer */
|
||||
/* = 0 to solve a*x = b , */
|
||||
/* = nonzero to solve trans(a)*x = b , where */
|
||||
/* trans(a) is the transpose. */
|
||||
|
||||
/* on return */
|
||||
|
||||
/* b the solution vector x . */
|
||||
|
||||
/* error condition */
|
||||
|
||||
/* a division by zero will occur if the input factor contains a */
|
||||
/* zero on the diagonal. technically this indicates singularity */
|
||||
/* but it is often caused by improper arguments or improper */
|
||||
/* setting of lda . it will not occur if the subroutines are */
|
||||
/* called correctly and if dgbco has set rcond .gt. 0.0 */
|
||||
/* or dgbfa has set info .eq. 0 . */
|
||||
|
||||
/* to compute inverse(a) * c where c is a matrix */
|
||||
/* with p columns */
|
||||
/* call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z) */
|
||||
/* if (rcond is too small) go to ... */
|
||||
/* do 10 j = 1, p */
|
||||
/* call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0) */
|
||||
/* 10 continue */
|
||||
|
||||
/* linpack. this version dated 08/14/78 . */
|
||||
/* cleve moler, university of new mexico, argonne national lab. */
|
||||
|
||||
/* subroutines and functions */
|
||||
|
||||
/* blas daxpy,ddot */
|
||||
/* fortran min0 */
|
||||
|
||||
/* internal variables */
|
||||
|
||||
|
||||
/* Parameter adjustments */
|
||||
abd_dim1 = *lda;
|
||||
abd_offset = 1 + abd_dim1;
|
||||
abd -= abd_offset;
|
||||
--ipvt;
|
||||
--b;
|
||||
|
||||
/* Function Body */
|
||||
m = *mu + *ml + 1;
|
||||
nm1 = *n - 1;
|
||||
if (*job != 0) {
|
||||
goto L50;
|
||||
}
|
||||
|
||||
/* job = 0 , solve a * x = b */
|
||||
/* first solve l*y = b */
|
||||
|
||||
if (*ml == 0) {
|
||||
goto L30;
|
||||
}
|
||||
if (nm1 < 1) {
|
||||
goto L30;
|
||||
}
|
||||
i__1 = nm1;
|
||||
for (k = 1; k <= i__1; ++k) {
|
||||
/* Computing MIN */
|
||||
i__2 = *ml, i__3 = *n - k;
|
||||
lm = min(i__2,i__3);
|
||||
l = ipvt[k];
|
||||
t = b[l];
|
||||
if (l == k) {
|
||||
goto L10;
|
||||
}
|
||||
b[l] = b[k];
|
||||
b[k] = t;
|
||||
L10:
|
||||
daxpy_(&lm, &t, &abd[m + 1 + k * abd_dim1], &c__1, &b[k + 1], &c__1);
|
||||
/* L20: */
|
||||
}
|
||||
L30:
|
||||
|
||||
/* now solve u*x = y */
|
||||
|
||||
i__1 = *n;
|
||||
for (kb = 1; kb <= i__1; ++kb) {
|
||||
k = *n + 1 - kb;
|
||||
b[k] /= abd[m + k * abd_dim1];
|
||||
lm = min(k,m) - 1;
|
||||
la = m - lm;
|
||||
lb = k - lm;
|
||||
t = -b[k];
|
||||
daxpy_(&lm, &t, &abd[la + k * abd_dim1], &c__1, &b[lb], &c__1);
|
||||
/* L40: */
|
||||
}
|
||||
goto L100;
|
||||
L50:
|
||||
|
||||
/* job = nonzero, solve trans(a) * x = b */
|
||||
/* first solve trans(u)*y = b */
|
||||
|
||||
i__1 = *n;
|
||||
for (k = 1; k <= i__1; ++k) {
|
||||
lm = min(k,m) - 1;
|
||||
la = m - lm;
|
||||
lb = k - lm;
|
||||
t = ddot_(&lm, &abd[la + k * abd_dim1], &c__1, &b[lb], &c__1);
|
||||
b[k] = (b[k] - t) / abd[m + k * abd_dim1];
|
||||
/* L60: */
|
||||
}
|
||||
|
||||
/* now solve trans(l)*x = y */
|
||||
|
||||
if (*ml == 0) {
|
||||
goto L90;
|
||||
}
|
||||
if (nm1 < 1) {
|
||||
goto L90;
|
||||
}
|
||||
i__1 = nm1;
|
||||
for (kb = 1; kb <= i__1; ++kb) {
|
||||
k = *n - kb;
|
||||
/* Computing MIN */
|
||||
i__2 = *ml, i__3 = *n - k;
|
||||
lm = min(i__2,i__3);
|
||||
b[k] += ddot_(&lm, &abd[m + 1 + k * abd_dim1], &c__1, &b[k + 1], &
|
||||
c__1);
|
||||
l = ipvt[k];
|
||||
if (l == k) {
|
||||
goto L70;
|
||||
}
|
||||
t = b[l];
|
||||
b[l] = b[k];
|
||||
b[k] = t;
|
||||
L70:
|
||||
/* L80: */
|
||||
;
|
||||
}
|
||||
L90:
|
||||
L100:
|
||||
return 0;
|
||||
} /* dgbsl_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
164
ext/f2c_math/dgefa.c
Normal file
164
ext/f2c_math/dgefa.c
Normal file
|
|
@ -0,0 +1,164 @@
|
|||
/* dgefa.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__1 = 1;
|
||||
|
||||
/* Subroutine */ int dgefa_(doublereal *a, integer *lda, integer *n, integer *
|
||||
ipvt, integer *info)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, i__1, i__2, i__3;
|
||||
|
||||
/* Local variables */
|
||||
static integer j, k, l;
|
||||
static doublereal t;
|
||||
static integer kp1, nm1;
|
||||
extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
|
||||
integer *), daxpy_(integer *, doublereal *, doublereal *, integer
|
||||
*, doublereal *, integer *);
|
||||
extern integer idamax_(integer *, doublereal *, integer *);
|
||||
|
||||
|
||||
/* dgefa factors a double precision matrix by gaussian elimination. */
|
||||
|
||||
/* dgefa is usually called by dgeco, but it can be called */
|
||||
/* directly with a saving in time if rcond is not needed. */
|
||||
/* (time for dgeco) = (1 + 9/n)*(time for dgefa) . */
|
||||
|
||||
/* on entry */
|
||||
|
||||
/* a double precision(lda, n) */
|
||||
/* the matrix to be factored. */
|
||||
|
||||
/* lda integer */
|
||||
/* the leading dimension of the array a . */
|
||||
|
||||
/* n integer */
|
||||
/* the order of the matrix a . */
|
||||
|
||||
/* on return */
|
||||
|
||||
/* a an upper triangular matrix and the multipliers */
|
||||
/* which were used to obtain it. */
|
||||
/* the factorization can be written a = l*u where */
|
||||
/* l is a product of permutation and unit lower */
|
||||
/* triangular matrices and u is upper triangular. */
|
||||
|
||||
/* ipvt integer(n) */
|
||||
/* an integer vector of pivot indices. */
|
||||
|
||||
/* info integer */
|
||||
/* = 0 normal value. */
|
||||
/* = k if u(k,k) .eq. 0.0 . this is not an error */
|
||||
/* condition for this subroutine, but it does */
|
||||
/* indicate that dgesl or dgedi will divide by zero */
|
||||
/* if called. use rcond in dgeco for a reliable */
|
||||
/* indication of singularity. */
|
||||
|
||||
/* linpack. this version dated 08/14/78 . */
|
||||
/* cleve moler, university of new mexico, argonne national lab. */
|
||||
|
||||
/* subroutines and functions */
|
||||
|
||||
/* blas daxpy,dscal,idamax */
|
||||
|
||||
/* internal variables */
|
||||
|
||||
|
||||
|
||||
/* gaussian elimination with partial pivoting */
|
||||
|
||||
/* Parameter adjustments */
|
||||
a_dim1 = *lda;
|
||||
a_offset = 1 + a_dim1;
|
||||
a -= a_offset;
|
||||
--ipvt;
|
||||
|
||||
/* Function Body */
|
||||
*info = 0;
|
||||
nm1 = *n - 1;
|
||||
if (nm1 < 1) {
|
||||
goto L70;
|
||||
}
|
||||
i__1 = nm1;
|
||||
for (k = 1; k <= i__1; ++k) {
|
||||
kp1 = k + 1;
|
||||
|
||||
/* find l = pivot index */
|
||||
|
||||
i__2 = *n - k + 1;
|
||||
l = idamax_(&i__2, &a[k + k * a_dim1], &c__1) + k - 1;
|
||||
ipvt[k] = l;
|
||||
|
||||
/* zero pivot implies this column already triangularized */
|
||||
|
||||
if (a[l + k * a_dim1] == 0.) {
|
||||
goto L40;
|
||||
}
|
||||
|
||||
/* interchange if necessary */
|
||||
|
||||
if (l == k) {
|
||||
goto L10;
|
||||
}
|
||||
t = a[l + k * a_dim1];
|
||||
a[l + k * a_dim1] = a[k + k * a_dim1];
|
||||
a[k + k * a_dim1] = t;
|
||||
L10:
|
||||
|
||||
/* compute multipliers */
|
||||
|
||||
t = -1. / a[k + k * a_dim1];
|
||||
i__2 = *n - k;
|
||||
dscal_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1);
|
||||
|
||||
/* row elimination with column indexing */
|
||||
|
||||
i__2 = *n;
|
||||
for (j = kp1; j <= i__2; ++j) {
|
||||
t = a[l + j * a_dim1];
|
||||
if (l == k) {
|
||||
goto L20;
|
||||
}
|
||||
a[l + j * a_dim1] = a[k + j * a_dim1];
|
||||
a[k + j * a_dim1] = t;
|
||||
L20:
|
||||
i__3 = *n - k;
|
||||
daxpy_(&i__3, &t, &a[k + 1 + k * a_dim1], &c__1, &a[k + 1 + j *
|
||||
a_dim1], &c__1);
|
||||
/* L30: */
|
||||
}
|
||||
goto L50;
|
||||
L40:
|
||||
*info = k;
|
||||
L50:
|
||||
/* L60: */
|
||||
;
|
||||
}
|
||||
L70:
|
||||
ipvt[*n] = *n;
|
||||
if (a[*n + *n * a_dim1] == 0.) {
|
||||
*info = *n;
|
||||
}
|
||||
return 0;
|
||||
} /* dgefa_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
183
ext/f2c_math/dgesl.c
Normal file
183
ext/f2c_math/dgesl.c
Normal file
|
|
@ -0,0 +1,183 @@
|
|||
/* dgesl.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__1 = 1;
|
||||
|
||||
/* Subroutine */ int dgesl_(doublereal *a, integer *lda, integer *n, integer *
|
||||
ipvt, doublereal *b, integer *job)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, i__1, i__2;
|
||||
|
||||
/* Local variables */
|
||||
static integer k, l;
|
||||
static doublereal t;
|
||||
static integer kb, nm1;
|
||||
extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
|
||||
integer *);
|
||||
extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
|
||||
integer *, doublereal *, integer *);
|
||||
|
||||
|
||||
/* dgesl solves the double precision system */
|
||||
/* a * x = b or trans(a) * x = b */
|
||||
/* using the factors computed by dgeco or dgefa. */
|
||||
|
||||
/* on entry */
|
||||
|
||||
/* a double precision(lda, n) */
|
||||
/* the output from dgeco or dgefa. */
|
||||
|
||||
/* lda integer */
|
||||
/* the leading dimension of the array a . */
|
||||
|
||||
/* n integer */
|
||||
/* the order of the matrix a . */
|
||||
|
||||
/* ipvt integer(n) */
|
||||
/* the pivot vector from dgeco or dgefa. */
|
||||
|
||||
/* b double precision(n) */
|
||||
/* the right hand side vector. */
|
||||
|
||||
/* job integer */
|
||||
/* = 0 to solve a*x = b , */
|
||||
/* = nonzero to solve trans(a)*x = b where */
|
||||
/* trans(a) is the transpose. */
|
||||
|
||||
/* on return */
|
||||
|
||||
/* b the solution vector x . */
|
||||
|
||||
/* error condition */
|
||||
|
||||
/* a division by zero will occur if the input factor contains a */
|
||||
/* zero on the diagonal. technically this indicates singularity */
|
||||
/* but it is often caused by improper arguments or improper */
|
||||
/* setting of lda . it will not occur if the subroutines are */
|
||||
/* called correctly and if dgeco has set rcond .gt. 0.0 */
|
||||
/* or dgefa has set info .eq. 0 . */
|
||||
|
||||
/* to compute inverse(a) * c where c is a matrix */
|
||||
/* with p columns */
|
||||
/* call dgeco(a,lda,n,ipvt,rcond,z) */
|
||||
/* if (rcond is too small) go to ... */
|
||||
/* do 10 j = 1, p */
|
||||
/* call dgesl(a,lda,n,ipvt,c(1,j),0) */
|
||||
/* 10 continue */
|
||||
|
||||
/* linpack. this version dated 08/14/78 . */
|
||||
/* cleve moler, university of new mexico, argonne national lab. */
|
||||
|
||||
/* subroutines and functions */
|
||||
|
||||
/* blas daxpy,ddot */
|
||||
|
||||
/* internal variables */
|
||||
|
||||
|
||||
/* Parameter adjustments */
|
||||
a_dim1 = *lda;
|
||||
a_offset = 1 + a_dim1;
|
||||
a -= a_offset;
|
||||
--ipvt;
|
||||
--b;
|
||||
|
||||
/* Function Body */
|
||||
nm1 = *n - 1;
|
||||
if (*job != 0) {
|
||||
goto L50;
|
||||
}
|
||||
|
||||
/* job = 0 , solve a * x = b */
|
||||
/* first solve l*y = b */
|
||||
|
||||
if (nm1 < 1) {
|
||||
goto L30;
|
||||
}
|
||||
i__1 = nm1;
|
||||
for (k = 1; k <= i__1; ++k) {
|
||||
l = ipvt[k];
|
||||
t = b[l];
|
||||
if (l == k) {
|
||||
goto L10;
|
||||
}
|
||||
b[l] = b[k];
|
||||
b[k] = t;
|
||||
L10:
|
||||
i__2 = *n - k;
|
||||
daxpy_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
|
||||
/* L20: */
|
||||
}
|
||||
L30:
|
||||
|
||||
/* now solve u*x = y */
|
||||
|
||||
i__1 = *n;
|
||||
for (kb = 1; kb <= i__1; ++kb) {
|
||||
k = *n + 1 - kb;
|
||||
b[k] /= a[k + k * a_dim1];
|
||||
t = -b[k];
|
||||
i__2 = k - 1;
|
||||
daxpy_(&i__2, &t, &a[k * a_dim1 + 1], &c__1, &b[1], &c__1);
|
||||
/* L40: */
|
||||
}
|
||||
goto L100;
|
||||
L50:
|
||||
|
||||
/* job = nonzero, solve trans(a) * x = b */
|
||||
/* first solve trans(u)*y = b */
|
||||
|
||||
i__1 = *n;
|
||||
for (k = 1; k <= i__1; ++k) {
|
||||
i__2 = k - 1;
|
||||
t = ddot_(&i__2, &a[k * a_dim1 + 1], &c__1, &b[1], &c__1);
|
||||
b[k] = (b[k] - t) / a[k + k * a_dim1];
|
||||
/* L60: */
|
||||
}
|
||||
|
||||
/* now solve trans(l)*x = y */
|
||||
|
||||
if (nm1 < 1) {
|
||||
goto L90;
|
||||
}
|
||||
i__1 = nm1;
|
||||
for (kb = 1; kb <= i__1; ++kb) {
|
||||
k = *n - kb;
|
||||
i__2 = *n - k;
|
||||
b[k] += ddot_(&i__2, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
|
||||
l = ipvt[k];
|
||||
if (l == k) {
|
||||
goto L70;
|
||||
}
|
||||
t = b[l];
|
||||
b[l] = b[k];
|
||||
b[k] = t;
|
||||
L70:
|
||||
/* L80: */
|
||||
;
|
||||
}
|
||||
L90:
|
||||
L100:
|
||||
return 0;
|
||||
} /* dgesl_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
255
ext/f2c_math/dp1vlu.c
Normal file
255
ext/f2c_math/dp1vlu.c
Normal file
|
|
@ -0,0 +1,255 @@
|
|||
/* dp1vlu.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__1 = 1;
|
||||
static integer c__8 = 8;
|
||||
static integer c__2 = 2;
|
||||
static integer c__5 = 5;
|
||||
|
||||
/* DECK DP1VLU */
|
||||
/* Subroutine */ int dp1vlu_(integer *l, integer *nder, doublereal *x,
|
||||
doublereal *yfit, doublereal *yp, doublereal *a)
|
||||
{
|
||||
/* System generated locals */
|
||||
address a__1[5];
|
||||
integer i__1, i__2, i__3[5];
|
||||
char ch__1[150];
|
||||
|
||||
/* Builtin functions */
|
||||
integer s_wsfi(icilist *), do_fio(integer *, char *, ftnlen), e_wsfi(void)
|
||||
;
|
||||
/* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, n, k1, k2, k3, k4;
|
||||
static doublereal cc;
|
||||
static integer ic, kc, in, k1i, lm1, lp1;
|
||||
static doublereal dif;
|
||||
static integer k3p1, k4p1, ndo;
|
||||
static doublereal val;
|
||||
static integer ilo, iup, ndp1, inp1, k3pn, k4pn, nord;
|
||||
static char xern1[8], xern2[8];
|
||||
static integer maxord;
|
||||
extern /* Subroutine */ int xermsg_(char *, char *, char *, integer *,
|
||||
integer *, ftnlen, ftnlen, ftnlen);
|
||||
|
||||
/* Fortran I/O blocks */
|
||||
static icilist io___28 = { 0, xern1, 0, "(I8)", 8, 1 };
|
||||
static icilist io___30 = { 0, xern2, 0, "(I8)", 8, 1 };
|
||||
|
||||
|
||||
/* ***BEGIN PROLOGUE DP1VLU */
|
||||
/* ***PURPOSE Use the coefficients generated by DPOLFT to evaluate the */
|
||||
/* polynomial fit of degree L, along with the first NDER of */
|
||||
/* its derivatives, at a specified point. */
|
||||
/* ***LIBRARY SLATEC */
|
||||
/* ***CATEGORY K6 */
|
||||
/* ***TYPE DOUBLE PRECISION (PVALUE-S, DP1VLU-D) */
|
||||
/* ***KEYWORDS CURVE FITTING, LEAST SQUARES, POLYNOMIAL APPROXIMATION */
|
||||
/* ***AUTHOR Shampine, L. F., (SNLA) */
|
||||
/* Davenport, S. M., (SNLA) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* Abstract */
|
||||
|
||||
/* The subroutine DP1VLU uses the coefficients generated by DPOLFT */
|
||||
/* to evaluate the polynomial fit of degree L , along with the first */
|
||||
/* NDER of its derivatives, at a specified point. Computationally */
|
||||
/* stable recurrence relations are used to perform this task. */
|
||||
|
||||
/* The parameters for DP1VLU are */
|
||||
|
||||
/* Input -- ALL TYPE REAL variables are DOUBLE PRECISION */
|
||||
/* L - the degree of polynomial to be evaluated. L may be */
|
||||
/* any non-negative integer which is less than or equal */
|
||||
/* to NDEG , the highest degree polynomial provided */
|
||||
/* by DPOLFT . */
|
||||
/* NDER - the number of derivatives to be evaluated. NDER */
|
||||
/* may be 0 or any positive value. If NDER is less */
|
||||
/* than 0, it will be treated as 0. */
|
||||
/* X - the argument at which the polynomial and its */
|
||||
/* derivatives are to be evaluated. */
|
||||
/* A - work and output array containing values from last */
|
||||
/* call to DPOLFT . */
|
||||
|
||||
/* Output -- ALL TYPE REAL variables are DOUBLE PRECISION */
|
||||
/* YFIT - value of the fitting polynomial of degree L at X */
|
||||
/* YP - array containing the first through NDER derivatives */
|
||||
/* of the polynomial of degree L . YP must be */
|
||||
/* dimensioned at least NDER in the calling program. */
|
||||
|
||||
/* ***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, */
|
||||
/* Curve fitting by polynomials in one variable, Report */
|
||||
/* SLA-74-0270, Sandia Laboratories, June 1974. */
|
||||
/* ***ROUTINES CALLED XERMSG */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 740601 DATE WRITTEN */
|
||||
/* 890531 Changed all specific intrinsics to generic. (WRB) */
|
||||
/* 890911 Removed unnecessary intrinsics. (WRB) */
|
||||
/* 891006 Cosmetic changes to prologue. (WRB) */
|
||||
/* 891006 REVISION DATE from Version 3.2 */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
|
||||
/* 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) */
|
||||
/* 900510 Convert XERRWV calls to XERMSG calls. (RWC) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* ***END PROLOGUE DP1VLU */
|
||||
/* ***FIRST EXECUTABLE STATEMENT DP1VLU */
|
||||
/* Parameter adjustments */
|
||||
--a;
|
||||
--yp;
|
||||
|
||||
/* Function Body */
|
||||
if (*l < 0) {
|
||||
goto L12;
|
||||
}
|
||||
ndo = max(*nder,0);
|
||||
ndo = min(ndo,*l);
|
||||
maxord = (integer) (a[1] + .5);
|
||||
k1 = maxord + 1;
|
||||
k2 = k1 + maxord;
|
||||
k3 = k2 + maxord + 2;
|
||||
nord = (integer) (a[k3] + .5);
|
||||
if (*l > nord) {
|
||||
goto L11;
|
||||
}
|
||||
k4 = k3 + *l + 1;
|
||||
if (*nder < 1) {
|
||||
goto L2;
|
||||
}
|
||||
i__1 = *nder;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
/* L1: */
|
||||
yp[i__] = 0.;
|
||||
}
|
||||
L2:
|
||||
if (*l >= 2) {
|
||||
goto L4;
|
||||
}
|
||||
if (*l == 1) {
|
||||
goto L3;
|
||||
}
|
||||
|
||||
/* L IS 0 */
|
||||
|
||||
val = a[k2 + 1];
|
||||
goto L10;
|
||||
|
||||
/* L IS 1 */
|
||||
|
||||
L3:
|
||||
cc = a[k2 + 2];
|
||||
val = a[k2 + 1] + (*x - a[2]) * cc;
|
||||
if (*nder >= 1) {
|
||||
yp[1] = cc;
|
||||
}
|
||||
goto L10;
|
||||
|
||||
/* L IS GREATER THAN 1 */
|
||||
|
||||
L4:
|
||||
ndp1 = ndo + 1;
|
||||
k3p1 = k3 + 1;
|
||||
k4p1 = k4 + 1;
|
||||
lp1 = *l + 1;
|
||||
lm1 = *l - 1;
|
||||
ilo = k3 + 3;
|
||||
iup = k4 + ndp1;
|
||||
i__1 = iup;
|
||||
for (i__ = ilo; i__ <= i__1; ++i__) {
|
||||
/* L5: */
|
||||
a[i__] = 0.;
|
||||
}
|
||||
dif = *x - a[lp1];
|
||||
kc = k2 + lp1;
|
||||
a[k4p1] = a[kc];
|
||||
a[k3p1] = a[kc - 1] + dif * a[k4p1];
|
||||
a[k3 + 2] = a[k4p1];
|
||||
|
||||
/* EVALUATE RECURRENCE RELATIONS FOR FUNCTION VALUE AND DERIVATIVES */
|
||||
|
||||
i__1 = lm1;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
in = *l - i__;
|
||||
inp1 = in + 1;
|
||||
k1i = k1 + inp1;
|
||||
ic = k2 + in;
|
||||
dif = *x - a[inp1];
|
||||
val = a[ic] + dif * a[k3p1] - a[k1i] * a[k4p1];
|
||||
if (ndo <= 0) {
|
||||
goto L8;
|
||||
}
|
||||
i__2 = ndo;
|
||||
for (n = 1; n <= i__2; ++n) {
|
||||
k3pn = k3p1 + n;
|
||||
k4pn = k4p1 + n;
|
||||
/* L6: */
|
||||
yp[n] = dif * a[k3pn] + n * a[k3pn - 1] - a[k1i] * a[k4pn];
|
||||
}
|
||||
|
||||
/* SAVE VALUES NEEDED FOR NEXT EVALUATION OF RECURRENCE RELATIONS */
|
||||
|
||||
i__2 = ndo;
|
||||
for (n = 1; n <= i__2; ++n) {
|
||||
k3pn = k3p1 + n;
|
||||
k4pn = k4p1 + n;
|
||||
a[k4pn] = a[k3pn];
|
||||
/* L7: */
|
||||
a[k3pn] = yp[n];
|
||||
}
|
||||
L8:
|
||||
a[k4p1] = a[k3p1];
|
||||
/* L9: */
|
||||
a[k3p1] = val;
|
||||
}
|
||||
|
||||
/* NORMAL RETURN OR ABORT DUE TO ERROR */
|
||||
|
||||
L10:
|
||||
*yfit = val;
|
||||
return 0;
|
||||
|
||||
L11:
|
||||
s_wsfi(&io___28);
|
||||
do_fio(&c__1, (char *)&(*l), (ftnlen)sizeof(integer));
|
||||
e_wsfi();
|
||||
s_wsfi(&io___30);
|
||||
do_fio(&c__1, (char *)&nord, (ftnlen)sizeof(integer));
|
||||
e_wsfi();
|
||||
/* Writing concatenation */
|
||||
i__3[0] = 40, a__1[0] = "THE ORDER OF POLYNOMIAL EVALUATION, L = ";
|
||||
i__3[1] = 8, a__1[1] = xern1;
|
||||
i__3[2] = 49, a__1[2] = " REQUESTED EXCEEDS THE HIGHEST ORDER FIT, NORD "
|
||||
"= ";
|
||||
i__3[3] = 8, a__1[3] = xern2;
|
||||
i__3[4] = 45, a__1[4] = ", COMPUTED BY DPOLFT -- EXECUTION TERMINATED.";
|
||||
s_cat(ch__1, a__1, i__3, &c__5, (ftnlen)150);
|
||||
xermsg_("SLATEC", "DP1VLU", ch__1, &c__8, &c__2, (ftnlen)6, (ftnlen)6, (
|
||||
ftnlen)150);
|
||||
return 0;
|
||||
|
||||
L12:
|
||||
xermsg_("SLATEC", "DP1VLU", "INVALID INPUT PARAMETER. ORDER OF POLYNOMI"
|
||||
"AL EVALUATION REQUESTED IS NEGATIVE.", &c__2, &c__2, (ftnlen)6, (
|
||||
ftnlen)6, (ftnlen)79);
|
||||
return 0;
|
||||
} /* dp1vlu_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
127
ext/f2c_math/dpcoef.c
Normal file
127
ext/f2c_math/dpcoef.c
Normal file
|
|
@ -0,0 +1,127 @@
|
|||
/* dpcoef.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* DECK DPCOEF */
|
||||
/* Subroutine */ int dpcoef_(integer *l, doublereal *c__, doublereal *tc,
|
||||
doublereal *a)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer i__1;
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, ll, nr;
|
||||
static doublereal fac;
|
||||
static integer new__, llp1, llp2;
|
||||
static doublereal save;
|
||||
extern /* Subroutine */ int dp1vlu_(integer *, integer *, doublereal *,
|
||||
doublereal *, doublereal *, doublereal *);
|
||||
|
||||
/* ***BEGIN PROLOGUE DPCOEF */
|
||||
/* ***PURPOSE Convert the DPOLFT coefficients to Taylor series form. */
|
||||
/* ***LIBRARY SLATEC */
|
||||
/* ***CATEGORY K1A1A2 */
|
||||
/* ***TYPE DOUBLE PRECISION (PCOEF-S, DPCOEF-D) */
|
||||
/* ***KEYWORDS CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT */
|
||||
/* ***AUTHOR Shampine, L. F., (SNLA) */
|
||||
/* Davenport, S. M., (SNLA) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* Abstract */
|
||||
|
||||
/* DPOLFT computes the least squares polynomial fit of degree L as */
|
||||
/* a sum of orthogonal polynomials. DPCOEF changes this fit to its */
|
||||
/* Taylor expansion about any point C , i.e. writes the polynomial */
|
||||
/* as a sum of powers of (X-C). Taking C=0. gives the polynomial */
|
||||
/* in powers of X, but a suitable non-zero C often leads to */
|
||||
/* polynomials which are better scaled and more accurately evaluated. */
|
||||
|
||||
/* The parameters for DPCOEF are */
|
||||
|
||||
/* INPUT -- All TYPE REAL variables are DOUBLE PRECISION */
|
||||
/* L - Indicates the degree of polynomial to be changed to */
|
||||
/* its Taylor expansion. To obtain the Taylor */
|
||||
/* coefficients in reverse order, input L as the */
|
||||
/* negative of the degree desired. The absolute value */
|
||||
/* of L must be less than or equal to NDEG, the highest */
|
||||
/* degree polynomial fitted by DPOLFT . */
|
||||
/* C - The point about which the Taylor expansion is to be */
|
||||
/* made. */
|
||||
/* A - Work and output array containing values from last */
|
||||
/* call to DPOLFT . */
|
||||
|
||||
/* OUTPUT -- All TYPE REAL variables are DOUBLE PRECISION */
|
||||
/* TC - Vector containing the first LL+1 Taylor coefficients */
|
||||
/* where LL=ABS(L). If L.GT.0 , the coefficients are */
|
||||
/* in the usual Taylor series order, i.e. */
|
||||
/* P(X) = TC(1) + TC(2)*(X-C) + ... + TC(N+1)*(X-C)**N */
|
||||
/* If L .LT. 0, the coefficients are in reverse order, */
|
||||
/* i.e. */
|
||||
/* P(X) = TC(1)*(X-C)**N + ... + TC(N)*(X-C) + TC(N+1) */
|
||||
|
||||
/* ***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, */
|
||||
/* Curve fitting by polynomials in one variable, Report */
|
||||
/* SLA-74-0270, Sandia Laboratories, June 1974. */
|
||||
/* ***ROUTINES CALLED DP1VLU */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 740601 DATE WRITTEN */
|
||||
/* 890531 Changed all specific intrinsics to generic. (WRB) */
|
||||
/* 891006 Cosmetic changes to prologue. (WRB) */
|
||||
/* 891006 REVISION DATE from Version 3.2 */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* ***END PROLOGUE DPCOEF */
|
||||
|
||||
/* ***FIRST EXECUTABLE STATEMENT DPCOEF */
|
||||
/* Parameter adjustments */
|
||||
--a;
|
||||
--tc;
|
||||
|
||||
/* Function Body */
|
||||
ll = abs(*l);
|
||||
llp1 = ll + 1;
|
||||
dp1vlu_(&ll, &ll, c__, &tc[1], &tc[2], &a[1]);
|
||||
if (ll < 2) {
|
||||
goto L2;
|
||||
}
|
||||
fac = 1.;
|
||||
i__1 = llp1;
|
||||
for (i__ = 3; i__ <= i__1; ++i__) {
|
||||
fac *= i__ - 1;
|
||||
/* L1: */
|
||||
tc[i__] /= fac;
|
||||
}
|
||||
L2:
|
||||
if (*l >= 0) {
|
||||
goto L4;
|
||||
}
|
||||
nr = llp1 / 2;
|
||||
llp2 = ll + 2;
|
||||
i__1 = nr;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
save = tc[i__];
|
||||
new__ = llp2 - i__;
|
||||
tc[i__] = tc[new__];
|
||||
/* L3: */
|
||||
tc[new__] = save;
|
||||
}
|
||||
L4:
|
||||
return 0;
|
||||
} /* dpcoef_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
539
ext/f2c_math/dpolft.c
Normal file
539
ext/f2c_math/dpolft.c
Normal file
|
|
@ -0,0 +1,539 @@
|
|||
/* dpolft.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__2 = 2;
|
||||
static integer c__1 = 1;
|
||||
|
||||
/* DECK DPOLFT */
|
||||
/* Subroutine */ int dpolft_(integer *n, doublereal *x, doublereal *y,
|
||||
doublereal *w, integer *maxdeg, integer *ndeg, doublereal *eps,
|
||||
doublereal *r__, integer *ierr, doublereal *a)
|
||||
{
|
||||
/* Initialized data */
|
||||
|
||||
static doublereal co[12] /* was [4][3] */ = { -13.08685,-2.4648165,
|
||||
-3.3846535,-1.2973162,-3.3381146,-1.7812271,-3.2578406,-1.6589279,
|
||||
-1.6282703,-1.3152745,-3.2640179,-1.9829776 };
|
||||
|
||||
/* System generated locals */
|
||||
integer i__1;
|
||||
doublereal d__1;
|
||||
|
||||
/* Builtin functions */
|
||||
double sqrt(doublereal);
|
||||
|
||||
/* Local variables */
|
||||
static doublereal f;
|
||||
static integer i__, j, m, k1, k2, k3, k4, k5;
|
||||
static doublereal w1, w11, xm, yp;
|
||||
static integer jp1;
|
||||
static doublereal den, sig;
|
||||
static integer k1pj, k2pj, k3pi, k4pi, k5pi, mop1;
|
||||
static doublereal degf;
|
||||
static integer nder;
|
||||
static doublereal sigj;
|
||||
static integer jpas, ksig;
|
||||
static doublereal temp, etst, temd1, temd2;
|
||||
static integer idegf, nfail;
|
||||
static doublereal fcrit, sigjm1;
|
||||
extern /* Subroutine */ int dp1vlu_(integer *, integer *, doublereal *,
|
||||
doublereal *, doublereal *, doublereal *);
|
||||
static doublereal sigpas;
|
||||
extern /* Subroutine */ int xermsg_(char *, char *, char *, integer *,
|
||||
integer *, ftnlen, ftnlen, ftnlen);
|
||||
|
||||
/* ***BEGIN PROLOGUE DPOLFT */
|
||||
/* ***PURPOSE Fit discrete data in a least squares sense by polynomials */
|
||||
/* in one variable. */
|
||||
/* ***LIBRARY SLATEC */
|
||||
/* ***CATEGORY K1A1A2 */
|
||||
/* ***TYPE DOUBLE PRECISION (POLFIT-S, DPOLFT-D) */
|
||||
/* ***KEYWORDS CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT */
|
||||
/* ***AUTHOR Shampine, L. F., (SNLA) */
|
||||
/* Davenport, S. M., (SNLA) */
|
||||
/* Huddleston, R. E., (SNLL) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* Abstract */
|
||||
|
||||
/* 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. */
|
||||
|
||||
/* The parameters for DPOLFT are */
|
||||
|
||||
/* Input -- All TYPE REAL variables are DOUBLE PRECISION */
|
||||
/* N - the number of data points. The arrays X, Y and W */
|
||||
/* must be dimensioned at least N (N .GE. 1). */
|
||||
/* X - array of values of the independent variable. These */
|
||||
/* values may appear in any order and need not all be */
|
||||
/* distinct. */
|
||||
/* Y - array of corresponding function values. */
|
||||
/* W - array of positive values to be used as weights. If */
|
||||
/* W(1) 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 . */
|
||||
/* 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. */
|
||||
/* 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. */
|
||||
|
||||
/* Output -- All TYPE REAL variables are DOUBLE PRECISION */
|
||||
/* NDEG - degree of the highest degree fit computed. */
|
||||
/* EPS - RMS error of the polynomial of degree NDEG . */
|
||||
/* R - vector of dimension at least NDEG containing values */
|
||||
/* of the fit of degree NDEG at each of the X(I) . */
|
||||
/* Except when the statistical test is used, these */
|
||||
/* values are more accurate than results from subroutine */
|
||||
/* DP1VLU normally are. */
|
||||
/* IERR - error flag with the following possible values. */
|
||||
/* 1 -- indicates normal execution, i.e., either */
|
||||
/* (1) the input value of EPS was negative, and the */
|
||||
/* computed polynomial fit of degree NDEG */
|
||||
/* satisfies the specified F test, or */
|
||||
/* (2) the input value of EPS was 0., and the fits of */
|
||||
/* all degrees up to MAXDEG are complete, or */
|
||||
/* (3) the input value of EPS was positive, and the */
|
||||
/* polynomial of degree NDEG satisfies the RMS */
|
||||
/* error requirement. */
|
||||
/* 2 -- invalid input parameter. At least one of the input */
|
||||
/* parameters has an illegal value and must be corrected */
|
||||
/* before DPOLFT can proceed. Valid input results */
|
||||
/* when the following restrictions are observed */
|
||||
/* N .GE. 1 */
|
||||
/* 0 .LE. MAXDEG .LE. N-1 for EPS .GE. 0. */
|
||||
/* 0 .LE. MAXDEG .LE. N-2 for EPS .LT. 0. */
|
||||
/* W(1)=-1.0 or W(I) .GT. 0., I=1,...,N . */
|
||||
/* 3 -- cannot satisfy the RMS error requirement with a */
|
||||
/* polynomial of degree no greater than MAXDEG . Best */
|
||||
/* fit found is of degree MAXDEG . */
|
||||
/* 4 -- cannot satisfy the test for significance using */
|
||||
/* current value of MAXDEG . Statistically, the */
|
||||
/* best fit found is of order NORD . (In this case, */
|
||||
/* NDEG will have one of the values: MAXDEG-2, */
|
||||
/* MAXDEG-1, or MAXDEG). Using a higher value of */
|
||||
/* MAXDEG may result in passing the test. */
|
||||
/* A - work and output array having at least 3N+3MAXDEG+3 */
|
||||
/* locations */
|
||||
|
||||
/* Note - DPOLFT calculates all fits of degrees up to and including */
|
||||
/* NDEG . Any or all of these fits can be evaluated or */
|
||||
/* expressed as powers of (X-C) using DP1VLU and DPCOEF */
|
||||
/* after just one call to DPOLFT . */
|
||||
|
||||
/* ***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, */
|
||||
/* Curve fitting by polynomials in one variable, Report */
|
||||
/* SLA-74-0270, Sandia Laboratories, June 1974. */
|
||||
/* ***ROUTINES CALLED DP1VLU, XERMSG */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 740601 DATE WRITTEN */
|
||||
/* 890531 Changed all specific intrinsics to generic. (WRB) */
|
||||
/* 891006 Cosmetic changes to prologue. (WRB) */
|
||||
/* 891006 REVISION DATE from Version 3.2 */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
|
||||
/* 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) */
|
||||
/* 900911 Added variable YP to DOUBLE PRECISION declaration. (WRB) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* 920527 Corrected erroneous statements in DESCRIPTION. (WRB) */
|
||||
/* ***END PROLOGUE DPOLFT */
|
||||
/* Parameter adjustments */
|
||||
--a;
|
||||
--r__;
|
||||
--w;
|
||||
--y;
|
||||
--x;
|
||||
|
||||
/* Function Body */
|
||||
/* ***FIRST EXECUTABLE STATEMENT DPOLFT */
|
||||
/* write(*,*) 'DPOLFT n = ',n */
|
||||
/* do ii = 1,n */
|
||||
/* write(*,*) x(ii), y(ii), w(ii) */
|
||||
/* end do */
|
||||
/* write(*,*) ' maxdeg, eps = ',maxdeg,eps */
|
||||
m = abs(*n);
|
||||
if (m == 0) {
|
||||
goto L30;
|
||||
}
|
||||
if (*maxdeg < 0) {
|
||||
goto L30;
|
||||
}
|
||||
a[1] = (doublereal) (*maxdeg);
|
||||
mop1 = *maxdeg + 1;
|
||||
if (m < mop1) {
|
||||
goto L30;
|
||||
}
|
||||
if (*eps < 0. && m == mop1) {
|
||||
goto L30;
|
||||
}
|
||||
xm = (doublereal) m;
|
||||
etst = *eps * *eps * xm;
|
||||
if (w[1] < 0.) {
|
||||
goto L2;
|
||||
}
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
if (w[i__] <= 0.) {
|
||||
goto L30;
|
||||
}
|
||||
/* L1: */
|
||||
}
|
||||
goto L4;
|
||||
L2:
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
/* L3: */
|
||||
w[i__] = 1.;
|
||||
}
|
||||
L4:
|
||||
if (*eps >= 0.) {
|
||||
goto L8;
|
||||
}
|
||||
|
||||
/* DETERMINE SIGNIFICANCE LEVEL INDEX TO BE USED IN STATISTICAL TEST FOR */
|
||||
/* CHOOSING DEGREE OF POLYNOMIAL FIT */
|
||||
|
||||
if (*eps > -.55) {
|
||||
goto L5;
|
||||
}
|
||||
idegf = m - *maxdeg - 1;
|
||||
ksig = 1;
|
||||
if (idegf < 10) {
|
||||
ksig = 2;
|
||||
}
|
||||
if (idegf < 5) {
|
||||
ksig = 3;
|
||||
}
|
||||
goto L8;
|
||||
L5:
|
||||
ksig = 1;
|
||||
if (*eps < -.03) {
|
||||
ksig = 2;
|
||||
}
|
||||
if (*eps < -.07) {
|
||||
ksig = 3;
|
||||
}
|
||||
|
||||
/* INITIALIZE INDEXES AND COEFFICIENTS FOR FITTING */
|
||||
|
||||
L8:
|
||||
k1 = *maxdeg + 1;
|
||||
k2 = k1 + *maxdeg;
|
||||
k3 = k2 + *maxdeg + 2;
|
||||
k4 = k3 + m;
|
||||
k5 = k4 + m;
|
||||
i__1 = k4;
|
||||
for (i__ = 2; i__ <= i__1; ++i__) {
|
||||
/* L9: */
|
||||
a[i__] = 0.;
|
||||
}
|
||||
w11 = 0.;
|
||||
if (*n < 0) {
|
||||
goto L11;
|
||||
}
|
||||
|
||||
/* UNCONSTRAINED CASE */
|
||||
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
a[k4pi] = 1.;
|
||||
/* L10: */
|
||||
w11 += w[i__];
|
||||
}
|
||||
goto L13;
|
||||
|
||||
/* CONSTRAINED CASE */
|
||||
|
||||
L11:
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
/* L12: */
|
||||
/* Computing 2nd power */
|
||||
d__1 = a[k4pi];
|
||||
w11 += w[i__] * (d__1 * d__1);
|
||||
}
|
||||
|
||||
/* COMPUTE FIT OF DEGREE ZERO */
|
||||
|
||||
L13:
|
||||
temd1 = 0.;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
temd1 += w[i__] * y[i__] * a[k4pi];
|
||||
/* L14: */
|
||||
}
|
||||
temd1 /= w11;
|
||||
a[k2 + 1] = temd1;
|
||||
sigj = 0.;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
k5pi = k5 + i__;
|
||||
temd2 = temd1 * a[k4pi];
|
||||
r__[i__] = temd2;
|
||||
a[k5pi] = temd2 - r__[i__];
|
||||
/* L15: */
|
||||
/* Computing 2nd power */
|
||||
d__1 = y[i__] - r__[i__] - a[k5pi];
|
||||
sigj += w[i__] * (d__1 * d__1);
|
||||
}
|
||||
j = 0;
|
||||
|
||||
/* SEE IF POLYNOMIAL OF DEGREE 0 SATISFIES THE DEGREE SELECTION CRITERION */
|
||||
|
||||
if (*eps < 0.) {
|
||||
goto L24;
|
||||
} else if (*eps == 0) {
|
||||
goto L26;
|
||||
} else {
|
||||
goto L27;
|
||||
}
|
||||
|
||||
/* INCREMENT DEGREE */
|
||||
|
||||
L16:
|
||||
++j;
|
||||
jp1 = j + 1;
|
||||
k1pj = k1 + j;
|
||||
k2pj = k2 + j;
|
||||
sigjm1 = sigj;
|
||||
|
||||
/* COMPUTE NEW B COEFFICIENT EXCEPT WHEN J = 1 */
|
||||
|
||||
if (j > 1) {
|
||||
a[k1pj] = w11 / w1;
|
||||
}
|
||||
|
||||
/* COMPUTE NEW A COEFFICIENT */
|
||||
|
||||
temd1 = 0.;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
temd2 = a[k4pi];
|
||||
temd1 += x[i__] * w[i__] * temd2 * temd2;
|
||||
/* L18: */
|
||||
}
|
||||
a[jp1] = temd1 / w11;
|
||||
|
||||
/* EVALUATE ORTHOGONAL POLYNOMIAL AT DATA POINTS */
|
||||
|
||||
w1 = w11;
|
||||
w11 = 0.;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k3pi = k3 + i__;
|
||||
k4pi = k4 + i__;
|
||||
temp = a[k3pi];
|
||||
a[k3pi] = a[k4pi];
|
||||
a[k4pi] = (x[i__] - a[jp1]) * a[k3pi] - a[k1pj] * temp;
|
||||
/* L19: */
|
||||
/* Computing 2nd power */
|
||||
d__1 = a[k4pi];
|
||||
w11 += w[i__] * (d__1 * d__1);
|
||||
}
|
||||
|
||||
/* GET NEW ORTHOGONAL POLYNOMIAL COEFFICIENT USING PARTIAL DOUBLE */
|
||||
/* PRECISION */
|
||||
|
||||
temd1 = 0.;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
k5pi = k5 + i__;
|
||||
temd2 = w[i__] * (y[i__] - r__[i__] - a[k5pi]) * a[k4pi];
|
||||
/* L20: */
|
||||
temd1 += temd2;
|
||||
}
|
||||
temd1 /= w11;
|
||||
a[k2pj + 1] = temd1;
|
||||
|
||||
/* UPDATE POLYNOMIAL EVALUATIONS AT EACH OF THE DATA POINTS, AND */
|
||||
/* ACCUMULATE SUM OF SQUARES OF ERRORS. THE POLYNOMIAL EVALUATIONS ARE */
|
||||
/* COMPUTED AND STORED IN EXTENDED PRECISION. FOR THE I-TH DATA POINT, */
|
||||
/* THE MOST SIGNIFICANT BITS ARE STORED IN R(I) , AND THE LEAST */
|
||||
/* SIGNIFICANT BITS ARE IN A(K5PI) . */
|
||||
|
||||
sigj = 0.;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
k5pi = k5 + i__;
|
||||
temd2 = r__[i__] + a[k5pi] + temd1 * a[k4pi];
|
||||
r__[i__] = temd2;
|
||||
a[k5pi] = temd2 - r__[i__];
|
||||
/* L21: */
|
||||
/* Computing 2nd power */
|
||||
d__1 = y[i__] - r__[i__] - a[k5pi];
|
||||
sigj += w[i__] * (d__1 * d__1);
|
||||
}
|
||||
|
||||
/* SEE IF DEGREE SELECTION CRITERION HAS BEEN SATISFIED OR IF DEGREE */
|
||||
/* MAXDEG HAS BEEN REACHED */
|
||||
|
||||
if (*eps < 0.) {
|
||||
goto L23;
|
||||
} else if (*eps == 0) {
|
||||
goto L26;
|
||||
} else {
|
||||
goto L27;
|
||||
}
|
||||
|
||||
/* COMPUTE F STATISTICS (INPUT EPS .LT. 0.) */
|
||||
|
||||
L23:
|
||||
if (sigj == 0.) {
|
||||
goto L29;
|
||||
}
|
||||
degf = (doublereal) (m - j - 1);
|
||||
den = (co[(ksig << 2) - 1] * degf + 1.) * degf;
|
||||
fcrit = ((co[(ksig << 2) - 2] * degf + co[(ksig << 2) - 3]) * degf + co[(
|
||||
ksig << 2) - 4]) / den;
|
||||
fcrit *= fcrit;
|
||||
f = (sigjm1 - sigj) * degf / sigj;
|
||||
if (f < fcrit) {
|
||||
goto L25;
|
||||
}
|
||||
|
||||
/* POLYNOMIAL OF DEGREE J SATISFIES F TEST */
|
||||
|
||||
L24:
|
||||
sigpas = sigj;
|
||||
jpas = j;
|
||||
nfail = 0;
|
||||
if (*maxdeg == j) {
|
||||
goto L32;
|
||||
}
|
||||
goto L16;
|
||||
|
||||
/* POLYNOMIAL OF DEGREE J FAILS F TEST. IF THERE HAVE BEEN THREE */
|
||||
/* SUCCESSIVE FAILURES, A STATISTICALLY BEST DEGREE HAS BEEN FOUND. */
|
||||
|
||||
L25:
|
||||
++nfail;
|
||||
if (nfail >= 3) {
|
||||
goto L29;
|
||||
}
|
||||
if (*maxdeg == j) {
|
||||
goto L32;
|
||||
}
|
||||
goto L16;
|
||||
|
||||
/* RAISE THE DEGREE IF DEGREE MAXDEG HAS NOT YET BEEN REACHED (INPUT */
|
||||
/* EPS = 0.) */
|
||||
|
||||
L26:
|
||||
if (*maxdeg == j) {
|
||||
goto L28;
|
||||
}
|
||||
goto L16;
|
||||
|
||||
/* SEE IF RMS ERROR CRITERION IS SATISFIED (INPUT EPS .GT. 0.) */
|
||||
|
||||
L27:
|
||||
if (sigj <= etst) {
|
||||
goto L28;
|
||||
}
|
||||
if (*maxdeg == j) {
|
||||
goto L31;
|
||||
}
|
||||
goto L16;
|
||||
|
||||
/* RETURNS */
|
||||
|
||||
L28:
|
||||
*ierr = 1;
|
||||
*ndeg = j;
|
||||
sig = sigj;
|
||||
goto L33;
|
||||
L29:
|
||||
*ierr = 1;
|
||||
*ndeg = jpas;
|
||||
sig = sigpas;
|
||||
goto L33;
|
||||
L30:
|
||||
*ierr = 2;
|
||||
xermsg_("SLATEC", "DPOLFT", "INVALID INPUT PARAMETER.", &c__2, &c__1, (
|
||||
ftnlen)6, (ftnlen)6, (ftnlen)24);
|
||||
goto L37;
|
||||
L31:
|
||||
*ierr = 3;
|
||||
*ndeg = *maxdeg;
|
||||
sig = sigj;
|
||||
goto L33;
|
||||
L32:
|
||||
*ierr = 4;
|
||||
*ndeg = jpas;
|
||||
sig = sigpas;
|
||||
|
||||
L33:
|
||||
a[k3] = (doublereal) (*ndeg);
|
||||
|
||||
/* WHEN STATISTICAL TEST HAS BEEN USED, EVALUATE THE BEST POLYNOMIAL AT */
|
||||
/* ALL THE DATA POINTS IF R DOES NOT ALREADY CONTAIN THESE VALUES */
|
||||
|
||||
if (*eps >= 0.f || *ndeg == *maxdeg) {
|
||||
goto L36;
|
||||
}
|
||||
nder = 0;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
dp1vlu_(ndeg, &nder, &x[i__], &r__[i__], &yp, &a[1]);
|
||||
/* L35: */
|
||||
}
|
||||
L36:
|
||||
*eps = sqrt(sig / xm);
|
||||
L37:
|
||||
return 0;
|
||||
} /* dpolft_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
120
ext/f2c_math/fdump.c
Normal file
120
ext/f2c_math/fdump.c
Normal file
|
|
@ -0,0 +1,120 @@
|
|||
/* fdump.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* DECK FDUMP */
|
||||
/* Subroutine */ int fdump_(void)
|
||||
{
|
||||
/* ***BEGIN PROLOGUE FDUMP */
|
||||
/* ***PURPOSE Symbolic dump (should be locally written). */
|
||||
/* ***LIBRARY SLATEC (XERMSG) */
|
||||
/* ***CATEGORY R3 */
|
||||
/* ***TYPE ALL (FDUMP-A) */
|
||||
/* ***KEYWORDS ERROR, XERMSG */
|
||||
/* ***AUTHOR Jones, R. E., (SNLA) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* ***Note*** Machine Dependent Routine */
|
||||
/* FDUMP is intended to be replaced by a locally written */
|
||||
/* version which produces a symbolic dump. Failing this, */
|
||||
/* it should be replaced by a version which prints the */
|
||||
/* subprogram nesting list. Note that this dump must be */
|
||||
/* printed on each of up to five files, as indicated by the */
|
||||
/* XGETUA routine. See XSETUA and XGETUA for details. */
|
||||
|
||||
/* Written by Ron Jones, with SLATEC Common Math Library Subcommittee */
|
||||
|
||||
/* ***REFERENCES (NONE) */
|
||||
/* ***ROUTINES CALLED (NONE) */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 790801 DATE WRITTEN */
|
||||
/* 861211 REVISION DATE from Version 3.2 */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
|
||||
/* ***END PROLOGUE FDUMP */
|
||||
/* ***FIRST EXECUTABLE STATEMENT FDUMP */
|
||||
return 0;
|
||||
} /* fdump_ */
|
||||
|
||||
integer isamax_(integer *n, real *sx, integer *incx)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer ret_val, i__1;
|
||||
real r__1;
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, ix;
|
||||
static real smax;
|
||||
|
||||
|
||||
/* finds the index of element having max. absolute value. */
|
||||
/* jack dongarra, linpack, 3/11/78. */
|
||||
/* modified 3/93 to return if incx .le. 0. */
|
||||
|
||||
|
||||
/* Parameter adjustments */
|
||||
--sx;
|
||||
|
||||
/* Function Body */
|
||||
ret_val = 0;
|
||||
if (*n < 1 || *incx <= 0) {
|
||||
return ret_val;
|
||||
}
|
||||
ret_val = 1;
|
||||
if (*n == 1) {
|
||||
return ret_val;
|
||||
}
|
||||
if (*incx == 1) {
|
||||
goto L20;
|
||||
}
|
||||
|
||||
/* code for increment not equal to 1 */
|
||||
|
||||
ix = 1;
|
||||
smax = dabs(sx[1]);
|
||||
ix += *incx;
|
||||
i__1 = *n;
|
||||
for (i__ = 2; i__ <= i__1; ++i__) {
|
||||
if ((r__1 = sx[ix], dabs(r__1)) <= smax) {
|
||||
goto L5;
|
||||
}
|
||||
ret_val = i__;
|
||||
smax = (r__1 = sx[ix], dabs(r__1));
|
||||
L5:
|
||||
ix += *incx;
|
||||
/* L10: */
|
||||
}
|
||||
return ret_val;
|
||||
|
||||
/* code for increment equal to 1 */
|
||||
|
||||
L20:
|
||||
smax = dabs(sx[1]);
|
||||
i__1 = *n;
|
||||
for (i__ = 2; i__ <= i__1; ++i__) {
|
||||
if ((r__1 = sx[i__], dabs(r__1)) <= smax) {
|
||||
goto L30;
|
||||
}
|
||||
ret_val = i__;
|
||||
smax = (r__1 = sx[i__], dabs(r__1));
|
||||
L30:
|
||||
;
|
||||
}
|
||||
return ret_val;
|
||||
} /* isamax_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
144
ext/f2c_math/gmres.h
Normal file
144
ext/f2c_math/gmres.h
Normal file
|
|
@ -0,0 +1,144 @@
|
|||
// -*- C++ -*-
|
||||
|
||||
#ifndef GMRES_BLAS_H
|
||||
#define GMRES_BLAS_H
|
||||
|
||||
// ============================================================================
|
||||
//
|
||||
// GMRES nach Saad, Schultz
|
||||
// GMRES: a generalized minimal residual algorithm for solving nonsymmetric
|
||||
// linear systems
|
||||
// SIAM J Sci Stat Comput 7, 856-869 (1986)
|
||||
//
|
||||
// ----------------------------
|
||||
// Christian Badura, Mai 1998
|
||||
//
|
||||
// ============================================================================
|
||||
|
||||
|
||||
template< class Matrix >
|
||||
inline int
|
||||
gmres( int m, int N, const Matrix &A, const doublereal *b, doublereal *x, doublereal eps );
|
||||
|
||||
|
||||
template< class Matrix >
|
||||
inline int
|
||||
gmres( int m, int N, const Matrix &A, const doublereal *b, doublereal *x, doublereal eps,
|
||||
bool detailed );
|
||||
|
||||
|
||||
// ============================================================================
|
||||
|
||||
// #include "../../Cantera/src/blas.h"
|
||||
|
||||
#include "cblas.h"
|
||||
#include "../../Cantera/src/ctlapack.h"
|
||||
using namespace Cantera;
|
||||
|
||||
template< class Matrix >
|
||||
inline int
|
||||
gmres( int m, int n, const Matrix &A, const doublereal *b, doublereal *x, doublereal eps,
|
||||
bool detailed ) {
|
||||
if ( n<=0 )
|
||||
return -1;
|
||||
typedef doublereal *doublerealP;
|
||||
doublereal *V = new doublereal[n*(m+1)];
|
||||
doublereal *U = new doublereal[m*(m+1)/2];
|
||||
doublereal *r = new doublereal[n];
|
||||
doublereal *y = new doublereal[m+1];
|
||||
doublereal *c = new doublereal[m];
|
||||
doublereal *s = new doublereal[m];
|
||||
doublereal **v = new doublerealP[m+1];
|
||||
for ( int i=0; i<=m; ++i ) v[i]=V+i*n;
|
||||
int its=-1;
|
||||
{
|
||||
doublereal beta, h, rd, dd, nrm2b;
|
||||
int j, io, uij, u0j;
|
||||
nrm2b=dnrm2(n,b,1);
|
||||
cout << " norm = " << nrm2b << endl;
|
||||
io=0;
|
||||
do { // "aussere Iteration
|
||||
++io;
|
||||
//mult(A,x,r);
|
||||
A.mult(x,r);
|
||||
daxpy(n,-1.,b,1,r,1);
|
||||
beta=dnrm2(n,r,1);
|
||||
dcopy(n,r,1,v[0],1);
|
||||
dscal(n,1./beta,v[0],1);
|
||||
|
||||
y[0]=beta;
|
||||
j=0;
|
||||
uij=0;
|
||||
do { // innere Iteration j=0,...,m-1
|
||||
u0j=uij;
|
||||
//mult(A,v[j],v[j+1]);
|
||||
A.mult(v[j],v[j+1]);
|
||||
|
||||
ct_dgemv(ctlapack::ColMajor, ctlapack::Transpose, n, j+1, 1.0, V, n, v[j+1], 1, 0.0, U+u0j, 1);
|
||||
ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, n, j+1, -1.0, V, n, U+u0j, 1, 1.0, v[j+1], 1);
|
||||
|
||||
//dgemv(Transpose,n,j+1,1.,V,n,v[j+1],1,0.,U+u0j,1);
|
||||
//dgemv(NoTranspose,n,j+1,-1.,V,n,U+u0j,1,1.,v[j+1],1);
|
||||
|
||||
h=dnrm2(n,v[j+1],1);
|
||||
|
||||
dscal(n,1./h,v[j+1],1);
|
||||
|
||||
for ( int i=0; i<j; ++i ) { // rotiere neue Spalte
|
||||
doublereal tmp = c[i]*U[uij]-s[i]*U[uij+1];
|
||||
U[uij+1] = s[i]*U[uij]+c[i]*U[uij+1];
|
||||
U[uij] = tmp;
|
||||
++uij;
|
||||
}
|
||||
{ // berechne neue Rotation
|
||||
rd = U[uij];
|
||||
dd = sqrt(rd*rd+h*h);
|
||||
c[j] = rd/dd;
|
||||
s[j] = -h/dd;
|
||||
U[uij] = dd;
|
||||
++uij;
|
||||
}
|
||||
{ // rotiere rechte Seite y (vorher: y[j+1]=0)
|
||||
y[j+1] = s[j]*y[j];
|
||||
y[j] = c[j]*y[j];
|
||||
}
|
||||
++j;
|
||||
if ( detailed ) {
|
||||
cout<<"gmres("<<m<<")\t"<<io<<"\t"<<j<<"\t"<<y[j]<< endl;
|
||||
}
|
||||
} while ( j<m && fabs(y[j])>=eps*nrm2b );
|
||||
{ // minimiere bzgl Y
|
||||
dtpsv(UpperTriangle,NoTranspose,NotUnitTriangular,j,U,y,1);
|
||||
// korrigiere X
|
||||
dgemv(NoTranspose,n,j,-1.,V,n,y,1,1.,x,1);
|
||||
}
|
||||
} while ( fabs(y[j])>=eps*nrm2b );
|
||||
|
||||
// R"uckgabe: Zahl der inneren Iterationen
|
||||
its = m*(io-1)+j;
|
||||
}
|
||||
|
||||
delete[] V;
|
||||
delete[] U;
|
||||
delete[] r;
|
||||
delete[] y;
|
||||
delete[] c;
|
||||
delete[] s;
|
||||
delete[] v;
|
||||
return its;
|
||||
}
|
||||
|
||||
|
||||
// ============================================================================
|
||||
|
||||
|
||||
template< class Matrix >
|
||||
inline int
|
||||
gmres( int m, int n, const Matrix &A, const doublereal *b, doublereal *x, doublereal eps ){
|
||||
return gmres(m,n,A,b,x,eps,false);
|
||||
}
|
||||
|
||||
// ============================================================================
|
||||
|
||||
|
||||
#endif // GMRES_BLAS_H
|
||||
87
ext/f2c_math/idamax.c
Normal file
87
ext/f2c_math/idamax.c
Normal file
|
|
@ -0,0 +1,87 @@
|
|||
/* idamax.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
integer idamax_(integer *n, doublereal *dx, integer *incx)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer ret_val, i__1;
|
||||
doublereal d__1;
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, ix;
|
||||
static doublereal dmax__;
|
||||
|
||||
|
||||
/* finds the index of element having max. absolute value. */
|
||||
/* jack dongarra, linpack, 3/11/78. */
|
||||
/* modified 3/93 to return if incx .le. 0. */
|
||||
|
||||
|
||||
/* Parameter adjustments */
|
||||
--dx;
|
||||
|
||||
/* Function Body */
|
||||
ret_val = 0;
|
||||
if (*n < 1 || *incx <= 0) {
|
||||
return ret_val;
|
||||
}
|
||||
ret_val = 1;
|
||||
if (*n == 1) {
|
||||
return ret_val;
|
||||
}
|
||||
if (*incx == 1) {
|
||||
goto L20;
|
||||
}
|
||||
|
||||
/* code for increment not equal to 1 */
|
||||
|
||||
ix = 1;
|
||||
dmax__ = abs(dx[1]);
|
||||
ix += *incx;
|
||||
i__1 = *n;
|
||||
for (i__ = 2; i__ <= i__1; ++i__) {
|
||||
if ((d__1 = dx[ix], abs(d__1)) <= dmax__) {
|
||||
goto L5;
|
||||
}
|
||||
ret_val = i__;
|
||||
dmax__ = (d__1 = dx[ix], abs(d__1));
|
||||
L5:
|
||||
ix += *incx;
|
||||
/* L10: */
|
||||
}
|
||||
return ret_val;
|
||||
|
||||
/* code for increment equal to 1 */
|
||||
|
||||
L20:
|
||||
dmax__ = abs(dx[1]);
|
||||
i__1 = *n;
|
||||
for (i__ = 2; i__ <= i__1; ++i__) {
|
||||
if ((d__1 = dx[i__], abs(d__1)) <= dmax__) {
|
||||
goto L30;
|
||||
}
|
||||
ret_val = i__;
|
||||
dmax__ = (d__1 = dx[i__], abs(d__1));
|
||||
L30:
|
||||
;
|
||||
}
|
||||
return ret_val;
|
||||
} /* idamax_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
90
ext/f2c_math/j4save.c
Normal file
90
ext/f2c_math/j4save.c
Normal file
|
|
@ -0,0 +1,90 @@
|
|||
/* j4save.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* DECK J4SAVE */
|
||||
integer j4save_(integer *iwhich, integer *ivalue, logical *iset)
|
||||
{
|
||||
/* Initialized data */
|
||||
|
||||
static integer iparam[9] = { 0,2,0,10,1,0,0,0,0 };
|
||||
|
||||
/* System generated locals */
|
||||
integer ret_val;
|
||||
|
||||
/* ***BEGIN PROLOGUE J4SAVE */
|
||||
/* ***SUBSIDIARY */
|
||||
/* ***PURPOSE Save or recall global variables needed by error */
|
||||
/* handling routines. */
|
||||
/* ***LIBRARY SLATEC (XERROR) */
|
||||
/* ***TYPE INTEGER (J4SAVE-I) */
|
||||
/* ***KEYWORDS ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR */
|
||||
/* ***AUTHOR Jones, R. E., (SNLA) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* Abstract */
|
||||
/* J4SAVE saves and recalls several global variables needed */
|
||||
/* by the library error handling routines. */
|
||||
|
||||
/* Description of Parameters */
|
||||
/* --Input-- */
|
||||
/* IWHICH - Index of item desired. */
|
||||
/* = 1 Refers to current error number. */
|
||||
/* = 2 Refers to current error control flag. */
|
||||
/* = 3 Refers to current unit number to which error */
|
||||
/* messages are to be sent. (0 means use standard.) */
|
||||
/* = 4 Refers to the maximum number of times any */
|
||||
/* message is to be printed (as set by XERMAX). */
|
||||
/* = 5 Refers to the total number of units to which */
|
||||
/* each error message is to be written. */
|
||||
/* = 6 Refers to the 2nd unit for error messages */
|
||||
/* = 7 Refers to the 3rd unit for error messages */
|
||||
/* = 8 Refers to the 4th unit for error messages */
|
||||
/* = 9 Refers to the 5th unit for error messages */
|
||||
/* IVALUE - The value to be set for the IWHICH-th parameter, */
|
||||
/* if ISET is .TRUE. . */
|
||||
/* ISET - If ISET=.TRUE., the IWHICH-th parameter will BE */
|
||||
/* given the value, IVALUE. If ISET=.FALSE., the */
|
||||
/* IWHICH-th parameter will be unchanged, and IVALUE */
|
||||
/* is a dummy parameter. */
|
||||
/* --Output-- */
|
||||
/* The (old) value of the IWHICH-th parameter will be returned */
|
||||
/* in the function value, J4SAVE. */
|
||||
|
||||
/* ***SEE ALSO XERMSG */
|
||||
/* ***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC */
|
||||
/* Error-handling Package, SAND82-0800, Sandia */
|
||||
/* Laboratories, 1982. */
|
||||
/* ***ROUTINES CALLED (NONE) */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 790801 DATE WRITTEN */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
|
||||
/* 900205 Minor modifications to prologue. (WRB) */
|
||||
/* 900402 Added TYPE section. (WRB) */
|
||||
/* 910411 Added KEYWORDS section. (WRB) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* ***END PROLOGUE J4SAVE */
|
||||
/* ***FIRST EXECUTABLE STATEMENT J4SAVE */
|
||||
ret_val = iparam[(0 + (0 + (*iwhich - 1 << 2))) / 4];
|
||||
if (*iset) {
|
||||
iparam[*iwhich - 1] = *ivalue;
|
||||
}
|
||||
return ret_val;
|
||||
} /* j4save_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
60
ext/f2c_math/mach.cpp
Normal file
60
ext/f2c_math/mach.cpp
Normal file
|
|
@ -0,0 +1,60 @@
|
|||
|
||||
/* Standard C source for D1MACH -- remove the * in column 1 */
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
#include <float.h>
|
||||
#include <math.h>
|
||||
#include <limits.h>
|
||||
|
||||
extern "C" {
|
||||
|
||||
double d1mach_(long *i)
|
||||
{
|
||||
switch(*i){
|
||||
case 1: return DBL_MIN;
|
||||
case 2: return DBL_MAX;
|
||||
case 3: return DBL_EPSILON/FLT_RADIX;
|
||||
case 4: return DBL_EPSILON;
|
||||
case 5: return log10((double)FLT_RADIX);
|
||||
}
|
||||
fprintf(stderr, "invalid argument: d1mach(%ld)\n", *i);
|
||||
exit(1);
|
||||
return 0; /* some compilers demand return values */
|
||||
}
|
||||
|
||||
double d1mach(long *i) {return d1mach_(i);}
|
||||
|
||||
|
||||
long i1mach_(long *i)
|
||||
{
|
||||
switch(*i){
|
||||
case 1: return 5; /* standard input */
|
||||
case 2: return 6; /* standard output */
|
||||
case 3: return 7; /* standard punch */
|
||||
case 4: return 0; /* standard error */
|
||||
case 5: return 32; /* bits per integer */
|
||||
case 6: return sizeof(int);
|
||||
case 7: return 2; /* base for integers */
|
||||
case 8: return 31; /* digits of integer base */
|
||||
case 9: return LONG_MAX;
|
||||
case 10: return FLT_RADIX;
|
||||
case 11: return FLT_MANT_DIG;
|
||||
case 12: return FLT_MIN_EXP;
|
||||
case 13: return FLT_MAX_EXP;
|
||||
case 14: return DBL_MANT_DIG;
|
||||
case 15: return DBL_MIN_EXP;
|
||||
case 16: return DBL_MAX_EXP;
|
||||
}
|
||||
fprintf(stderr, "invalid argument: i1mach(%ld)\n", *i);
|
||||
exit(1);
|
||||
return 0; /* some compilers demand return values */
|
||||
}
|
||||
|
||||
long i1mach(long *i) { return i1mach_(i); }
|
||||
|
||||
long _i1mach_(long *i) {
|
||||
return i1mach_(i);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
644
ext/f2c_math/mkl_cblas.h
Normal file
644
ext/f2c_math/mkl_cblas.h
Normal file
|
|
@ -0,0 +1,644 @@
|
|||
/*
|
||||
// INTEL CORPORATION PROPRIETARY INFORMATION
|
||||
// This software is supplied under the terms of a license agreement or
|
||||
// nondisclosure agreement with Intel Corporation and may not be copied
|
||||
// or disclosed except in accordance with the terms of that agreement.
|
||||
// Copyright 1999, 2000 Intel Corporation. All Rights Reserved.
|
||||
//
|
||||
// File : mkl_cblas.h
|
||||
// Purpose : MKL CBLAS interface
|
||||
// Author : Shemyakin Andrey
|
||||
*/
|
||||
|
||||
#ifndef __MKL_CBLAS_H__
|
||||
#define __MKL_CBLAS_H__
|
||||
#include <stddef.h>
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" { /* Assume C declarations for C++ */
|
||||
#endif /* __cplusplus */
|
||||
|
||||
/*
|
||||
* Enumerated and derived types
|
||||
*/
|
||||
#define CBLAS_INDEX size_t /* this may vary between platforms */
|
||||
|
||||
typedef enum {CblasRowMajor=101, CblasColMajor=102} CBLAS_ORDER;
|
||||
typedef enum {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113} CBLAS_TRANSPOSE;
|
||||
typedef enum {CblasUpper=121, CblasLower=122} CBLAS_UPLO;
|
||||
typedef enum {CblasNonUnit=131, CblasUnit=132} CBLAS_DIAG;
|
||||
typedef enum {CblasLeft=141, CblasRight=142} CBLAS_SIDE;
|
||||
|
||||
/*
|
||||
* ===========================================================================
|
||||
* Prototypes for level 1 BLAS functions (complex are recast as routines)
|
||||
* ===========================================================================
|
||||
*/
|
||||
|
||||
float cblas_sdot(const int N, const float *X, const int incX,
|
||||
const float *Y, const int incY);
|
||||
float cblas_sdoti(const int N, const float *X, const int *indx,
|
||||
const float *Y);
|
||||
double cblas_ddot(const int N, const double *X, const int incX,
|
||||
const double *Y, const int incY);
|
||||
double cblas_ddoti(const int N, const double *X, const int *indx,
|
||||
const double *Y);
|
||||
|
||||
/*
|
||||
* Functions having prefixes Z and C only
|
||||
*/
|
||||
void cblas_cdotu_sub(const int N, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *dotu);
|
||||
void cblas_cdotui_sub(const int N, const void *X, const int *indx,
|
||||
const void *Y, void *dotui);
|
||||
void cblas_cdotc_sub(const int N, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *dotc);
|
||||
void cblas_cdotci_sub(const int N, const void *X, const int *indx,
|
||||
const void *Y, void *dotui);
|
||||
|
||||
void cblas_zdotu_sub(const int N, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *dotu);
|
||||
void cblas_zdotui_sub(const int N, const void *X, const int *indx,
|
||||
const void *Y, void *dotui);
|
||||
void cblas_zdotc_sub(const int N, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *dotc);
|
||||
void cblas_zdotci_sub(const int N, const void *X, const int *indx,
|
||||
const void *Y, void *dotui);
|
||||
|
||||
|
||||
/*
|
||||
* Functions having prefixes S D SC DZ
|
||||
*/
|
||||
float cblas_snrm2(const int N, const float *X, const int incX);
|
||||
float cblas_sasum(const int N, const float *X, const int incX);
|
||||
|
||||
double cblas_dnrm2(const int N, const double *X, const int incX);
|
||||
double cblas_dasum(const int N, const double *X, const int incX);
|
||||
|
||||
float cblas_scnrm2(const int N, const void *X, const int incX);
|
||||
float cblas_scasum(const int N, const void *X, const int incX);
|
||||
|
||||
double cblas_dznrm2(const int N, const void *X, const int incX);
|
||||
double cblas_dzasum(const int N, const void *X, const int incX);
|
||||
|
||||
|
||||
/*
|
||||
* Functions having standard 4 prefixes (S D C Z)
|
||||
*/
|
||||
CBLAS_INDEX cblas_isamax(const int N, const float *X, const int incX);
|
||||
CBLAS_INDEX cblas_idamax(const int N, const double *X, const int incX);
|
||||
CBLAS_INDEX cblas_icamax(const int N, const void *X, const int incX);
|
||||
CBLAS_INDEX cblas_izamax(const int N, const void *X, const int incX);
|
||||
CBLAS_INDEX cblas_isamin(const int N, const float *X, const int incX);
|
||||
CBLAS_INDEX cblas_idamin(const int N, const double *X, const int incX);
|
||||
CBLAS_INDEX cblas_icamin(const int N, const void *X, const int incX);
|
||||
CBLAS_INDEX cblas_izamin(const int N, const void *X, const int incX);
|
||||
|
||||
/*
|
||||
* ===========================================================================
|
||||
* Prototypes for level 1 BLAS routines
|
||||
* ===========================================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
* Routines with standard 4 prefixes (s, d, c, z)
|
||||
*/
|
||||
void cblas_sswap(const int N, float *X, const int incX,
|
||||
float *Y, const int incY);
|
||||
void cblas_scopy(const int N, const float *X, const int incX,
|
||||
float *Y, const int incY);
|
||||
void cblas_saxpy(const int N, const float alpha, const float *X,
|
||||
const int incX, float *Y, const int incY);
|
||||
void cblas_saxpyi(const int N, const float alpha, const float *X,
|
||||
const int *indx, float *Y);
|
||||
void cblas_sgthr(const int N, const float *Y, float *X,
|
||||
const int *indx);
|
||||
void cblas_sgthrz(const int N, float *Y, float *X,
|
||||
const int *indx);
|
||||
void cblas_ssctr(const int N, const float *X, const int *indx,
|
||||
float *Y);
|
||||
void cblas_srotg(float *a, float *b, float *c, float *s);
|
||||
|
||||
void cblas_dswap(const int N, double *X, const int incX,
|
||||
double *Y, const int incY);
|
||||
void cblas_dcopy(const int N, const double *X, const int incX,
|
||||
double *Y, const int incY);
|
||||
void cblas_daxpy(const int N, const double alpha, const double *X,
|
||||
const int incX, double *Y, const int incY);
|
||||
void cblas_daxpyi(const int N, const double alpha, const double *X,
|
||||
const int *indx, double *Y);
|
||||
void cblas_dgthr(const int N, const double *Y, double *X,
|
||||
const int *indx);
|
||||
void cblas_dgthrz(const int N, double *Y, double *X,
|
||||
const int *indx);
|
||||
void cblas_dsctr(const int N, const double *X, const int *indx,
|
||||
double *Y);
|
||||
void cblas_drotg(double *a, double *b, double *c, double *s);
|
||||
|
||||
void cblas_cswap(const int N, void *X, const int incX,
|
||||
void *Y, const int incY);
|
||||
void cblas_ccopy(const int N, const void *X, const int incX,
|
||||
void *Y, const int incY);
|
||||
void cblas_caxpy(const int N, const void *alpha, const void *X,
|
||||
const int incX, void *Y, const int incY);
|
||||
void cblas_caxpyi(const int N, const void *alpha, const void *X,
|
||||
const int *indx, void *Y);
|
||||
void cblas_cgthr(const int N, const void *Y, void *X,
|
||||
const int *indx);
|
||||
void cblas_cgthrz(const int N, void *Y, void *X,
|
||||
const int *indx);
|
||||
void cblas_csctr(const int N, const void *X, const int *indx,
|
||||
void *Y);
|
||||
void cblas_crotg(void *a, const void *b, float *c, void *s);
|
||||
|
||||
void cblas_zswap(const int N, void *X, const int incX,
|
||||
void *Y, const int incY);
|
||||
void cblas_zcopy(const int N, const void *X, const int incX,
|
||||
void *Y, const int incY);
|
||||
void cblas_zaxpy(const int N, const void *alpha, const void *X,
|
||||
const int incX, void *Y, const int incY);
|
||||
void cblas_zaxpyi(const int N, const void *alpha, const void *X,
|
||||
const int *indx, void *Y);
|
||||
void cblas_zgthr(const int N, const void *Y, void *X,
|
||||
const int *indx);
|
||||
void cblas_zgthrz(const int N, void *Y, void *X,
|
||||
const int *indx);
|
||||
void cblas_zsctr(const int N, const void *X, const int *indx,
|
||||
void *Y);
|
||||
void cblas_zrotg(void *a, const void *b, double *c, void *s);
|
||||
|
||||
/*
|
||||
* Routines with S and D prefix only
|
||||
*/
|
||||
void cblas_srotmg(float *d1, float *d2, float *b1, const float *b2, float *P);
|
||||
void cblas_srot(const int N, float *X, const int incX,
|
||||
float *Y, const int incY, const float c, const float s);
|
||||
void cblas_sroti(const int N, float *X, const int *indx,
|
||||
float *Y, const float c, const float s);
|
||||
void cblas_srotm(const int N, float *X, const int incX,
|
||||
float *Y, const int incY, const float *P);
|
||||
|
||||
void cblas_drotmg(double *d1, double *d2, double *b1, const double *b2, double *P);
|
||||
void cblas_drot(const int N, double *X, const int incX,
|
||||
double *Y, const int incY, const double c, const double s);
|
||||
void cblas_drotm(const int N, double *X, const int incX,
|
||||
double *Y, const int incY, const double *P);
|
||||
void cblas_droti(const int N, double *X, const int *indx,
|
||||
double *Y, const double c, const double s);
|
||||
|
||||
/*
|
||||
* Routines with CS and ZD prefix only
|
||||
*/
|
||||
void cblas_csrot(const int N, void *X, const int incX,
|
||||
void *Y, const int incY, const float c, const float s);
|
||||
void cblas_zdrot(const int N, void *X, const int incX,
|
||||
void *Y, const int incY, const double c, const double s);
|
||||
|
||||
/*
|
||||
* Routines with S D C Z CS and ZD prefixes
|
||||
*/
|
||||
void cblas_sscal(const int N, const float alpha, float *X, const int incX);
|
||||
void cblas_dscal(const int N, const double alpha, double *X, const int incX);
|
||||
void cblas_cscal(const int N, const void *alpha, void *X, const int incX);
|
||||
void cblas_zscal(const int N, const void *alpha, void *X, const int incX);
|
||||
void cblas_csscal(const int N, const float alpha, void *X, const int incX);
|
||||
void cblas_zdscal(const int N, const double alpha, void *X, const int incX);
|
||||
|
||||
/*
|
||||
* ===========================================================================
|
||||
* Prototypes for level 2 BLAS
|
||||
* ===========================================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
* Routines with standard 4 prefixes (S, D, C, Z)
|
||||
*/
|
||||
void cblas_sgemv(const CBLAS_ORDER order,
|
||||
const CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const float alpha, const float *A, const int lda,
|
||||
const float *X, const int incX, const float beta,
|
||||
float *Y, const int incY);
|
||||
void cblas_sgbmv(const CBLAS_ORDER order,
|
||||
const CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const int KL, const int KU, const float alpha,
|
||||
const float *A, const int lda, const float *X,
|
||||
const int incX, const float beta, float *Y, const int incY);
|
||||
void cblas_strmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const float *A, const int lda,
|
||||
float *X, const int incX);
|
||||
void cblas_stbmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const int K, const float *A, const int lda,
|
||||
float *X, const int incX);
|
||||
void cblas_stpmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const float *Ap, float *X, const int incX);
|
||||
void cblas_strsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const float *A, const int lda, float *X,
|
||||
const int incX);
|
||||
void cblas_stbsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const int K, const float *A, const int lda,
|
||||
float *X, const int incX);
|
||||
void cblas_stpsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const float *Ap, float *X, const int incX);
|
||||
|
||||
void cblas_dgemv(const CBLAS_ORDER order,
|
||||
const CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const double alpha, const double *A, const int lda,
|
||||
const double *X, const int incX, const double beta,
|
||||
double *Y, const int incY);
|
||||
void cblas_dgbmv(const CBLAS_ORDER order,
|
||||
const CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const int KL, const int KU, const double alpha,
|
||||
const double *A, const int lda, const double *X,
|
||||
const int incX, const double beta, double *Y, const int incY);
|
||||
void cblas_dtrmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const double *A, const int lda,
|
||||
double *X, const int incX);
|
||||
void cblas_dtbmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const int K, const double *A, const int lda,
|
||||
double *X, const int incX);
|
||||
void cblas_dtpmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const double *Ap, double *X, const int incX);
|
||||
void cblas_dtrsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const double *A, const int lda, double *X,
|
||||
const int incX);
|
||||
void cblas_dtbsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const int K, const double *A, const int lda,
|
||||
double *X, const int incX);
|
||||
void cblas_dtpsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const double *Ap, double *X, const int incX);
|
||||
|
||||
void cblas_cgemv(const CBLAS_ORDER order,
|
||||
const CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *X, const int incX, const void *beta,
|
||||
void *Y, const int incY);
|
||||
void cblas_cgbmv(const CBLAS_ORDER order,
|
||||
const CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const int KL, const int KU, const void *alpha,
|
||||
const void *A, const int lda, const void *X,
|
||||
const int incX, const void *beta, void *Y, const int incY);
|
||||
void cblas_ctrmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ctbmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const int K, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ctpmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const void *Ap, void *X, const int incX);
|
||||
void cblas_ctrsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const void *A, const int lda, void *X,
|
||||
const int incX);
|
||||
void cblas_ctbsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const int K, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ctpsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const void *Ap, void *X, const int incX);
|
||||
|
||||
void cblas_zgemv(const CBLAS_ORDER order,
|
||||
const CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *X, const int incX, const void *beta,
|
||||
void *Y, const int incY);
|
||||
void cblas_zgbmv(const CBLAS_ORDER order,
|
||||
const CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const int KL, const int KU, const void *alpha,
|
||||
const void *A, const int lda, const void *X,
|
||||
const int incX, const void *beta, void *Y, const int incY);
|
||||
void cblas_ztrmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ztbmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const int K, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ztpmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const void *Ap, void *X, const int incX);
|
||||
void cblas_ztrsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const void *A, const int lda, void *X,
|
||||
const int incX);
|
||||
void cblas_ztbsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const int K, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ztpsv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE TransA, const CBLAS_DIAG Diag,
|
||||
const int N, const void *Ap, void *X, const int incX);
|
||||
|
||||
|
||||
/*
|
||||
* Routines with S and D prefixes only
|
||||
*/
|
||||
void cblas_ssymv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *A,
|
||||
const int lda, const float *X, const int incX,
|
||||
const float beta, float *Y, const int incY);
|
||||
void cblas_ssbmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const int K, const float alpha, const float *A,
|
||||
const int lda, const float *X, const int incX,
|
||||
const float beta, float *Y, const int incY);
|
||||
void cblas_sspmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *Ap,
|
||||
const float *X, const int incX,
|
||||
const float beta, float *Y, const int incY);
|
||||
void cblas_sger(const CBLAS_ORDER order, const int M, const int N,
|
||||
const float alpha, const float *X, const int incX,
|
||||
const float *Y, const int incY, float *A, const int lda);
|
||||
void cblas_ssyr(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *X,
|
||||
const int incX, float *A, const int lda);
|
||||
void cblas_sspr(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *X,
|
||||
const int incX, float *Ap);
|
||||
void cblas_ssyr2(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *X,
|
||||
const int incX, const float *Y, const int incY, float *A,
|
||||
const int lda);
|
||||
void cblas_sspr2(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *X,
|
||||
const int incX, const float *Y, const int incY, float *A);
|
||||
|
||||
void cblas_dsymv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *A,
|
||||
const int lda, const double *X, const int incX,
|
||||
const double beta, double *Y, const int incY);
|
||||
void cblas_dsbmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const int K, const double alpha, const double *A,
|
||||
const int lda, const double *X, const int incX,
|
||||
const double beta, double *Y, const int incY);
|
||||
void cblas_dspmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *Ap,
|
||||
const double *X, const int incX,
|
||||
const double beta, double *Y, const int incY);
|
||||
void cblas_dger(const CBLAS_ORDER order, const int M, const int N,
|
||||
const double alpha, const double *X, const int incX,
|
||||
const double *Y, const int incY, double *A, const int lda);
|
||||
void cblas_dsyr(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *X,
|
||||
const int incX, double *A, const int lda);
|
||||
void cblas_dspr(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *X,
|
||||
const int incX, double *Ap);
|
||||
void cblas_dsyr2(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *X,
|
||||
const int incX, const double *Y, const int incY, double *A,
|
||||
const int lda);
|
||||
void cblas_dspr2(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *X,
|
||||
const int incX, const double *Y, const int incY, double *A);
|
||||
|
||||
|
||||
/*
|
||||
* Routines with C and Z prefixes only
|
||||
*/
|
||||
void cblas_chemv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const void *alpha, const void *A,
|
||||
const int lda, const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_chbmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const int K, const void *alpha, const void *A,
|
||||
const int lda, const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_chpmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const void *alpha, const void *Ap,
|
||||
const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_cgeru(const CBLAS_ORDER order, const int M, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_cgerc(const CBLAS_ORDER order, const int M, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_cher(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const void *X, const int incX,
|
||||
void *A, const int lda);
|
||||
void cblas_chpr(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const void *X,
|
||||
const int incX, void *A);
|
||||
void cblas_cher2(const CBLAS_ORDER order, const CBLAS_UPLO Uplo, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_chpr2(const CBLAS_ORDER order, const CBLAS_UPLO Uplo, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *Ap);
|
||||
|
||||
void cblas_zhemv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const void *alpha, const void *A,
|
||||
const int lda, const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_zhbmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const int K, const void *alpha, const void *A,
|
||||
const int lda, const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_zhpmv(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const void *alpha, const void *Ap,
|
||||
const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_zgeru(const CBLAS_ORDER order, const int M, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_zgerc(const CBLAS_ORDER order, const int M, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_zher(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const void *X, const int incX,
|
||||
void *A, const int lda);
|
||||
void cblas_zhpr(const CBLAS_ORDER order, const CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const void *X,
|
||||
const int incX, void *A);
|
||||
void cblas_zher2(const CBLAS_ORDER order, const CBLAS_UPLO Uplo, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_zhpr2(const CBLAS_ORDER order, const CBLAS_UPLO Uplo, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *Ap);
|
||||
|
||||
/*
|
||||
* ===========================================================================
|
||||
* Prototypes for level 3 BLAS
|
||||
* ===========================================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
* Routines with standard 4 prefixes (S, D, C, Z)
|
||||
*/
|
||||
void cblas_sgemm(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_TRANSPOSE TransB, const int M, const int N,
|
||||
const int K, const float alpha, const float *A,
|
||||
const int lda, const float *B, const int ldb,
|
||||
const float beta, float *C, const int ldc);
|
||||
void cblas_ssymm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const float alpha, const float *A, const int lda,
|
||||
const float *B, const int ldb, const float beta,
|
||||
float *C, const int ldc);
|
||||
void cblas_ssyrk(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const float alpha, const float *A, const int lda,
|
||||
const float beta, float *C, const int ldc);
|
||||
void cblas_ssyr2k(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const float alpha, const float *A, const int lda,
|
||||
const float *B, const int ldb, const float beta,
|
||||
float *C, const int ldc);
|
||||
void cblas_strmm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_DIAG Diag, const int M, const int N,
|
||||
const float alpha, const float *A, const int lda,
|
||||
float *B, const int ldb);
|
||||
void cblas_strsm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_DIAG Diag, const int M, const int N,
|
||||
const float alpha, const float *A, const int lda,
|
||||
float *B, const int ldb);
|
||||
|
||||
void cblas_dgemm(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_TRANSPOSE TransB, const int M, const int N,
|
||||
const int K, const double alpha, const double *A,
|
||||
const int lda, const double *B, const int ldb,
|
||||
const double beta, double *C, const int ldc);
|
||||
void cblas_dsymm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const double alpha, const double *A, const int lda,
|
||||
const double *B, const int ldb, const double beta,
|
||||
double *C, const int ldc);
|
||||
void cblas_dsyrk(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const double alpha, const double *A, const int lda,
|
||||
const double beta, double *C, const int ldc);
|
||||
void cblas_dsyr2k(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const double alpha, const double *A, const int lda,
|
||||
const double *B, const int ldb, const double beta,
|
||||
double *C, const int ldc);
|
||||
void cblas_dtrmm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_DIAG Diag, const int M, const int N,
|
||||
const double alpha, const double *A, const int lda,
|
||||
double *B, const int ldb);
|
||||
void cblas_dtrsm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_DIAG Diag, const int M, const int N,
|
||||
const double alpha, const double *A, const int lda,
|
||||
double *B, const int ldb);
|
||||
|
||||
void cblas_cgemm(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_TRANSPOSE TransB, const int M, const int N,
|
||||
const int K, const void *alpha, const void *A,
|
||||
const int lda, const void *B, const int ldb,
|
||||
const void *beta, void *C, const int ldc);
|
||||
void cblas_csymm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_csyrk(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *beta, void *C, const int ldc);
|
||||
void cblas_csyr2k(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_ctrmm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_DIAG Diag, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
void *B, const int ldb);
|
||||
void cblas_ctrsm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_DIAG Diag, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
void *B, const int ldb);
|
||||
|
||||
void cblas_zgemm(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_TRANSPOSE TransB, const int M, const int N,
|
||||
const int K, const void *alpha, const void *A,
|
||||
const int lda, const void *B, const int ldb,
|
||||
const void *beta, void *C, const int ldc);
|
||||
void cblas_zsymm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_zsyrk(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *beta, void *C, const int ldc);
|
||||
void cblas_zsyr2k(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_ztrmm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_DIAG Diag, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
void *B, const int ldb);
|
||||
void cblas_ztrsm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const CBLAS_TRANSPOSE TransA,
|
||||
const CBLAS_DIAG Diag, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
void *B, const int ldb);
|
||||
|
||||
|
||||
/*
|
||||
* Routines with prefixes C and Z only
|
||||
*/
|
||||
void cblas_chemm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_cherk(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const float alpha, const void *A, const int lda,
|
||||
const float beta, void *C, const int ldc);
|
||||
void cblas_cher2k(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const float beta,
|
||||
void *C, const int ldc);
|
||||
|
||||
void cblas_zhemm(const CBLAS_ORDER Order, const CBLAS_SIDE Side,
|
||||
const CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_zherk(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const double alpha, const void *A, const int lda,
|
||||
const double beta, void *C, const int ldc);
|
||||
void cblas_zher2k(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
|
||||
const CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const double beta,
|
||||
void *C, const int ldc);
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif /* __cplusplus */
|
||||
|
||||
#endif /* __MKL_CBLAS_H__ */
|
||||
1006
ext/f2c_math/pcoef.c
Normal file
1006
ext/f2c_math/pcoef.c
Normal file
File diff suppressed because it is too large
Load diff
533
ext/f2c_math/polfit.c
Normal file
533
ext/f2c_math/polfit.c
Normal file
|
|
@ -0,0 +1,533 @@
|
|||
/* polfit.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__2 = 2;
|
||||
static integer c__1 = 1;
|
||||
|
||||
/* DECK POLFIT */
|
||||
/* Subroutine */ int polfit_(integer *n, real *x, real *y, real *w, integer *
|
||||
maxdeg, integer *ndeg, real *eps, real *r__, integer *ierr, real *a)
|
||||
{
|
||||
/* Initialized data */
|
||||
|
||||
static real co[12] /* was [4][3] */ = { -13.08685f,-2.4648165f,
|
||||
-3.3846535f,-1.2973162f,-3.3381146f,-1.7812271f,-3.2578406f,
|
||||
-1.6589279f,-1.6282703f,-1.3152745f,-3.2640179f,-1.9829776f };
|
||||
|
||||
/* System generated locals */
|
||||
integer i__1;
|
||||
real r__1;
|
||||
|
||||
/* Builtin functions */
|
||||
double sqrt(doublereal);
|
||||
|
||||
/* Local variables */
|
||||
static real f;
|
||||
static integer i__, j, m, k1, k2, k3, k4, k5;
|
||||
static real w1, w11, xm, yp;
|
||||
static integer jp1;
|
||||
static real den, sig;
|
||||
static integer k1pj, k2pj, k4pi, k5pi, k3pi, mop1;
|
||||
static real degf;
|
||||
static integer nder;
|
||||
static real sigj;
|
||||
static integer ksig, jpas;
|
||||
static real temp, etst;
|
||||
static doublereal temd1, temd2;
|
||||
static integer idegf, nfail;
|
||||
static real fcrit, sigjm1, sigpas;
|
||||
extern /* Subroutine */ int pvalue_(integer *, integer *, real *, real *,
|
||||
real *, real *), xermsg_(char *, char *, char *, integer *,
|
||||
integer *, ftnlen, ftnlen, ftnlen);
|
||||
|
||||
/* ***BEGIN PROLOGUE POLFIT */
|
||||
/* ***PURPOSE Fit discrete data in a least squares sense by polynomials */
|
||||
/* in one variable. */
|
||||
/* ***LIBRARY SLATEC */
|
||||
/* ***CATEGORY K1A1A2 */
|
||||
/* ***TYPE SINGLE PRECISION (POLFIT-S, DPOLFT-D) */
|
||||
/* ***KEYWORDS CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT */
|
||||
/* ***AUTHOR Shampine, L. F., (SNLA) */
|
||||
/* Davenport, S. M., (SNLA) */
|
||||
/* Huddleston, R. E., (SNLL) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* Abstract */
|
||||
|
||||
/* 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 POLFIT 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 PVALUE may then be */
|
||||
/* called to evaluate the fitted polynomials and any of their */
|
||||
/* derivatives at any point. The subroutine PCOEF may be used to */
|
||||
/* express the polynomial fits as powers of (X-C) for any specified */
|
||||
/* point C. */
|
||||
|
||||
/* The parameters for POLFIT are */
|
||||
|
||||
/* Input -- */
|
||||
/* N - the number of data points. The arrays X, Y and W */
|
||||
/* must be dimensioned at least N (N .GE. 1). */
|
||||
/* X - array of values of the independent variable. These */
|
||||
/* values may appear in any order and need not all be */
|
||||
/* distinct. */
|
||||
/* Y - array of corresponding function values. */
|
||||
/* W - array of positive values to be used as weights. If */
|
||||
/* W(1) is negative, POLFIT 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 . */
|
||||
/* 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. */
|
||||
/* EPS - specifies the criterion to be used in determining */
|
||||
/* the degree of fit to be computed. */
|
||||
/* (1) If EPS is input negative, POLFIT 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., POLFIT 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. POLFIT will increase the */
|
||||
/* degree of fit until this criterion is met or */
|
||||
/* until the maximum degree is reached. */
|
||||
|
||||
/* Output -- */
|
||||
/* NDEG - degree of the highest degree fit computed. */
|
||||
/* EPS - RMS error of the polynomial of degree NDEG . */
|
||||
/* R - vector of dimension at least NDEG containing values */
|
||||
/* of the fit of degree NDEG at each of the X(I) . */
|
||||
/* Except when the statistical test is used, these */
|
||||
/* values are more accurate than results from subroutine */
|
||||
/* PVALUE normally are. */
|
||||
/* IERR - error flag with the following possible values. */
|
||||
/* 1 -- indicates normal execution, i.e., either */
|
||||
/* (1) the input value of EPS was negative, and the */
|
||||
/* computed polynomial fit of degree NDEG */
|
||||
/* satisfies the specified F test, or */
|
||||
/* (2) the input value of EPS was 0., and the fits of */
|
||||
/* all degrees up to MAXDEG are complete, or */
|
||||
/* (3) the input value of EPS was positive, and the */
|
||||
/* polynomial of degree NDEG satisfies the RMS */
|
||||
/* error requirement. */
|
||||
/* 2 -- invalid input parameter. At least one of the input */
|
||||
/* parameters has an illegal value and must be corrected */
|
||||
/* before POLFIT can proceed. Valid input results */
|
||||
/* when the following restrictions are observed */
|
||||
/* N .GE. 1 */
|
||||
/* 0 .LE. MAXDEG .LE. N-1 for EPS .GE. 0. */
|
||||
/* 0 .LE. MAXDEG .LE. N-2 for EPS .LT. 0. */
|
||||
/* W(1)=-1.0 or W(I) .GT. 0., I=1,...,N . */
|
||||
/* 3 -- cannot satisfy the RMS error requirement with a */
|
||||
/* polynomial of degree no greater than MAXDEG . Best */
|
||||
/* fit found is of degree MAXDEG . */
|
||||
/* 4 -- cannot satisfy the test for significance using */
|
||||
/* current value of MAXDEG . Statistically, the */
|
||||
/* best fit found is of order NORD . (In this case, */
|
||||
/* NDEG will have one of the values: MAXDEG-2, */
|
||||
/* MAXDEG-1, or MAXDEG). Using a higher value of */
|
||||
/* MAXDEG may result in passing the test. */
|
||||
/* A - work and output array having at least 3N+3MAXDEG+3 */
|
||||
/* locations */
|
||||
|
||||
/* Note - POLFIT calculates all fits of degrees up to and including */
|
||||
/* NDEG . Any or all of these fits can be evaluated or */
|
||||
/* expressed as powers of (X-C) using PVALUE and PCOEF */
|
||||
/* after just one call to POLFIT . */
|
||||
|
||||
/* ***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, */
|
||||
/* Curve fitting by polynomials in one variable, Report */
|
||||
/* SLA-74-0270, Sandia Laboratories, June 1974. */
|
||||
/* ***ROUTINES CALLED PVALUE, XERMSG */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 740601 DATE WRITTEN */
|
||||
/* 890531 Changed all specific intrinsics to generic. (WRB) */
|
||||
/* 890531 REVISION DATE from Version 3.2 */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
|
||||
/* 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* 920527 Corrected erroneous statements in DESCRIPTION. (WRB) */
|
||||
/* ***END PROLOGUE POLFIT */
|
||||
/* Parameter adjustments */
|
||||
--a;
|
||||
--r__;
|
||||
--w;
|
||||
--y;
|
||||
--x;
|
||||
|
||||
/* Function Body */
|
||||
/* ***FIRST EXECUTABLE STATEMENT POLFIT */
|
||||
m = abs(*n);
|
||||
if (m == 0) {
|
||||
goto L30;
|
||||
}
|
||||
if (*maxdeg < 0) {
|
||||
goto L30;
|
||||
}
|
||||
a[1] = (real) (*maxdeg);
|
||||
mop1 = *maxdeg + 1;
|
||||
if (m < mop1) {
|
||||
goto L30;
|
||||
}
|
||||
if (*eps < 0.f && m == mop1) {
|
||||
goto L30;
|
||||
}
|
||||
xm = (real) m;
|
||||
etst = *eps * *eps * xm;
|
||||
if (w[1] < 0.f) {
|
||||
goto L2;
|
||||
}
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
if (w[i__] <= 0.f) {
|
||||
goto L30;
|
||||
}
|
||||
/* L1: */
|
||||
}
|
||||
goto L4;
|
||||
L2:
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
/* L3: */
|
||||
w[i__] = 1.f;
|
||||
}
|
||||
L4:
|
||||
if (*eps >= 0.f) {
|
||||
goto L8;
|
||||
}
|
||||
|
||||
/* DETERMINE SIGNIFICANCE LEVEL INDEX TO BE USED IN STATISTICAL TEST FOR */
|
||||
/* CHOOSING DEGREE OF POLYNOMIAL FIT */
|
||||
|
||||
if (*eps > -.55f) {
|
||||
goto L5;
|
||||
}
|
||||
idegf = m - *maxdeg - 1;
|
||||
ksig = 1;
|
||||
if (idegf < 10) {
|
||||
ksig = 2;
|
||||
}
|
||||
if (idegf < 5) {
|
||||
ksig = 3;
|
||||
}
|
||||
goto L8;
|
||||
L5:
|
||||
ksig = 1;
|
||||
if (*eps < -.03f) {
|
||||
ksig = 2;
|
||||
}
|
||||
if (*eps < -.07f) {
|
||||
ksig = 3;
|
||||
}
|
||||
|
||||
/* INITIALIZE INDEXES AND COEFFICIENTS FOR FITTING */
|
||||
|
||||
L8:
|
||||
k1 = *maxdeg + 1;
|
||||
k2 = k1 + *maxdeg;
|
||||
k3 = k2 + *maxdeg + 2;
|
||||
k4 = k3 + m;
|
||||
k5 = k4 + m;
|
||||
i__1 = k4;
|
||||
for (i__ = 2; i__ <= i__1; ++i__) {
|
||||
/* L9: */
|
||||
a[i__] = 0.f;
|
||||
}
|
||||
w11 = 0.f;
|
||||
if (*n < 0) {
|
||||
goto L11;
|
||||
}
|
||||
|
||||
/* UNCONSTRAINED CASE */
|
||||
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
a[k4pi] = 1.f;
|
||||
/* L10: */
|
||||
w11 += w[i__];
|
||||
}
|
||||
goto L13;
|
||||
|
||||
/* CONSTRAINED CASE */
|
||||
|
||||
L11:
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
/* L12: */
|
||||
/* Computing 2nd power */
|
||||
r__1 = a[k4pi];
|
||||
w11 += w[i__] * (r__1 * r__1);
|
||||
}
|
||||
|
||||
/* COMPUTE FIT OF DEGREE ZERO */
|
||||
|
||||
L13:
|
||||
temd1 = 0.;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
temd1 += (doublereal) w[i__] * (doublereal) y[i__] * (doublereal) a[
|
||||
k4pi];
|
||||
/* L14: */
|
||||
}
|
||||
temd1 /= (doublereal) w11;
|
||||
a[k2 + 1] = temd1;
|
||||
sigj = 0.f;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
k5pi = k5 + i__;
|
||||
temd2 = temd1 * (doublereal) a[k4pi];
|
||||
r__[i__] = temd2;
|
||||
a[k5pi] = temd2 - (doublereal) r__[i__];
|
||||
/* L15: */
|
||||
/* Computing 2nd power */
|
||||
r__1 = y[i__] - r__[i__] - a[k5pi];
|
||||
sigj += w[i__] * (r__1 * r__1);
|
||||
}
|
||||
j = 0;
|
||||
|
||||
/* SEE IF POLYNOMIAL OF DEGREE 0 SATISFIES THE DEGREE SELECTION CRITERION */
|
||||
|
||||
if (*eps < 0.f) {
|
||||
goto L24;
|
||||
} else if (*eps == 0) {
|
||||
goto L26;
|
||||
} else {
|
||||
goto L27;
|
||||
}
|
||||
|
||||
/* INCREMENT DEGREE */
|
||||
|
||||
L16:
|
||||
++j;
|
||||
jp1 = j + 1;
|
||||
k1pj = k1 + j;
|
||||
k2pj = k2 + j;
|
||||
sigjm1 = sigj;
|
||||
|
||||
/* COMPUTE NEW B COEFFICIENT EXCEPT WHEN J = 1 */
|
||||
|
||||
if (j > 1) {
|
||||
a[k1pj] = w11 / w1;
|
||||
}
|
||||
|
||||
/* COMPUTE NEW A COEFFICIENT */
|
||||
|
||||
temd1 = 0.;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
temd2 = a[k4pi];
|
||||
temd1 += (doublereal) x[i__] * (doublereal) w[i__] * temd2 * temd2;
|
||||
/* L18: */
|
||||
}
|
||||
a[jp1] = temd1 / (doublereal) w11;
|
||||
|
||||
/* EVALUATE ORTHOGONAL POLYNOMIAL AT DATA POINTS */
|
||||
|
||||
w1 = w11;
|
||||
w11 = 0.f;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k3pi = k3 + i__;
|
||||
k4pi = k4 + i__;
|
||||
temp = a[k3pi];
|
||||
a[k3pi] = a[k4pi];
|
||||
a[k4pi] = (x[i__] - a[jp1]) * a[k3pi] - a[k1pj] * temp;
|
||||
/* L19: */
|
||||
/* Computing 2nd power */
|
||||
r__1 = a[k4pi];
|
||||
w11 += w[i__] * (r__1 * r__1);
|
||||
}
|
||||
|
||||
/* GET NEW ORTHOGONAL POLYNOMIAL COEFFICIENT USING PARTIAL DOUBLE */
|
||||
/* PRECISION */
|
||||
|
||||
temd1 = 0.;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
k5pi = k5 + i__;
|
||||
temd2 = (doublereal) w[i__] * (doublereal) (y[i__] - r__[i__] - a[
|
||||
k5pi]) * (doublereal) a[k4pi];
|
||||
/* L20: */
|
||||
temd1 += temd2;
|
||||
}
|
||||
temd1 /= (doublereal) w11;
|
||||
a[k2pj + 1] = temd1;
|
||||
|
||||
/* UPDATE POLYNOMIAL EVALUATIONS AT EACH OF THE DATA POINTS, AND */
|
||||
/* ACCUMULATE SUM OF SQUARES OF ERRORS. THE POLYNOMIAL EVALUATIONS ARE */
|
||||
/* COMPUTED AND STORED IN EXTENDED PRECISION. FOR THE I-TH DATA POINT, */
|
||||
/* THE MOST SIGNIFICANT BITS ARE STORED IN R(I) , AND THE LEAST */
|
||||
/* SIGNIFICANT BITS ARE IN A(K5PI) . */
|
||||
|
||||
sigj = 0.f;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
k4pi = k4 + i__;
|
||||
k5pi = k5 + i__;
|
||||
temd2 = (doublereal) r__[i__] + (doublereal) a[k5pi] + temd1 * (
|
||||
doublereal) a[k4pi];
|
||||
r__[i__] = temd2;
|
||||
a[k5pi] = temd2 - (doublereal) r__[i__];
|
||||
/* L21: */
|
||||
/* Computing 2nd power */
|
||||
r__1 = y[i__] - r__[i__] - a[k5pi];
|
||||
sigj += w[i__] * (r__1 * r__1);
|
||||
}
|
||||
|
||||
/* SEE IF DEGREE SELECTION CRITERION HAS BEEN SATISFIED OR IF DEGREE */
|
||||
/* MAXDEG HAS BEEN REACHED */
|
||||
|
||||
if (*eps < 0.f) {
|
||||
goto L23;
|
||||
} else if (*eps == 0) {
|
||||
goto L26;
|
||||
} else {
|
||||
goto L27;
|
||||
}
|
||||
|
||||
/* COMPUTE F STATISTICS (INPUT EPS .LT. 0.) */
|
||||
|
||||
L23:
|
||||
if (sigj == 0.f) {
|
||||
goto L29;
|
||||
}
|
||||
degf = (real) (m - j - 1);
|
||||
den = (co[(ksig << 2) - 1] * degf + 1.f) * degf;
|
||||
fcrit = ((co[(ksig << 2) - 2] * degf + co[(ksig << 2) - 3]) * degf + co[(
|
||||
ksig << 2) - 4]) / den;
|
||||
fcrit *= fcrit;
|
||||
f = (sigjm1 - sigj) * degf / sigj;
|
||||
if (f < fcrit) {
|
||||
goto L25;
|
||||
}
|
||||
|
||||
/* POLYNOMIAL OF DEGREE J SATISFIES F TEST */
|
||||
|
||||
L24:
|
||||
sigpas = sigj;
|
||||
jpas = j;
|
||||
nfail = 0;
|
||||
if (*maxdeg == j) {
|
||||
goto L32;
|
||||
}
|
||||
goto L16;
|
||||
|
||||
/* POLYNOMIAL OF DEGREE J FAILS F TEST. IF THERE HAVE BEEN THREE */
|
||||
/* SUCCESSIVE FAILURES, A STATISTICALLY BEST DEGREE HAS BEEN FOUND. */
|
||||
|
||||
L25:
|
||||
++nfail;
|
||||
if (nfail >= 3) {
|
||||
goto L29;
|
||||
}
|
||||
if (*maxdeg == j) {
|
||||
goto L32;
|
||||
}
|
||||
goto L16;
|
||||
|
||||
/* RAISE THE DEGREE IF DEGREE MAXDEG HAS NOT YET BEEN REACHED (INPUT */
|
||||
/* EPS = 0.) */
|
||||
|
||||
L26:
|
||||
if (*maxdeg == j) {
|
||||
goto L28;
|
||||
}
|
||||
goto L16;
|
||||
|
||||
/* SEE IF RMS ERROR CRITERION IS SATISFIED (INPUT EPS .GT. 0.) */
|
||||
|
||||
L27:
|
||||
if (sigj <= etst) {
|
||||
goto L28;
|
||||
}
|
||||
if (*maxdeg == j) {
|
||||
goto L31;
|
||||
}
|
||||
goto L16;
|
||||
|
||||
/* RETURNS */
|
||||
|
||||
L28:
|
||||
*ierr = 1;
|
||||
*ndeg = j;
|
||||
sig = sigj;
|
||||
goto L33;
|
||||
L29:
|
||||
*ierr = 1;
|
||||
*ndeg = jpas;
|
||||
sig = sigpas;
|
||||
goto L33;
|
||||
L30:
|
||||
*ierr = 2;
|
||||
xermsg_("SLATEC", "POLFIT", "INVALID INPUT PARAMETER.", &c__2, &c__1, (
|
||||
ftnlen)6, (ftnlen)6, (ftnlen)24);
|
||||
goto L37;
|
||||
L31:
|
||||
*ierr = 3;
|
||||
*ndeg = *maxdeg;
|
||||
sig = sigj;
|
||||
goto L33;
|
||||
L32:
|
||||
*ierr = 4;
|
||||
*ndeg = jpas;
|
||||
sig = sigpas;
|
||||
|
||||
L33:
|
||||
a[k3] = (real) (*ndeg);
|
||||
|
||||
/* WHEN STATISTICAL TEST HAS BEEN USED, EVALUATE THE BEST POLYNOMIAL AT */
|
||||
/* ALL THE DATA POINTS IF R DOES NOT ALREADY CONTAIN THESE VALUES */
|
||||
|
||||
if (*eps >= 0.f || *ndeg == *maxdeg) {
|
||||
goto L36;
|
||||
}
|
||||
nder = 0;
|
||||
i__1 = m;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
pvalue_(ndeg, &nder, &x[i__], &r__[i__], &yp, &a[1]);
|
||||
/* L35: */
|
||||
}
|
||||
L36:
|
||||
*eps = sqrt(sig / xm);
|
||||
L37:
|
||||
return 0;
|
||||
} /* polfit_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
255
ext/f2c_math/pvalue.c
Normal file
255
ext/f2c_math/pvalue.c
Normal file
|
|
@ -0,0 +1,255 @@
|
|||
/* pvalue.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__1 = 1;
|
||||
static integer c__8 = 8;
|
||||
static integer c__2 = 2;
|
||||
static integer c__5 = 5;
|
||||
|
||||
/* DECK PVALUE */
|
||||
/* Subroutine */ int pvalue_(integer *l, integer *nder, real *x, real *yfit,
|
||||
real *yp, real *a)
|
||||
{
|
||||
/* System generated locals */
|
||||
address a__1[5];
|
||||
integer i__1, i__2, i__3[5];
|
||||
char ch__1[150];
|
||||
|
||||
/* Builtin functions */
|
||||
integer s_wsfi(icilist *), do_fio(integer *, char *, ftnlen), e_wsfi(void)
|
||||
;
|
||||
/* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, n, k1, k2, k3, k4;
|
||||
static real cc;
|
||||
static integer ic, kc, in, k1i, lm1, lp1;
|
||||
static real dif;
|
||||
static integer k3p1, k4p1, ndo;
|
||||
static real val;
|
||||
static integer ilo, iup, ndp1, inp1, k3pn, k4pn, nord;
|
||||
static char xern1[8], xern2[8];
|
||||
static integer maxord;
|
||||
extern /* Subroutine */ int xermsg_(char *, char *, char *, integer *,
|
||||
integer *, ftnlen, ftnlen, ftnlen);
|
||||
|
||||
/* Fortran I/O blocks */
|
||||
static icilist io___28 = { 0, xern1, 0, "(I8)", 8, 1 };
|
||||
static icilist io___30 = { 0, xern2, 0, "(I8)", 8, 1 };
|
||||
|
||||
|
||||
/* ***BEGIN PROLOGUE PVALUE */
|
||||
/* ***PURPOSE Use the coefficients generated by POLFIT to evaluate the */
|
||||
/* polynomial fit of degree L, along with the first NDER of */
|
||||
/* its derivatives, at a specified point. */
|
||||
/* ***LIBRARY SLATEC */
|
||||
/* ***CATEGORY K6 */
|
||||
/* ***TYPE SINGLE PRECISION (PVALUE-S, DP1VLU-D) */
|
||||
/* ***KEYWORDS CURVE FITTING, LEAST SQUARES, POLYNOMIAL APPROXIMATION */
|
||||
/* ***AUTHOR Shampine, L. F., (SNLA) */
|
||||
/* Davenport, S. M., (SNLA) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* Written by L. F. Shampine and S. M. Davenport. */
|
||||
|
||||
/* Abstract */
|
||||
|
||||
/* The subroutine PVALUE uses the coefficients generated by POLFIT */
|
||||
/* to evaluate the polynomial fit of degree L , along with the first */
|
||||
/* NDER of its derivatives, at a specified point. Computationally */
|
||||
/* stable recurrence relations are used to perform this task. */
|
||||
|
||||
/* The parameters for PVALUE are */
|
||||
|
||||
/* Input -- */
|
||||
/* L - the degree of polynomial to be evaluated. L may be */
|
||||
/* any non-negative integer which is less than or equal */
|
||||
/* to NDEG , the highest degree polynomial provided */
|
||||
/* by POLFIT . */
|
||||
/* NDER - the number of derivatives to be evaluated. NDER */
|
||||
/* may be 0 or any positive value. If NDER is less */
|
||||
/* than 0, it will be treated as 0. */
|
||||
/* X - the argument at which the polynomial and its */
|
||||
/* derivatives are to be evaluated. */
|
||||
/* A - work and output array containing values from last */
|
||||
/* call to POLFIT . */
|
||||
|
||||
/* Output -- */
|
||||
/* YFIT - value of the fitting polynomial of degree L at X */
|
||||
/* YP - array containing the first through NDER derivatives */
|
||||
/* of the polynomial of degree L . YP must be */
|
||||
/* dimensioned at least NDER in the calling program. */
|
||||
|
||||
/* ***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, */
|
||||
/* Curve fitting by polynomials in one variable, Report */
|
||||
/* SLA-74-0270, Sandia Laboratories, June 1974. */
|
||||
/* ***ROUTINES CALLED XERMSG */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 740601 DATE WRITTEN */
|
||||
/* 890531 Changed all specific intrinsics to generic. (WRB) */
|
||||
/* 890531 REVISION DATE from Version 3.2 */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
|
||||
/* 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) */
|
||||
/* 900510 Convert XERRWV calls to XERMSG calls. (RWC) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* ***END PROLOGUE PVALUE */
|
||||
/* ***FIRST EXECUTABLE STATEMENT PVALUE */
|
||||
/* Parameter adjustments */
|
||||
--a;
|
||||
--yp;
|
||||
|
||||
/* Function Body */
|
||||
if (*l < 0) {
|
||||
goto L12;
|
||||
}
|
||||
ndo = max(*nder,0);
|
||||
ndo = min(ndo,*l);
|
||||
maxord = a[1] + .5f;
|
||||
k1 = maxord + 1;
|
||||
k2 = k1 + maxord;
|
||||
k3 = k2 + maxord + 2;
|
||||
nord = a[k3] + .5f;
|
||||
if (*l > nord) {
|
||||
goto L11;
|
||||
}
|
||||
k4 = k3 + *l + 1;
|
||||
if (*nder < 1) {
|
||||
goto L2;
|
||||
}
|
||||
i__1 = *nder;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
/* L1: */
|
||||
yp[i__] = 0.f;
|
||||
}
|
||||
L2:
|
||||
if (*l >= 2) {
|
||||
goto L4;
|
||||
}
|
||||
if (*l == 1) {
|
||||
goto L3;
|
||||
}
|
||||
|
||||
/* L IS 0 */
|
||||
|
||||
val = a[k2 + 1];
|
||||
goto L10;
|
||||
|
||||
/* L IS 1 */
|
||||
|
||||
L3:
|
||||
cc = a[k2 + 2];
|
||||
val = a[k2 + 1] + (*x - a[2]) * cc;
|
||||
if (*nder >= 1) {
|
||||
yp[1] = cc;
|
||||
}
|
||||
goto L10;
|
||||
|
||||
/* L IS GREATER THAN 1 */
|
||||
|
||||
L4:
|
||||
ndp1 = ndo + 1;
|
||||
k3p1 = k3 + 1;
|
||||
k4p1 = k4 + 1;
|
||||
lp1 = *l + 1;
|
||||
lm1 = *l - 1;
|
||||
ilo = k3 + 3;
|
||||
iup = k4 + ndp1;
|
||||
i__1 = iup;
|
||||
for (i__ = ilo; i__ <= i__1; ++i__) {
|
||||
/* L5: */
|
||||
a[i__] = 0.f;
|
||||
}
|
||||
dif = *x - a[lp1];
|
||||
kc = k2 + lp1;
|
||||
a[k4p1] = a[kc];
|
||||
a[k3p1] = a[kc - 1] + dif * a[k4p1];
|
||||
a[k3 + 2] = a[k4p1];
|
||||
|
||||
/* EVALUATE RECURRENCE RELATIONS FOR FUNCTION VALUE AND DERIVATIVES */
|
||||
|
||||
i__1 = lm1;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
in = *l - i__;
|
||||
inp1 = in + 1;
|
||||
k1i = k1 + inp1;
|
||||
ic = k2 + in;
|
||||
dif = *x - a[inp1];
|
||||
val = a[ic] + dif * a[k3p1] - a[k1i] * a[k4p1];
|
||||
if (ndo <= 0) {
|
||||
goto L8;
|
||||
}
|
||||
i__2 = ndo;
|
||||
for (n = 1; n <= i__2; ++n) {
|
||||
k3pn = k3p1 + n;
|
||||
k4pn = k4p1 + n;
|
||||
/* L6: */
|
||||
yp[n] = dif * a[k3pn] + n * a[k3pn - 1] - a[k1i] * a[k4pn];
|
||||
}
|
||||
|
||||
/* SAVE VALUES NEEDED FOR NEXT EVALUATION OF RECURRENCE RELATIONS */
|
||||
|
||||
i__2 = ndo;
|
||||
for (n = 1; n <= i__2; ++n) {
|
||||
k3pn = k3p1 + n;
|
||||
k4pn = k4p1 + n;
|
||||
a[k4pn] = a[k3pn];
|
||||
/* L7: */
|
||||
a[k3pn] = yp[n];
|
||||
}
|
||||
L8:
|
||||
a[k4p1] = a[k3p1];
|
||||
/* L9: */
|
||||
a[k3p1] = val;
|
||||
}
|
||||
|
||||
/* NORMAL RETURN OR ABORT DUE TO ERROR */
|
||||
|
||||
L10:
|
||||
*yfit = val;
|
||||
return 0;
|
||||
|
||||
L11:
|
||||
s_wsfi(&io___28);
|
||||
do_fio(&c__1, (char *)&(*l), (ftnlen)sizeof(integer));
|
||||
e_wsfi();
|
||||
s_wsfi(&io___30);
|
||||
do_fio(&c__1, (char *)&nord, (ftnlen)sizeof(integer));
|
||||
e_wsfi();
|
||||
/* Writing concatenation */
|
||||
i__3[0] = 40, a__1[0] = "THE ORDER OF POLYNOMIAL EVALUATION, L = ";
|
||||
i__3[1] = 8, a__1[1] = xern1;
|
||||
i__3[2] = 49, a__1[2] = " REQUESTED EXCEEDS THE HIGHEST ORDER FIT, NORD "
|
||||
"= ";
|
||||
i__3[3] = 8, a__1[3] = xern2;
|
||||
i__3[4] = 45, a__1[4] = ", COMPUTED BY POLFIT -- EXECUTION TERMINATED.";
|
||||
s_cat(ch__1, a__1, i__3, &c__5, (ftnlen)150);
|
||||
xermsg_("SLATEC", "PVALUE", ch__1, &c__8, &c__2, (ftnlen)6, (ftnlen)6, (
|
||||
ftnlen)150);
|
||||
return 0;
|
||||
|
||||
L12:
|
||||
xermsg_("SLATEC", "PVALUE", "INVALID INPUT PARAMETER. ORDER OF POLYNOMI"
|
||||
"AL EVALUATION REQUESTED IS NEGATIVE -- EXECUTION TERMINATED.", &
|
||||
c__2, &c__2, (ftnlen)6, (ftnlen)6, (ftnlen)103);
|
||||
return 0;
|
||||
} /* pvalue_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
83
ext/f2c_math/xercnt.c
Normal file
83
ext/f2c_math/xercnt.c
Normal file
|
|
@ -0,0 +1,83 @@
|
|||
/* xercnt.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* DECK XERCNT */
|
||||
/* Subroutine */ int xercnt_(char *librar, char *subrou, char *messg, integer
|
||||
*nerr, integer *level, integer *kontrl, ftnlen librar_len, ftnlen
|
||||
subrou_len, ftnlen messg_len)
|
||||
{
|
||||
/* ***BEGIN PROLOGUE XERCNT */
|
||||
/* ***SUBSIDIARY */
|
||||
/* ***PURPOSE Allow user control over handling of errors. */
|
||||
/* ***LIBRARY SLATEC (XERROR) */
|
||||
/* ***CATEGORY R3C */
|
||||
/* ***TYPE ALL (XERCNT-A) */
|
||||
/* ***KEYWORDS ERROR, XERROR */
|
||||
/* ***AUTHOR Jones, R. E., (SNLA) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* Abstract */
|
||||
/* Allows user control over handling of individual errors. */
|
||||
/* Just after each message is recorded, but before it is */
|
||||
/* processed any further (i.e., before it is printed or */
|
||||
/* a decision to abort is made), a call is made to XERCNT. */
|
||||
/* If the user has provided his own version of XERCNT, he */
|
||||
/* can then override the value of KONTROL used in processing */
|
||||
/* this message by redefining its value. */
|
||||
/* KONTRL may be set to any value from -2 to 2. */
|
||||
/* The meanings for KONTRL are the same as in XSETF, except */
|
||||
/* that the value of KONTRL changes only for this message. */
|
||||
/* If KONTRL is set to a value outside the range from -2 to 2, */
|
||||
/* it will be moved back into that range. */
|
||||
|
||||
/* Description of Parameters */
|
||||
|
||||
/* --Input-- */
|
||||
/* LIBRAR - the library that the routine is in. */
|
||||
/* SUBROU - the subroutine that XERMSG is being called from */
|
||||
/* MESSG - the first 20 characters of the error message. */
|
||||
/* NERR - same as in the call to XERMSG. */
|
||||
/* LEVEL - same as in the call to XERMSG. */
|
||||
/* KONTRL - the current value of the control flag as set */
|
||||
/* by a call to XSETF. */
|
||||
|
||||
/* --Output-- */
|
||||
/* KONTRL - the new value of KONTRL. If KONTRL is not */
|
||||
/* defined, it will remain at its original value. */
|
||||
/* This changed value of control affects only */
|
||||
/* the current occurrence of the current message. */
|
||||
|
||||
/* ***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC */
|
||||
/* Error-handling Package, SAND82-0800, Sandia */
|
||||
/* Laboratories, 1982. */
|
||||
/* ***ROUTINES CALLED (NONE) */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 790801 DATE WRITTEN */
|
||||
/* 861211 REVISION DATE from Version 3.2 */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
|
||||
/* 900206 Routine changed from user-callable to subsidiary. (WRB) */
|
||||
/* 900510 Changed calling sequence to include LIBRARY and SUBROUTINE */
|
||||
/* names, changed routine name from XERCTL to XERCNT. (RWC) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* ***END PROLOGUE XERCNT */
|
||||
/* ***FIRST EXECUTABLE STATEMENT XERCNT */
|
||||
return 0;
|
||||
} /* xercnt_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
78
ext/f2c_math/xerhlt.c
Normal file
78
ext/f2c_math/xerhlt.c
Normal file
|
|
@ -0,0 +1,78 @@
|
|||
/* xerhlt.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__9 = 9;
|
||||
static integer c__1 = 1;
|
||||
|
||||
/* DECK XERHLT */
|
||||
/* Subroutine */ int xerhlt_(char *messg, ftnlen messg_len)
|
||||
{
|
||||
/* Builtin functions */
|
||||
integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
|
||||
e_wsle(void);
|
||||
/* Subroutine */ int s_stop(char *, ftnlen);
|
||||
|
||||
/* Fortran I/O blocks */
|
||||
static cilist io___1 = { 0, 6, 0, 0, 0 };
|
||||
|
||||
|
||||
/* ***BEGIN PROLOGUE XERHLT */
|
||||
/* ***SUBSIDIARY */
|
||||
/* ***PURPOSE Abort program execution and print error message. */
|
||||
/* ***LIBRARY SLATEC (XERROR) */
|
||||
/* ***CATEGORY R3C */
|
||||
/* ***TYPE ALL (XERHLT-A) */
|
||||
/* ***KEYWORDS ABORT PROGRAM EXECUTION, ERROR, XERROR */
|
||||
/* ***AUTHOR Jones, R. E., (SNLA) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* Abstract */
|
||||
/* ***Note*** machine dependent routine */
|
||||
/* XERHLT aborts the execution of the program. */
|
||||
/* The error message causing the abort is given in the calling */
|
||||
/* sequence, in case one needs it for printing on a dayfile, */
|
||||
/* for example. */
|
||||
|
||||
/* Description of Parameters */
|
||||
/* MESSG is as in XERMSG. */
|
||||
|
||||
/* ***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC */
|
||||
/* Error-handling Package, SAND82-0800, Sandia */
|
||||
/* Laboratories, 1982. */
|
||||
/* ***ROUTINES CALLED (NONE) */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 790801 DATE WRITTEN */
|
||||
/* 861211 REVISION DATE from Version 3.2 */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
|
||||
/* 900206 Routine changed from user-callable to subsidiary. (WRB) */
|
||||
/* 900510 Changed calling sequence to delete length of character */
|
||||
/* and changed routine name from XERABT to XERHLT. (RWC) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* ***END PROLOGUE XERHLT */
|
||||
/* ***FIRST EXECUTABLE STATEMENT XERHLT */
|
||||
s_wsle(&io___1);
|
||||
do_lio(&c__9, &c__1, "stopping...", (ftnlen)11);
|
||||
e_wsle();
|
||||
s_stop("", (ftnlen)0);
|
||||
return 0;
|
||||
} /* xerhlt_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
475
ext/f2c_math/xermsg.c
Normal file
475
ext/f2c_math/xermsg.c
Normal file
|
|
@ -0,0 +1,475 @@
|
|||
/* xermsg.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__2 = 2;
|
||||
static integer c__0 = 0;
|
||||
static logical c_false = FALSE_;
|
||||
static integer c__4 = 4;
|
||||
static integer c_n1 = -1;
|
||||
static integer c__72 = 72;
|
||||
static integer c__1 = 1;
|
||||
static logical c_true = TRUE_;
|
||||
|
||||
/* DECK XERMSG */
|
||||
/* Subroutine */ int xermsg_(char *librar, char *subrou, char *messg, integer
|
||||
*nerr, integer *level, ftnlen librar_len, ftnlen subrou_len, ftnlen
|
||||
messg_len)
|
||||
{
|
||||
/* System generated locals */
|
||||
address a__1[2];
|
||||
integer i__1, i__2, i__3[2];
|
||||
char ch__1[87];
|
||||
|
||||
/* Builtin functions */
|
||||
/* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
|
||||
integer i_len(char *, ftnlen), s_wsfi(icilist *), do_fio(integer *, char *
|
||||
, ftnlen), e_wsfi(void);
|
||||
/* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, lerr;
|
||||
static char temp[72];
|
||||
extern /* Subroutine */ int fdump_(void);
|
||||
static char xlibr[8];
|
||||
static integer ltemp, kount;
|
||||
static char xsubr[8];
|
||||
extern integer j4save_(integer *, integer *, logical *);
|
||||
static integer llevel, maxmes;
|
||||
static char lfirst[20];
|
||||
extern /* Subroutine */ int xercnt_(char *, char *, char *, integer *,
|
||||
integer *, integer *, ftnlen, ftnlen, ftnlen);
|
||||
static integer lkntrl, kdummy;
|
||||
extern /* Subroutine */ int xerhlt_(char *, ftnlen);
|
||||
static integer mkntrl;
|
||||
extern /* Subroutine */ int xersve_(char *, char *, char *, integer *,
|
||||
integer *, integer *, integer *, ftnlen, ftnlen, ftnlen), xerprn_(
|
||||
char *, integer *, char *, integer *, ftnlen, ftnlen);
|
||||
|
||||
/* Fortran I/O blocks */
|
||||
static icilist io___14 = { 0, temp, 0, "('ERROR NUMBER = ', I8)", 72, 1 };
|
||||
|
||||
|
||||
/* ***BEGIN PROLOGUE XERMSG */
|
||||
/* ***PURPOSE Process error messages for SLATEC and other libraries. */
|
||||
/* ***LIBRARY SLATEC (XERROR) */
|
||||
/* ***CATEGORY R3C */
|
||||
/* ***TYPE ALL (XERMSG-A) */
|
||||
/* ***KEYWORDS ERROR MESSAGE, XERROR */
|
||||
/* ***AUTHOR Fong, Kirby, (NMFECC at LLNL) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* XERMSG processes a diagnostic message in a manner determined by the */
|
||||
/* value of LEVEL and the current value of the library error control */
|
||||
/* flag, KONTRL. See subroutine XSETF for details. */
|
||||
|
||||
/* LIBRAR A character constant (or character variable) with the name */
|
||||
/* of the library. This will be 'SLATEC' for the SLATEC */
|
||||
/* Common Math Library. The error handling package is */
|
||||
/* general enough to be used by many libraries */
|
||||
/* simultaneously, so it is desirable for the routine that */
|
||||
/* detects and reports an error to identify the library name */
|
||||
/* as well as the routine name. */
|
||||
|
||||
/* SUBROU A character constant (or character variable) with the name */
|
||||
/* of the routine that detected the error. Usually it is the */
|
||||
/* name of the routine that is calling XERMSG. There are */
|
||||
/* some instances where a user callable library routine calls */
|
||||
/* lower level subsidiary routines where the error is */
|
||||
/* detected. In such cases it may be more informative to */
|
||||
/* supply the name of the routine the user called rather than */
|
||||
/* the name of the subsidiary routine that detected the */
|
||||
/* error. */
|
||||
|
||||
/* MESSG A character constant (or character variable) with the text */
|
||||
/* of the error or warning message. In the example below, */
|
||||
/* the message is a character constant that contains a */
|
||||
/* generic message. */
|
||||
|
||||
/* CALL XERMSG ('SLATEC', 'MMPY', */
|
||||
/* *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION', */
|
||||
/* *3, 1) */
|
||||
|
||||
/* It is possible (and is sometimes desirable) to generate a */
|
||||
/* specific message--e.g., one that contains actual numeric */
|
||||
/* values. Specific numeric values can be converted into */
|
||||
/* character strings using formatted WRITE statements into */
|
||||
/* character variables. This is called standard Fortran */
|
||||
/* internal file I/O and is exemplified in the first three */
|
||||
/* lines of the following example. You can also catenate */
|
||||
/* substrings of characters to construct the error message. */
|
||||
/* Here is an example showing the use of both writing to */
|
||||
/* an internal file and catenating character strings. */
|
||||
|
||||
/* CHARACTER*5 CHARN, CHARL */
|
||||
/* WRITE (CHARN,10) N */
|
||||
/* WRITE (CHARL,10) LDA */
|
||||
/* 10 FORMAT(I5) */
|
||||
/* CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN// */
|
||||
/* * ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'// */
|
||||
/* * CHARL, 3, 1) */
|
||||
|
||||
/* There are two subtleties worth mentioning. One is that */
|
||||
/* the // for character catenation is used to construct the */
|
||||
/* error message so that no single character constant is */
|
||||
/* continued to the next line. This avoids confusion as to */
|
||||
/* whether there are trailing blanks at the end of the line. */
|
||||
/* The second is that by catenating the parts of the message */
|
||||
/* as an actual argument rather than encoding the entire */
|
||||
/* message into one large character variable, we avoid */
|
||||
/* having to know how long the message will be in order to */
|
||||
/* declare an adequate length for that large character */
|
||||
/* variable. XERMSG calls XERPRN to print the message using */
|
||||
/* multiple lines if necessary. If the message is very long, */
|
||||
/* XERPRN will break it into pieces of 72 characters (as */
|
||||
/* requested by XERMSG) for printing on multiple lines. */
|
||||
/* Also, XERMSG asks XERPRN to prefix each line with ' * ' */
|
||||
/* so that the total line length could be 76 characters. */
|
||||
/* Note also that XERPRN scans the error message backwards */
|
||||
/* to ignore trailing blanks. Another feature is that */
|
||||
/* the substring '$$' is treated as a new line sentinel */
|
||||
/* by XERPRN. If you want to construct a multiline */
|
||||
/* message without having to count out multiples of 72 */
|
||||
/* characters, just use '$$' as a separator. '$$' */
|
||||
/* obviously must occur within 72 characters of the */
|
||||
/* start of each line to have its intended effect since */
|
||||
/* XERPRN is asked to wrap around at 72 characters in */
|
||||
/* addition to looking for '$$'. */
|
||||
|
||||
/* NERR An integer value that is chosen by the library routine's */
|
||||
/* author. It must be in the range -99 to 999 (three */
|
||||
/* printable digits). Each distinct error should have its */
|
||||
/* own error number. These error numbers should be described */
|
||||
/* in the machine readable documentation for the routine. */
|
||||
/* The error numbers need be unique only within each routine, */
|
||||
/* so it is reasonable for each routine to start enumerating */
|
||||
/* errors from 1 and proceeding to the next integer. */
|
||||
|
||||
/* LEVEL An integer value in the range 0 to 2 that indicates the */
|
||||
/* level (severity) of the error. Their meanings are */
|
||||
|
||||
/* -1 A warning message. This is used if it is not clear */
|
||||
/* that there really is an error, but the user's attention */
|
||||
/* may be needed. An attempt is made to only print this */
|
||||
/* message once. */
|
||||
|
||||
/* 0 A warning message. This is used if it is not clear */
|
||||
/* that there really is an error, but the user's attention */
|
||||
/* may be needed. */
|
||||
|
||||
/* 1 A recoverable error. This is used even if the error is */
|
||||
/* so serious that the routine cannot return any useful */
|
||||
/* answer. If the user has told the error package to */
|
||||
/* return after recoverable errors, then XERMSG will */
|
||||
/* return to the Library routine which can then return to */
|
||||
/* the user's routine. The user may also permit the error */
|
||||
/* package to terminate the program upon encountering a */
|
||||
/* recoverable error. */
|
||||
|
||||
/* 2 A fatal error. XERMSG will not return to its caller */
|
||||
/* after it receives a fatal error. This level should */
|
||||
/* hardly ever be used; it is much better to allow the */
|
||||
/* user a chance to recover. An example of one of the few */
|
||||
/* cases in which it is permissible to declare a level 2 */
|
||||
/* error is a reverse communication Library routine that */
|
||||
/* is likely to be called repeatedly until it integrates */
|
||||
/* across some interval. If there is a serious error in */
|
||||
/* the input such that another step cannot be taken and */
|
||||
/* the Library routine is called again without the input */
|
||||
/* error having been corrected by the caller, the Library */
|
||||
/* routine will probably be called forever with improper */
|
||||
/* input. In this case, it is reasonable to declare the */
|
||||
/* error to be fatal. */
|
||||
|
||||
/* Each of the arguments to XERMSG is input; none will be modified by */
|
||||
/* XERMSG. A routine may make multiple calls to XERMSG with warning */
|
||||
/* level messages; however, after a call to XERMSG with a recoverable */
|
||||
/* error, the routine should return to the user. Do not try to call */
|
||||
/* XERMSG with a second recoverable error after the first recoverable */
|
||||
/* error because the error package saves the error number. The user */
|
||||
/* can retrieve this error number by calling another entry point in */
|
||||
/* the error handling package and then clear the error number when */
|
||||
/* recovering from the error. Calling XERMSG in succession causes the */
|
||||
/* old error number to be overwritten by the latest error number. */
|
||||
/* This is considered harmless for error numbers associated with */
|
||||
/* warning messages but must not be done for error numbers of serious */
|
||||
/* errors. After a call to XERMSG with a recoverable error, the user */
|
||||
/* must be given a chance to call NUMXER or XERCLR to retrieve or */
|
||||
/* clear the error number. */
|
||||
/* ***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC */
|
||||
/* Error-handling Package, SAND82-0800, Sandia */
|
||||
/* Laboratories, 1982. */
|
||||
/* ***ROUTINES CALLED FDUMP, J4SAVE, XERCNT, XERHLT, XERPRN, XERSVE */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 880101 DATE WRITTEN */
|
||||
/* 880621 REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988. */
|
||||
/* THERE ARE TWO BASIC CHANGES. */
|
||||
/* 1. A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO */
|
||||
/* PRINT MESSAGES. THIS ROUTINE WILL BREAK LONG MESSAGES */
|
||||
/* INTO PIECES FOR PRINTING ON MULTIPLE LINES. '$$' IS */
|
||||
/* ACCEPTED AS A NEW LINE SENTINEL. A PREFIX CAN BE */
|
||||
/* ADDED TO EACH LINE TO BE PRINTED. XERMSG USES EITHER */
|
||||
/* ' ***' OR ' * ' AND LONG MESSAGES ARE BROKEN EVERY */
|
||||
/* 72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE */
|
||||
/* LENGTH OUTPUT CAN NOW BE AS GREAT AS 76. */
|
||||
/* 2. THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE */
|
||||
/* FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE */
|
||||
/* OF LOWER CASE. */
|
||||
/* 880708 REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30. */
|
||||
/* THE PRINCIPAL CHANGES ARE */
|
||||
/* 1. CLARIFY COMMENTS IN THE PROLOGUES */
|
||||
/* 2. RENAME XRPRNT TO XERPRN */
|
||||
/* 3. REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES */
|
||||
/* SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE / */
|
||||
/* CHARACTER FOR NEW RECORDS. */
|
||||
/* 890706 REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO */
|
||||
/* CLEAN UP THE CODING. */
|
||||
/* 890721 REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN */
|
||||
/* PREFIX. */
|
||||
/* 891013 REVISED TO CORRECT COMMENTS. */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (WRB) */
|
||||
/* 900510 Changed test on NERR to be -9999999 < NERR < 99999999, but */
|
||||
/* NERR .ne. 0, and on LEVEL to be -2 < LEVEL < 3. Added */
|
||||
/* LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and */
|
||||
/* XERCTL to XERCNT. (RWC) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* ***END PROLOGUE XERMSG */
|
||||
/* ***FIRST EXECUTABLE STATEMENT XERMSG */
|
||||
lkntrl = j4save_(&c__2, &c__0, &c_false);
|
||||
maxmes = j4save_(&c__4, &c__0, &c_false);
|
||||
|
||||
/* LKNTRL IS A LOCAL COPY OF THE CONTROL FLAG KONTRL. */
|
||||
/* MAXMES IS THE MAXIMUM NUMBER OF TIMES ANY PARTICULAR MESSAGE */
|
||||
/* SHOULD BE PRINTED. */
|
||||
|
||||
/* WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN */
|
||||
/* CALLING XERMSG. THE ERROR NUMBER SHOULD BE POSITIVE, */
|
||||
/* AND THE LEVEL SHOULD BE BETWEEN 0 AND 2. */
|
||||
|
||||
if (*nerr < -9999999 || *nerr > 99999999 || *nerr == 0 || *level < -1 || *
|
||||
level > 2) {
|
||||
xerprn_(" ***", &c_n1, "FATAL ERROR IN...$$ XERMSG -- INVALID ERROR "
|
||||
"NUMBER OR LEVEL$$ JOB ABORT DUE TO FATAL ERROR.", &c__72, (
|
||||
ftnlen)4, (ftnlen)91);
|
||||
xersve_(" ", " ", " ", &c__0, &c__0, &c__0, &kdummy, (ftnlen)1, (
|
||||
ftnlen)1, (ftnlen)1);
|
||||
xerhlt_(" ***XERMSG -- INVALID INPUT", (ftnlen)27);
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* RECORD THE MESSAGE. */
|
||||
|
||||
i__ = j4save_(&c__1, nerr, &c_true);
|
||||
xersve_(librar, subrou, messg, &c__1, nerr, level, &kount, librar_len,
|
||||
subrou_len, messg_len);
|
||||
|
||||
/* HANDLE PRINT-ONCE WARNING MESSAGES. */
|
||||
|
||||
if (*level == -1 && kount > 1) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* ALLOW TEMPORARY USER OVERRIDE OF THE CONTROL FLAG. */
|
||||
|
||||
s_copy(xlibr, librar, (ftnlen)8, librar_len);
|
||||
s_copy(xsubr, subrou, (ftnlen)8, subrou_len);
|
||||
s_copy(lfirst, messg, (ftnlen)20, messg_len);
|
||||
lerr = *nerr;
|
||||
llevel = *level;
|
||||
xercnt_(xlibr, xsubr, lfirst, &lerr, &llevel, &lkntrl, (ftnlen)8, (ftnlen)
|
||||
8, (ftnlen)20);
|
||||
|
||||
/* Computing MAX */
|
||||
i__1 = -2, i__2 = min(2,lkntrl);
|
||||
lkntrl = max(i__1,i__2);
|
||||
mkntrl = abs(lkntrl);
|
||||
|
||||
/* SKIP PRINTING IF THE CONTROL FLAG VALUE AS RESET IN XERCNT IS */
|
||||
/* ZERO AND THE ERROR IS NOT FATAL. */
|
||||
|
||||
if (*level < 2 && lkntrl == 0) {
|
||||
goto L30;
|
||||
}
|
||||
if (*level == 0 && kount > maxmes) {
|
||||
goto L30;
|
||||
}
|
||||
if (*level == 1 && kount > maxmes && mkntrl == 1) {
|
||||
goto L30;
|
||||
}
|
||||
if (*level == 2 && kount > max(1,maxmes)) {
|
||||
goto L30;
|
||||
}
|
||||
|
||||
/* ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A */
|
||||
/* MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS) */
|
||||
/* AND SENDING IT OUT VIA XERPRN. PRINT ONLY IF CONTROL FLAG */
|
||||
/* IS NOT ZERO. */
|
||||
|
||||
if (lkntrl != 0) {
|
||||
s_copy(temp, "MESSAGE FROM ROUTINE ", (ftnlen)21, (ftnlen)21);
|
||||
/* Computing MIN */
|
||||
i__1 = i_len(subrou, subrou_len);
|
||||
i__ = min(i__1,16);
|
||||
s_copy(temp + 21, subrou, i__, i__);
|
||||
i__1 = i__ + 21;
|
||||
s_copy(temp + i__1, " IN LIBRARY ", i__ + 33 - i__1, (ftnlen)12);
|
||||
ltemp = i__ + 33;
|
||||
/* Computing MIN */
|
||||
i__1 = i_len(librar, librar_len);
|
||||
i__ = min(i__1,16);
|
||||
i__1 = ltemp;
|
||||
s_copy(temp + i__1, librar, ltemp + i__ - i__1, i__);
|
||||
i__1 = ltemp + i__;
|
||||
s_copy(temp + i__1, ".", ltemp + i__ + 1 - i__1, (ftnlen)1);
|
||||
ltemp = ltemp + i__ + 1;
|
||||
xerprn_(" ***", &c_n1, temp, &c__72, (ftnlen)4, ltemp);
|
||||
}
|
||||
|
||||
/* IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE */
|
||||
/* PRINTING THE MESSAGE. THE INTRODUCTORY LINE TELLS THE CHOICE */
|
||||
/* FROM EACH OF THE FOLLOWING THREE OPTIONS. */
|
||||
/* 1. LEVEL OF THE MESSAGE */
|
||||
/* 'INFORMATIVE MESSAGE' */
|
||||
/* 'POTENTIALLY RECOVERABLE ERROR' */
|
||||
/* 'FATAL ERROR' */
|
||||
/* 2. WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE */
|
||||
/* 'PROG CONTINUES' */
|
||||
/* 'PROG ABORTED' */
|
||||
/* 3. WHETHER OR NOT A TRACEBACK WAS REQUESTED. (THE TRACEBACK */
|
||||
/* MAY NOT BE IMPLEMENTED AT SOME SITES, SO THIS ONLY TELLS */
|
||||
/* WHAT WAS REQUESTED, NOT WHAT WAS DELIVERED.) */
|
||||
/* 'TRACEBACK REQUESTED' */
|
||||
/* 'TRACEBACK NOT REQUESTED' */
|
||||
/* NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT */
|
||||
/* EXCEED 74 CHARACTERS. */
|
||||
/* WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED. */
|
||||
|
||||
if (lkntrl > 0) {
|
||||
|
||||
/* THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL. */
|
||||
|
||||
if (*level <= 0) {
|
||||
s_copy(temp, "INFORMATIVE MESSAGE,", (ftnlen)20, (ftnlen)20);
|
||||
ltemp = 20;
|
||||
} else if (*level == 1) {
|
||||
s_copy(temp, "POTENTIALLY RECOVERABLE ERROR,", (ftnlen)30, (
|
||||
ftnlen)30);
|
||||
ltemp = 30;
|
||||
} else {
|
||||
s_copy(temp, "FATAL ERROR,", (ftnlen)12, (ftnlen)12);
|
||||
ltemp = 12;
|
||||
}
|
||||
|
||||
/* THEN WHETHER THE PROGRAM WILL CONTINUE. */
|
||||
|
||||
if (mkntrl == 2 && *level >= 1 || mkntrl == 1 && *level == 2) {
|
||||
i__1 = ltemp;
|
||||
s_copy(temp + i__1, " PROG ABORTED,", ltemp + 14 - i__1, (ftnlen)
|
||||
14);
|
||||
ltemp += 14;
|
||||
} else {
|
||||
i__1 = ltemp;
|
||||
s_copy(temp + i__1, " PROG CONTINUES,", ltemp + 16 - i__1, (
|
||||
ftnlen)16);
|
||||
ltemp += 16;
|
||||
}
|
||||
|
||||
/* FINALLY TELL WHETHER THERE SHOULD BE A TRACEBACK. */
|
||||
|
||||
if (lkntrl > 0) {
|
||||
i__1 = ltemp;
|
||||
s_copy(temp + i__1, " TRACEBACK REQUESTED", ltemp + 20 - i__1, (
|
||||
ftnlen)20);
|
||||
ltemp += 20;
|
||||
} else {
|
||||
i__1 = ltemp;
|
||||
s_copy(temp + i__1, " TRACEBACK NOT REQUESTED", ltemp + 24 - i__1,
|
||||
(ftnlen)24);
|
||||
ltemp += 24;
|
||||
}
|
||||
xerprn_(" ***", &c_n1, temp, &c__72, (ftnlen)4, ltemp);
|
||||
}
|
||||
|
||||
/* NOW SEND OUT THE MESSAGE. */
|
||||
|
||||
xerprn_(" * ", &c_n1, messg, &c__72, (ftnlen)4, messg_len);
|
||||
|
||||
/* IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER AND REQUEST A */
|
||||
/* TRACEBACK. */
|
||||
|
||||
if (lkntrl > 0) {
|
||||
s_wsfi(&io___14);
|
||||
do_fio(&c__1, (char *)&(*nerr), (ftnlen)sizeof(integer));
|
||||
e_wsfi();
|
||||
for (i__ = 16; i__ <= 22; ++i__) {
|
||||
if (*(unsigned char *)&temp[i__ - 1] != ' ') {
|
||||
goto L20;
|
||||
}
|
||||
/* L10: */
|
||||
}
|
||||
|
||||
L20:
|
||||
/* Writing concatenation */
|
||||
i__3[0] = 15, a__1[0] = temp;
|
||||
i__3[1] = 23 - (i__ - 1), a__1[1] = temp + (i__ - 1);
|
||||
s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)87);
|
||||
xerprn_(" * ", &c_n1, ch__1, &c__72, (ftnlen)4, 23 - (i__ - 1) + 15);
|
||||
fdump_();
|
||||
}
|
||||
|
||||
/* IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE. */
|
||||
|
||||
if (lkntrl != 0) {
|
||||
xerprn_(" * ", &c_n1, " ", &c__72, (ftnlen)4, (ftnlen)1);
|
||||
xerprn_(" ***", &c_n1, "END OF MESSAGE", &c__72, (ftnlen)4, (ftnlen)
|
||||
14);
|
||||
xerprn_(" ", &c__0, " ", &c__72, (ftnlen)4, (ftnlen)1);
|
||||
}
|
||||
|
||||
/* IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE */
|
||||
/* CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN. */
|
||||
|
||||
L30:
|
||||
if (*level <= 0 || *level == 1 && mkntrl <= 1) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A */
|
||||
/* FATAL ERROR. PRINT THE REASON FOR THE ABORT AND THE ERROR */
|
||||
/* SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT. */
|
||||
|
||||
if (lkntrl > 0 && kount < max(1,maxmes)) {
|
||||
if (*level == 1) {
|
||||
xerprn_(" ***", &c_n1, "JOB ABORT DUE TO UNRECOVERED ERROR.", &
|
||||
c__72, (ftnlen)4, (ftnlen)35);
|
||||
} else {
|
||||
xerprn_(" ***", &c_n1, "JOB ABORT DUE TO FATAL ERROR.", &c__72, (
|
||||
ftnlen)4, (ftnlen)29);
|
||||
}
|
||||
xersve_(" ", " ", " ", &c_n1, &c__0, &c__0, &kdummy, (ftnlen)1, (
|
||||
ftnlen)1, (ftnlen)1);
|
||||
xerhlt_(" ", (ftnlen)1);
|
||||
} else {
|
||||
xerhlt_(messg, messg_len);
|
||||
}
|
||||
return 0;
|
||||
} /* xermsg_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
316
ext/f2c_math/xerprn.c
Normal file
316
ext/f2c_math/xerprn.c
Normal file
|
|
@ -0,0 +1,316 @@
|
|||
/* xerprn.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__4 = 4;
|
||||
static integer c__1 = 1;
|
||||
|
||||
/* DECK XERPRN */
|
||||
/* Subroutine */ int xerprn_(char *prefix, integer *npref, char *messg,
|
||||
integer *nwrap, ftnlen prefix_len, ftnlen messg_len)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer i__1, i__2;
|
||||
|
||||
/* Builtin functions */
|
||||
integer i_len(char *, ftnlen);
|
||||
/* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
|
||||
integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void),
|
||||
i_indx(char *, char *, ftnlen, ftnlen), s_cmp(char *, char *,
|
||||
ftnlen, ftnlen);
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, n, iu[5];
|
||||
static char cbuff[148];
|
||||
static integer lpref, nextc, lwrap, nunit;
|
||||
extern integer i1mach_(integer *);
|
||||
static integer lpiece, idelta, lenmsg;
|
||||
extern /* Subroutine */ int xgetua_(integer *, integer *);
|
||||
|
||||
/* Fortran I/O blocks */
|
||||
static cilist io___9 = { 0, 0, 0, "(A)", 0 };
|
||||
static cilist io___13 = { 0, 0, 0, "(A)", 0 };
|
||||
|
||||
|
||||
/* ***BEGIN PROLOGUE XERPRN */
|
||||
/* ***SUBSIDIARY */
|
||||
/* ***PURPOSE Print error messages processed by XERMSG. */
|
||||
/* ***LIBRARY SLATEC (XERROR) */
|
||||
/* ***CATEGORY R3C */
|
||||
/* ***TYPE ALL (XERPRN-A) */
|
||||
/* ***KEYWORDS ERROR MESSAGES, PRINTING, XERROR */
|
||||
/* ***AUTHOR Fong, Kirby, (NMFECC at LLNL) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* This routine sends one or more lines to each of the (up to five) */
|
||||
/* logical units to which error messages are to be sent. This routine */
|
||||
/* is called several times by XERMSG, sometimes with a single line to */
|
||||
/* print and sometimes with a (potentially very long) message that may */
|
||||
/* wrap around into multiple lines. */
|
||||
|
||||
/* PREFIX Input argument of type CHARACTER. This argument contains */
|
||||
/* characters to be put at the beginning of each line before */
|
||||
/* the body of the message. No more than 16 characters of */
|
||||
/* PREFIX will be used. */
|
||||
|
||||
/* NPREF Input argument of type INTEGER. This argument is the number */
|
||||
/* of characters to use from PREFIX. If it is negative, the */
|
||||
/* intrinsic function LEN is used to determine its length. If */
|
||||
/* it is zero, PREFIX is not used. If it exceeds 16 or if */
|
||||
/* LEN(PREFIX) exceeds 16, only the first 16 characters will be */
|
||||
/* used. If NPREF is positive and the length of PREFIX is less */
|
||||
/* than NPREF, a copy of PREFIX extended with blanks to length */
|
||||
/* NPREF will be used. */
|
||||
|
||||
/* MESSG Input argument of type CHARACTER. This is the text of a */
|
||||
/* message to be printed. If it is a long message, it will be */
|
||||
/* broken into pieces for printing on multiple lines. Each line */
|
||||
/* will start with the appropriate prefix and be followed by a */
|
||||
/* piece of the message. NWRAP is the number of characters per */
|
||||
/* piece; that is, after each NWRAP characters, we break and */
|
||||
/* start a new line. In addition the characters '$$' embedded */
|
||||
/* in MESSG are a sentinel for a new line. The counting of */
|
||||
/* characters up to NWRAP starts over for each new line. The */
|
||||
/* value of NWRAP typically used by XERMSG is 72 since many */
|
||||
/* older error messages in the SLATEC Library are laid out to */
|
||||
/* rely on wrap-around every 72 characters. */
|
||||
|
||||
/* NWRAP Input argument of type INTEGER. This gives the maximum size */
|
||||
/* piece into which to break MESSG for printing on multiple */
|
||||
/* lines. An embedded '$$' ends a line, and the count restarts */
|
||||
/* at the following character. If a line break does not occur */
|
||||
/* on a blank (it would split a word) that word is moved to the */
|
||||
/* next line. Values of NWRAP less than 16 will be treated as */
|
||||
/* 16. Values of NWRAP greater than 132 will be treated as 132. */
|
||||
/* The actual line length will be NPREF + NWRAP after NPREF has */
|
||||
/* been adjusted to fall between 0 and 16 and NWRAP has been */
|
||||
/* adjusted to fall between 16 and 132. */
|
||||
|
||||
/* ***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC */
|
||||
/* Error-handling Package, SAND82-0800, Sandia */
|
||||
/* Laboratories, 1982. */
|
||||
/* ***ROUTINES CALLED I1MACH, XGETUA */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 880621 DATE WRITTEN */
|
||||
/* 880708 REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF */
|
||||
/* JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK */
|
||||
/* THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE */
|
||||
/* SLASH CHARACTER IN FORMAT STATEMENTS. */
|
||||
/* 890706 REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO */
|
||||
/* STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK */
|
||||
/* LINES TO BE PRINTED. */
|
||||
/* 890721 REVISED TO ADD A NEW FEATURE. A NEGATIVE VALUE OF NPREF */
|
||||
/* CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH. */
|
||||
/* 891013 REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH. */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (WRB) */
|
||||
/* 900510 Added code to break messages between words. (RWC) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* ***END PROLOGUE XERPRN */
|
||||
/* ***FIRST EXECUTABLE STATEMENT XERPRN */
|
||||
xgetua_(iu, &nunit);
|
||||
|
||||
/* A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD */
|
||||
/* ERROR MESSAGE UNIT INSTEAD. I1MACH(4) RETRIEVES THE STANDARD */
|
||||
/* ERROR MESSAGE UNIT. */
|
||||
|
||||
n = i1mach_(&c__4);
|
||||
i__1 = nunit;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
if (iu[i__ - 1] == 0) {
|
||||
iu[i__ - 1] = n;
|
||||
}
|
||||
/* L10: */
|
||||
}
|
||||
|
||||
/* LPREF IS THE LENGTH OF THE PREFIX. THE PREFIX IS PLACED AT THE */
|
||||
/* BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING */
|
||||
/* THE REST OF THIS ROUTINE. */
|
||||
|
||||
if (*npref < 0) {
|
||||
lpref = i_len(prefix, prefix_len);
|
||||
} else {
|
||||
lpref = *npref;
|
||||
}
|
||||
lpref = min(16,lpref);
|
||||
if (lpref != 0) {
|
||||
s_copy(cbuff, prefix, lpref, prefix_len);
|
||||
}
|
||||
|
||||
/* LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE */
|
||||
/* TIME FROM MESSG TO PRINT ON ONE LINE. */
|
||||
|
||||
/* Computing MAX */
|
||||
i__1 = 16, i__2 = min(132,*nwrap);
|
||||
lwrap = max(i__1,i__2);
|
||||
|
||||
/* SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS. */
|
||||
|
||||
lenmsg = i_len(messg, messg_len);
|
||||
n = lenmsg;
|
||||
i__1 = n;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
if (*(unsigned char *)&messg[lenmsg - 1] != ' ') {
|
||||
goto L30;
|
||||
}
|
||||
--lenmsg;
|
||||
/* L20: */
|
||||
}
|
||||
L30:
|
||||
|
||||
/* IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE. */
|
||||
|
||||
if (lenmsg == 0) {
|
||||
i__1 = lpref;
|
||||
s_copy(cbuff + i__1, " ", lpref + 1 - i__1, (ftnlen)1);
|
||||
i__1 = nunit;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
io___9.ciunit = iu[i__ - 1];
|
||||
s_wsfe(&io___9);
|
||||
do_fio(&c__1, cbuff, lpref + 1);
|
||||
e_wsfe();
|
||||
/* L40: */
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING */
|
||||
/* STARTS. FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL. */
|
||||
/* WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT. */
|
||||
/* WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED. */
|
||||
|
||||
/* WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL. THE */
|
||||
/* INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE */
|
||||
/* OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH */
|
||||
/* OF THE SECOND ARGUMENT. */
|
||||
|
||||
/* THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE */
|
||||
/* FOLLOWING ORDER. WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER */
|
||||
/* OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT */
|
||||
/* POSITION NEXTC. */
|
||||
|
||||
/* LPIECE .EQ. 0 THE NEW LINE SENTINEL DOES NOT OCCUR IN THE */
|
||||
/* REMAINDER OF THE CHARACTER STRING. LPIECE */
|
||||
/* SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC, */
|
||||
/* WHICHEVER IS LESS. */
|
||||
|
||||
/* LPIECE .EQ. 1 THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC: */
|
||||
/* NEXTC). LPIECE IS EFFECTIVELY ZERO, AND WE */
|
||||
/* PRINT NOTHING TO AVOID PRODUCING UNNECESSARY */
|
||||
/* BLANK LINES. THIS TAKES CARE OF THE SITUATION */
|
||||
/* WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF */
|
||||
/* EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE */
|
||||
/* SENTINEL FOLLOWED BY MORE CHARACTERS. NEXTC */
|
||||
/* SHOULD BE INCREMENTED BY 2. */
|
||||
|
||||
/* LPIECE .GT. LWRAP+1 REDUCE LPIECE TO LWRAP. */
|
||||
|
||||
/* ELSE THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1 */
|
||||
/* RESET LPIECE = LPIECE-1. NOTE THAT THIS */
|
||||
/* PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ. */
|
||||
/* LWRAP+1. THAT IS, THE SENTINEL FALLS EXACTLY */
|
||||
/* AT THE END OF A LINE. */
|
||||
|
||||
nextc = 1;
|
||||
L50:
|
||||
lpiece = i_indx(messg + (nextc - 1), "$$", lenmsg - (nextc - 1), (ftnlen)
|
||||
2);
|
||||
if (lpiece == 0) {
|
||||
|
||||
/* THERE WAS NO NEW LINE SENTINEL FOUND. */
|
||||
|
||||
idelta = 0;
|
||||
/* Computing MIN */
|
||||
i__1 = lwrap, i__2 = lenmsg + 1 - nextc;
|
||||
lpiece = min(i__1,i__2);
|
||||
if (lpiece < lenmsg + 1 - nextc) {
|
||||
for (i__ = lpiece + 1; i__ >= 2; --i__) {
|
||||
i__1 = nextc + i__ - 2;
|
||||
if (s_cmp(messg + i__1, " ", nextc + i__ - 1 - i__1, (ftnlen)
|
||||
1) == 0) {
|
||||
lpiece = i__ - 1;
|
||||
idelta = 1;
|
||||
goto L54;
|
||||
}
|
||||
/* L52: */
|
||||
}
|
||||
}
|
||||
L54:
|
||||
i__1 = lpref;
|
||||
s_copy(cbuff + i__1, messg + (nextc - 1), lpref + lpiece - i__1,
|
||||
nextc + lpiece - 1 - (nextc - 1));
|
||||
nextc = nextc + lpiece + idelta;
|
||||
} else if (lpiece == 1) {
|
||||
|
||||
/* WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1). */
|
||||
/* DON'T PRINT A BLANK LINE. */
|
||||
|
||||
nextc += 2;
|
||||
goto L50;
|
||||
} else if (lpiece > lwrap + 1) {
|
||||
|
||||
/* LPIECE SHOULD BE SET DOWN TO LWRAP. */
|
||||
|
||||
idelta = 0;
|
||||
lpiece = lwrap;
|
||||
for (i__ = lpiece + 1; i__ >= 2; --i__) {
|
||||
i__1 = nextc + i__ - 2;
|
||||
if (s_cmp(messg + i__1, " ", nextc + i__ - 1 - i__1, (ftnlen)1) ==
|
||||
0) {
|
||||
lpiece = i__ - 1;
|
||||
idelta = 1;
|
||||
goto L58;
|
||||
}
|
||||
/* L56: */
|
||||
}
|
||||
L58:
|
||||
i__1 = lpref;
|
||||
s_copy(cbuff + i__1, messg + (nextc - 1), lpref + lpiece - i__1,
|
||||
nextc + lpiece - 1 - (nextc - 1));
|
||||
nextc = nextc + lpiece + idelta;
|
||||
} else {
|
||||
|
||||
/* IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1. */
|
||||
/* WE SHOULD DECREMENT LPIECE BY ONE. */
|
||||
|
||||
--lpiece;
|
||||
i__1 = lpref;
|
||||
s_copy(cbuff + i__1, messg + (nextc - 1), lpref + lpiece - i__1,
|
||||
nextc + lpiece - 1 - (nextc - 1));
|
||||
nextc = nextc + lpiece + 2;
|
||||
}
|
||||
|
||||
/* PRINT */
|
||||
|
||||
i__1 = nunit;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
io___13.ciunit = iu[i__ - 1];
|
||||
s_wsfe(&io___13);
|
||||
do_fio(&c__1, cbuff, lpref + lpiece);
|
||||
e_wsfe();
|
||||
/* L60: */
|
||||
}
|
||||
|
||||
if (nextc <= lenmsg) {
|
||||
goto L50;
|
||||
}
|
||||
return 0;
|
||||
} /* xerprn_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
239
ext/f2c_math/xersve.c
Normal file
239
ext/f2c_math/xersve.c
Normal file
|
|
@ -0,0 +1,239 @@
|
|||
/* xersve.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__4 = 4;
|
||||
static integer c__1 = 1;
|
||||
|
||||
/* DECK XERSVE */
|
||||
/* Subroutine */ int xersve_(char *librar, char *subrou, char *messg, integer
|
||||
*kflag, integer *nerr, integer *level, integer *icount, ftnlen
|
||||
librar_len, ftnlen subrou_len, ftnlen messg_len)
|
||||
{
|
||||
/* Initialized data */
|
||||
|
||||
static integer kountx = 0;
|
||||
static integer nmsg = 0;
|
||||
|
||||
/* Format strings */
|
||||
static char fmt_9000[] = "(\0020 ERROR MESSAGE SUMMARY\002/\002"
|
||||
" LIBRARY SUBROUTINE MESSAGE START NERR\002,\002 "
|
||||
" LEVEL COUNT\002)";
|
||||
static char fmt_9010[] = "(1x,a,3x,a,3x,a,3i10)";
|
||||
static char fmt_9020[] = "(\0020OTHER ERRORS NOT INDIVIDUALLY TABULATED "
|
||||
"= \002,i10)";
|
||||
static char fmt_9030[] = "(1x)";
|
||||
|
||||
/* System generated locals */
|
||||
integer i__1, i__2;
|
||||
|
||||
/* Builtin functions */
|
||||
integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);
|
||||
/* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
|
||||
integer s_cmp(char *, char *, ftnlen, ftnlen);
|
||||
|
||||
/* Local variables */
|
||||
static integer i__;
|
||||
static char lib[8], mes[20], sub[8];
|
||||
static integer lun[5], iunit, kunit, nunit, kount[10];
|
||||
extern integer i1mach_(integer *);
|
||||
static char libtab[8*10], mestab[20*10];
|
||||
static integer nertab[10], levtab[10];
|
||||
static char subtab[8*10];
|
||||
extern /* Subroutine */ int xgetua_(integer *, integer *);
|
||||
|
||||
/* Fortran I/O blocks */
|
||||
static cilist io___7 = { 0, 0, 0, fmt_9000, 0 };
|
||||
static cilist io___9 = { 0, 0, 0, fmt_9010, 0 };
|
||||
static cilist io___16 = { 0, 0, 0, fmt_9020, 0 };
|
||||
static cilist io___17 = { 0, 0, 0, fmt_9030, 0 };
|
||||
|
||||
|
||||
/* ***BEGIN PROLOGUE XERSVE */
|
||||
/* ***SUBSIDIARY */
|
||||
/* ***PURPOSE Record that an error has occurred. */
|
||||
/* ***LIBRARY SLATEC (XERROR) */
|
||||
/* ***CATEGORY R3 */
|
||||
/* ***TYPE ALL (XERSVE-A) */
|
||||
/* ***KEYWORDS ERROR, XERROR */
|
||||
/* ***AUTHOR Jones, R. E., (SNLA) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* *Usage: */
|
||||
|
||||
/* INTEGER KFLAG, NERR, LEVEL, ICOUNT */
|
||||
/* CHARACTER * (len) LIBRAR, SUBROU, MESSG */
|
||||
|
||||
/* CALL XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, ICOUNT) */
|
||||
|
||||
/* *Arguments: */
|
||||
|
||||
/* LIBRAR :IN is the library that the message is from. */
|
||||
/* SUBROU :IN is the subroutine that the message is from. */
|
||||
/* MESSG :IN is the message to be saved. */
|
||||
/* KFLAG :IN indicates the action to be performed. */
|
||||
/* when KFLAG > 0, the message in MESSG is saved. */
|
||||
/* when KFLAG=0 the tables will be dumped and */
|
||||
/* cleared. */
|
||||
/* when KFLAG < 0, the tables will be dumped and */
|
||||
/* not cleared. */
|
||||
/* NERR :IN is the error number. */
|
||||
/* LEVEL :IN is the error severity. */
|
||||
/* ICOUNT :OUT the number of times this message has been seen, */
|
||||
/* or zero if the table has overflowed and does not */
|
||||
/* contain this message specifically. When KFLAG=0, */
|
||||
/* ICOUNT will not be altered. */
|
||||
|
||||
/* *Description: */
|
||||
|
||||
/* Record that this error occurred and possibly dump and clear the */
|
||||
/* tables. */
|
||||
|
||||
/* ***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC */
|
||||
/* Error-handling Package, SAND82-0800, Sandia */
|
||||
/* Laboratories, 1982. */
|
||||
/* ***ROUTINES CALLED I1MACH, XGETUA */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 800319 DATE WRITTEN */
|
||||
/* 861211 REVISION DATE from Version 3.2 */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
|
||||
/* 900413 Routine modified to remove reference to KFLAG. (WRB) */
|
||||
/* 900510 Changed to add LIBRARY NAME and SUBROUTINE to calling */
|
||||
/* sequence, use IF-THEN-ELSE, make number of saved entries */
|
||||
/* easily changeable, changed routine name from XERSAV to */
|
||||
/* XERSVE. (RWC) */
|
||||
/* 910626 Added LIBTAB and SUBTAB to SAVE statement. (BKS) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* ***END PROLOGUE XERSVE */
|
||||
/* ***FIRST EXECUTABLE STATEMENT XERSVE */
|
||||
|
||||
if (*kflag <= 0) {
|
||||
|
||||
/* Dump the table. */
|
||||
|
||||
if (nmsg == 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Print to each unit. */
|
||||
|
||||
xgetua_(lun, &nunit);
|
||||
i__1 = nunit;
|
||||
for (kunit = 1; kunit <= i__1; ++kunit) {
|
||||
iunit = lun[kunit - 1];
|
||||
if (iunit == 0) {
|
||||
iunit = i1mach_(&c__4);
|
||||
}
|
||||
|
||||
/* Print the table header. */
|
||||
|
||||
io___7.ciunit = iunit;
|
||||
s_wsfe(&io___7);
|
||||
e_wsfe();
|
||||
|
||||
/* Print body of table. */
|
||||
|
||||
i__2 = nmsg;
|
||||
for (i__ = 1; i__ <= i__2; ++i__) {
|
||||
io___9.ciunit = iunit;
|
||||
s_wsfe(&io___9);
|
||||
do_fio(&c__1, libtab + (i__ - 1 << 3), (ftnlen)8);
|
||||
do_fio(&c__1, subtab + (i__ - 1 << 3), (ftnlen)8);
|
||||
do_fio(&c__1, mestab + (i__ - 1) * 20, (ftnlen)20);
|
||||
do_fio(&c__1, (char *)&nertab[i__ - 1], (ftnlen)sizeof(
|
||||
integer));
|
||||
do_fio(&c__1, (char *)&levtab[i__ - 1], (ftnlen)sizeof(
|
||||
integer));
|
||||
do_fio(&c__1, (char *)&kount[i__ - 1], (ftnlen)sizeof(integer)
|
||||
);
|
||||
e_wsfe();
|
||||
/* L10: */
|
||||
}
|
||||
|
||||
/* Print number of other errors. */
|
||||
|
||||
if (kountx != 0) {
|
||||
io___16.ciunit = iunit;
|
||||
s_wsfe(&io___16);
|
||||
do_fio(&c__1, (char *)&kountx, (ftnlen)sizeof(integer));
|
||||
e_wsfe();
|
||||
}
|
||||
io___17.ciunit = iunit;
|
||||
s_wsfe(&io___17);
|
||||
e_wsfe();
|
||||
/* L20: */
|
||||
}
|
||||
|
||||
/* Clear the error tables. */
|
||||
|
||||
if (*kflag == 0) {
|
||||
nmsg = 0;
|
||||
kountx = 0;
|
||||
}
|
||||
} else {
|
||||
|
||||
/* PROCESS A MESSAGE... */
|
||||
/* SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG, */
|
||||
/* OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL. */
|
||||
|
||||
s_copy(lib, librar, (ftnlen)8, librar_len);
|
||||
s_copy(sub, subrou, (ftnlen)8, subrou_len);
|
||||
s_copy(mes, messg, (ftnlen)20, messg_len);
|
||||
i__1 = nmsg;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
if (s_cmp(lib, libtab + (i__ - 1 << 3), (ftnlen)8, (ftnlen)8) ==
|
||||
0 && s_cmp(sub, subtab + (i__ - 1 << 3), (ftnlen)8, (
|
||||
ftnlen)8) == 0 && s_cmp(mes, mestab + (i__ - 1) * 20, (
|
||||
ftnlen)20, (ftnlen)20) == 0 && *nerr == nertab[i__ - 1] &&
|
||||
*level == levtab[i__ - 1]) {
|
||||
++kount[i__ - 1];
|
||||
*icount = kount[i__ - 1];
|
||||
return 0;
|
||||
}
|
||||
/* L30: */
|
||||
}
|
||||
|
||||
if (nmsg < 10) {
|
||||
|
||||
/* Empty slot found for new message. */
|
||||
|
||||
++nmsg;
|
||||
s_copy(libtab + (i__ - 1 << 3), lib, (ftnlen)8, (ftnlen)8);
|
||||
s_copy(subtab + (i__ - 1 << 3), sub, (ftnlen)8, (ftnlen)8);
|
||||
s_copy(mestab + (i__ - 1) * 20, mes, (ftnlen)20, (ftnlen)20);
|
||||
nertab[i__ - 1] = *nerr;
|
||||
levtab[i__ - 1] = *level;
|
||||
kount[i__ - 1] = 1;
|
||||
*icount = 1;
|
||||
} else {
|
||||
|
||||
/* Table is full. */
|
||||
|
||||
++kountx;
|
||||
*icount = 0;
|
||||
}
|
||||
}
|
||||
return 0;
|
||||
|
||||
/* Formats. */
|
||||
|
||||
} /* xersve_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
93
ext/f2c_math/xgetua.c
Normal file
93
ext/f2c_math/xgetua.c
Normal file
|
|
@ -0,0 +1,93 @@
|
|||
/* xgetua.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#ifdef _cpluscplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__5 = 5;
|
||||
static integer c__0 = 0;
|
||||
static logical c_false = FALSE_;
|
||||
|
||||
/* DECK XGETUA */
|
||||
/* Subroutine */ int xgetua_(integer *iunita, integer *n)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer i__1;
|
||||
|
||||
/* Local variables */
|
||||
static integer i__, index;
|
||||
extern integer j4save_(integer *, integer *, logical *);
|
||||
|
||||
/* ***BEGIN PROLOGUE XGETUA */
|
||||
/* ***PURPOSE Return unit number(s) to which error messages are being */
|
||||
/* sent. */
|
||||
/* ***LIBRARY SLATEC (XERROR) */
|
||||
/* ***CATEGORY R3C */
|
||||
/* ***TYPE ALL (XGETUA-A) */
|
||||
/* ***KEYWORDS ERROR, XERROR */
|
||||
/* ***AUTHOR Jones, R. E., (SNLA) */
|
||||
/* ***DESCRIPTION */
|
||||
|
||||
/* Abstract */
|
||||
/* XGETUA may be called to determine the unit number or numbers */
|
||||
/* to which error messages are being sent. */
|
||||
/* These unit numbers may have been set by a call to XSETUN, */
|
||||
/* or a call to XSETUA, or may be a default value. */
|
||||
|
||||
/* Description of Parameters */
|
||||
/* --Output-- */
|
||||
/* IUNIT - an array of one to five unit numbers, depending */
|
||||
/* on the value of N. A value of zero refers to the */
|
||||
/* default unit, as defined by the I1MACH machine */
|
||||
/* constant routine. Only IUNIT(1),...,IUNIT(N) are */
|
||||
/* defined by XGETUA. The values of IUNIT(N+1),..., */
|
||||
/* IUNIT(5) are not defined (for N .LT. 5) or altered */
|
||||
/* in any way by XGETUA. */
|
||||
/* N - the number of units to which copies of the */
|
||||
/* error messages are being sent. N will be in the */
|
||||
/* range from 1 to 5. */
|
||||
|
||||
/* ***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC */
|
||||
/* Error-handling Package, SAND82-0800, Sandia */
|
||||
/* Laboratories, 1982. */
|
||||
/* ***ROUTINES CALLED J4SAVE */
|
||||
/* ***REVISION HISTORY (YYMMDD) */
|
||||
/* 790801 DATE WRITTEN */
|
||||
/* 861211 REVISION DATE from Version 3.2 */
|
||||
/* 891214 Prologue converted to Version 4.0 format. (BAB) */
|
||||
/* 920501 Reformatted the REFERENCES section. (WRB) */
|
||||
/* ***END PROLOGUE XGETUA */
|
||||
/* ***FIRST EXECUTABLE STATEMENT XGETUA */
|
||||
/* Parameter adjustments */
|
||||
--iunita;
|
||||
|
||||
/* Function Body */
|
||||
*n = j4save_(&c__5, &c__0, &c_false);
|
||||
i__1 = *n;
|
||||
for (i__ = 1; i__ <= i__1; ++i__) {
|
||||
index = i__ + 4;
|
||||
if (i__ == 1) {
|
||||
index = 3;
|
||||
}
|
||||
iunita[i__] = j4save_(&index, &c__0, &c_false);
|
||||
/* L30: */
|
||||
}
|
||||
return 0;
|
||||
} /* xgetua_ */
|
||||
|
||||
#ifdef _cpluscplus
|
||||
}
|
||||
#endif
|
||||
Loading…
Add table
Reference in a new issue