e.cpp

一个非常有用的开源代码

/* ----------------------- MODULE feigen.cpp ------------------------ *
 * Reference: "Numerical Algorithms with C by G. Engeln-Mueller and   *
 *             F. Uhlig, Springer-Verlag, 1996"                       *
 * ------------------------------------------------------------------ */




// Header file of basis.cpp

/***********************************************************************
*                                                                      *
* Basic functions: Declarations file (with types and macros)           *
* ----------------------------------------------------------           *
*                                                                      *
* Programming language: ANSI C                                         *
* Compiler:             Turbo C 2.0                                    *
* Computer:             IBM PS/2 70 with 80387                         *
* Author:               Juergen Dietel, Computer Center, RWTH Aachen   *
* Date:                 9.30.1992                                      *
*                                                                      *
***********************************************************************/



/***********************************************************************
* Make preparations in case this declarations file is included more    *
* than once in a source code                                           *
***********************************************************************/

#ifndef BASIS_H_INCLUDED
#define BASIS_H_INCLUDED



/***********************************************************************
* Include a few other declarations file here.                          *
* We include the standard declarations files here mainly so that only  *
* this file needs to be included in the C modules of the Numerical     *
* Library in order to have access to all standard names and variables. *
* Moreover modifications for nonstandard compilers may then be limited *
* to this file only.                                                   *
***********************************************************************/

#ifdef __EMX__            /* emx 0.9b with GNU CC 2.7.2?              */
#if __GNUC__ == 2         /* To prevent a compiler crash in           */
#if __GNUC_MINOR__ == 7   /* `15\gax.c' while compiling with LDOUBLE, */
#define EMX09B            /* by defining the macro _WITH_UNDERSCORE   */
#define _WITH_UNDERSCORE  /* the special variants of the mathematical */
#define fabsl _fabsl      /* functions for `long double' in `math.h'  */
#define sqrtl _sqrtl      /* are activated and their use further      */
#define powl  _powl       /* below in shape of the macros FABS, SQRT, */
#define sinl  _sinl       /* ... is prepared here; which seems quite  */
#define cosl  _cosl       /* useful even without the affair with the  */
#define expl  _expl       /* crashing compiler.                       */
#define logl  _logl
#define atanl _atanl
#define acosl _acosl
#define coshl _coshl
#endif
#endif
#endif

#include <stdlib.h>
#include <stdio.h>
#include <math.h>     /*  for  fabs, sqrt, pow, exp, sin, cos, log,   */
                      /*       atan, acos                             */
#include <float.h>    /*  for  DBL_EPSILON, DBL_MAX */ 
#include <string.h>
                 
#ifdef sun            /*  for  GNU CC under SunOS                     */
#include <unistd.h>   /*  for  SEEK_END                               */
#else
#ifdef amigados       /*  for  GNU CC on an Amiga                     */
#include <unistd.h>   /*  for  SEEK_END                               */
#endif
#endif
#ifndef SEEK_END      /* SEEK_END not yet defined (such as for        */
                      /* GNU CC 2.1 for an i386 MS-DOS computer)?     */                            
#endif



/***********************************************************************
* Predefine the desired precision for floating point arithmetic:       *
* If the macro FLOAT has been defined, we compute in single precision  *
* (data type float), for LDOUBLE we compute in maximal precision       *
* (data type long double), otherwise we use double precision (data type*
* double). LDOUBLE only adds precision on those compilers for which    *
* long double is different from double such as on Turbo C, but not on  *
* QuickC compilers for example.                                        *
* For input and output of floating numbers the macros LZS (length      *
* prefix for scanf()) and LZP (length prefix for printf()) should be   *
* used (obviously two different ones), since according to the ANSI C   *
* Standard the output of double values should be done without the      *
* length declaration "l" in contrast to the input where this "l" is    *
* needed.                                                              *
* NOTE:    If the compiler used does not know the macros FLT_MAX,      *
* ====     LDBL_MAX, DBL_MAX or FLT_MAX_EXP, LDBL_MAX_EXP, DBL_MAX_EXP *
*          (which can usually be found in float.h), the user must      *
*          insert suitable values at the locations marked with !!!!.   *
***********************************************************************/

#ifdef FLOAT                     /* single precision ? ...............*/

typedef float     REAL;          /* Standard floating point type float*/
/*.IX{REAL}*/
typedef double    LONG_REAL;     /* more precise floating point type  */
/*.IX{LONG\unt REAL}*/

#ifdef FLT_EPSILON               /* ANSI C compiler ? ................*/
#define MACH_EPS  (REAL)FLT_EPSILON
/*.IX{MACH\unt EPS}*/
#else                            /* not an ANSI C compiler ? .........*/
#define MACH_EPS  mach_eps()     /* machine constant .................*/
#endif

#ifdef FLT_MAX_EXP               /* ANSI C compiler? .................*/
#define MAX_EXP   FLT_MAX_EXP    /* Binary exponent of POSMAX ........*/
/*.IX{MAX\unt EXP}*/
#else                            /* not an ANSI C compiler ? .........*/
#define MAX_EXP   128            /* make adjustments here !!!! .......*/
#endif

#ifdef FLT_MAX                   /* ANSI C compiler? ................ */
#define POSMAX    (REAL)FLT_MAX  /* largest floating point number ... */
/*.IX{POS\unt MAX}*/
#else                            /* not an ANSI C compiler ? ........ */
#define POSMAX    1e38f          /* adjust here !!!! ................ */
#endif

#ifdef FLT_MIN                   /* ANSI C compiler? .................*/
#define POSMIN    (REAL)FLT_MIN  /* smallest positive floating point  */
/*.IX{POS\unt MIN}*/
                                 /* number ...........................*/
#else                            /* not an ANSI C compiler ? .........*/
#define POSMIN    posmin()
#endif

#define LZS       ""             /* length prefix for formatted       */
/*.IX{LZS}*/
                                 /* input of floating point numbers ..*/
#define LZP       ""             /* length prefix for formatted       */
/*.IX{LZP}*/
                                 /* output of floating point numbers .*/


#else
#ifdef LDOUBLE                   /* maximal precision ? ..............*/

typedef long double  REAL;       /* Standard floating point type      */
                                 /* long double                       */
/*.IX{REAL}*/
typedef long double  LONG_REAL;  /* "more precise" floating point type*/
/*.IX{LONG\unt REAL}*/
#define LONG_DOUBLE_USED

#ifdef LDBL_EPSILON              /* ANSI C compiler? .................*/
#define MACH_EPS  (REAL)LDBL_EPSILON
/*.IX{MACH\unt EPS}*/
#else                            /* not an ANSI C compiler ? .........*/
#define MACH_EPS  mach_eps()     /* machine constant .................*/
#endif

#ifdef LDBL_MAX_EXP              /* ANSI C compiler? .................*/
#define MAX_EXP   LDBL_MAX_EXP   /* Binary exponent of POSMAX ........*/
/*.IX{MAX\unt EXP}*/
#else                            /* not an ANSI C compiler ? .........*/
#define MAX_EXP   1023           /* adjust here !!!! .................*/
#endif

#ifdef LDBL_MAX                  /* ANSI C compiler? .................*/
#define POSMAX    (REAL)LDBL_MAX /* largest floating point number ....*/
/*.IX{POS\unt MAX}*/
#else                            /* not an ANSI C compiler ? .........*/
#define POSMAX    1e100l         /* adjust here !!!! .................*/
#endif

#ifdef LDBL_MIN                  /* ANSI C compiler? .................*/
#define POSMIN    (REAL)LDBL_MIN /* smallest positive floating point  */
/*.IX{POS\unt MIN}*/                /* number ...........................*/
#else                            /* not an  ANSI C compiler ? ........*/
#define POSMIN    posmin()
#endif

#define LZS       "L"            /* length prefix for formatted       */
/*.IX{LZS}*/
                                 /* input of floating point numbers ..*/
#define LZP       "L"            /* length prefix for formatted       */
/*.IX{LZP}*/
                                 /* output of floating point numbers .*/


#else                            /* double precision ? ...............*/

typedef double       REAL;       /* Standard floating point type      */
                                 /* double                            */
/*.IX{REAL}*/
typedef long double  LONG_REAL;  /* more precise floating point type  */
/*.IX{LONG\unt REAL}*/

#ifdef DBL_EPSILON               /* ANSI C compiler? .................*/
#define MACH_EPS  (REAL)DBL_EPSILON
/*.IX{MACH\unt EPS}*/
#else                            /* not an ANSI C compiler ? .........*/
#define MACH_EPS  mach_eps()     /* machine constant .................*/
#endif

#ifdef DBL_MAX_EXP               /* ANSI C compiler? .................*/
#define MAX_EXP   DBL_MAX_EXP    /* Binary exponent of POSMAX ........*/
/*.IX{MAX\unt EXP}*/
#else                            /* not an ANSI C compiler ? .........*/
#define MAX_EXP   1023           /* adjust here !!!! .................*/
#endif

#ifdef DBL_MAX                   /* ANSI C compiler? .................*/
#define POSMAX    (REAL)DBL_MAX  /* largest floating point number ....*/
/*.IX{POS\unt MAX}*/
#else                            /* not an ANSI C compiler ? .........*/
#define POSMAX    1e100          /* adjust here !!!! .................*/
#endif

#ifdef DBL_MIN                   /* ANSI C compiler? .................*/
#ifdef __BORLANDC__
#if __BORLANDC__ <= 0x0200       /* Borland C++ 2.0 for DOS? .........*/
                                           /* because the value       */
#define POSMIN    2.2250738585072017E-308  /* 2.2250738585072014E-308 */
#else                                      /* from `float.h' is       */
                                           /* regarded as zero!       */
#define POSMIN    DBL_MIN        /* smallest positive floating point  */
#endif                           /* number ...........................*/
#else
#define POSMIN    DBL_MIN        /* smallest positive floating point  */
#endif
/*.IX{POS\unt MIN}*/                /* number ...........................*/
#else                            /* not an ANSI C compiler ? .........*/
#define POSMIN    posmin()
#endif

#define LZS       "l"            /* length prefix for formatted       */
/*.IX{LZS}*/
                                 /* input of floating point numbers ..*/
#define LZP       ""             /* length prefix for formatted       */
/*.IX{LZP}*/
                                 /* output of floating point numbers .*/
#endif
#endif



/***********************************************************************
* declare several important data types                                 *
***********************************************************************/

#undef boolean
#undef FALSE
#undef TRUE
typedef enum {FALSE, TRUE} boolean;
/*.IX{FALSE}*/
/*.IX{TRUE}*/
/*.IX{boolean}*/

/* function pointer types for approximation in chapter 8 .............*/
typedef REAL (*ansatzfnk) (int i, REAL x);
/*.IX{ansatzfnk}*/
typedef REAL (*approxfnk) (REAL c[], REAL x);
/*.IX{approxfnk}*/
typedef void (*ableitfnk) (REAL x, REAL c[], REAL *d);
/*.IX{ableitfnk}*/

/* Type of function, which evaluates the right hand side of an .......*/
/* explicit ordinary differential equation  y' = f(x,y)        .......*/
typedef REAL (*dglfnk)(REAL x, REAL y);
/*.IX{dglfnk}*/


/* Type of function, which evaluates the right hand side .............*/
/* of a system of differential equations  y' = f(x,y)    .............*/
typedef void (*dglsysfnk)(REAL x, REAL y[], REAL f[]);
/*.IX{dglsysfnk}*/

/* Type of function, which computes the boundary value r(ya, yb) .....*/
/* of a two-point first order boundary value problem             .....*/
typedef void (*rndbedfnk)(REAL ya[], REAL yb[], REAL r[]);
/*.IX{rndbedfnk}*/

/* enumeration type for classification of error codes that are  ......*/
/* returned from most of the functions that realize a numerical ......*/
/* method                                                       ......*/
typedef enum { KEIN_FEHLER, WARNUNG, UNBEKANNT, FATAL } fehler_t;
/*.IX{KEIN\unt FEHLER}*/
/*.IX{WARNUNG}*/
/*.IX{UNBEKANNT}*/
/*.IX{FATAL}*/
/*.IX{fehler\unt t}*/

typedef REAL abl_mat1[4][2];     /* used to evaluate splinefunctions  */
/*.IX{abl\unt mat1}*/