ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng)/ProfitPress-MegaCDROM2.B6I/PROG/MISC/LIBRY51.ZIP/LIBRY4B.DOC

Text File | 1989-11-10 | 20.0 KB | 480 lines

.pa
                 QUICK LIST OF MATHEMATICAL SUBROUTINES

BFNLQ.... brute force method for nonlinear simultaneous equations (see BRYDN)
BFNLQD... double precision version of BFNLQ (see BRYDN)
BRYDN.... modified Broyden's method for nonlinear simultaneous equations
BRYDND... double precision version of BRYDN
CONJG.... conjugate gradient method for nonlinear simult. equns. (see BRYDN)
CONJGD... double precision version of CONJG (see BRYDN)
CUBIC.... solve a cubic equation
CUBICD... double precision form of CUBIC
FMIN..... find the minimum of a function within an interval
FMIND.... double precision version of FMIN
FZERO.... find the zero (root) of a function within an interval
FZEROD... double precision version of FZERO
GAUSP.... Gauss-Jordan elimination of simultaneous equations
GAUSPD... double precision form of GAUSP
JLDAY.... convert month/day/year to Julian (annual) day
LLSLC.... linear least-squares with linear constraints
LNREG.... linear regression
MADD..... matrix addition
MADDD.... double precision matrix addition
MCPY..... matrix copy
MCPYD.... double precision matrix copy
MDAY..... convert Julian (annual) day to month/day/year
MINV..... matrix inversion
MINVD.... double precision form of MINV
MPRD..... matrix product (multiplication)
MPRDD.... double precision matrix product (multiplication)
MSUB..... matrix subtraction
MSUBD.... double precision matrix subtraction
MTRA..... matrix transpose
MTRAD.... double precision matrix transpose
MTRS..... in-place transpose of a square matrix
MTRSD.... in-place transpose of a double precision square matrix
NTNLQ.... Newton's method for nonlinear simultaneous equations (see BRYDN)
NTNLQD... double precision version of NTNLQ (see BRYDN)
PROOT.... polynomial root finder
QUAD..... solve a quadratic
QUADD.... double precision form of QUAD
RAND..... random number generator (real)
RK4...... 4th order Runge-Kutta (solution of ordinary differential eqns.)
RK6D..... 6th order Runge-Kutta (solution of ordinary differential eqns.)
SPLNE.... determine coefficients for a cubic spline
SVD...... singular value decomposition using Householder transformations
TDIAG.... tridiagonal matrix solver
TDIAGD... double precision form of TDIAG
TREND.... compute trend of periodic function

NOTE: A word of caution about matrix operation subroutines.  All of these
      routines (MADD, MSUB, MCPY, MPRD, MSUB, MTRA, MTRS and their double
      precision  counterparts)  use vector emulations.  Subsequently,  no
      arrays  can  cross  a segment boundary.  See  Chapter  7  for  more
      cautions.
.pa
                           MATHEMATICAL SUBROUTINES


NAME:     BRYDN
PURPOSE:  modified Broyden's method for nonlinear simultaneous equations
TYPE:     subroutine (far external)
SYNTAX:   CALL BRYDN(USER,XMIN,XMAX,X,F,N,M,WORK,MW,MCALL,IPRT,IER)
IMPUT:    USER (far external) a subroutine that YOU MUST PROVIDE and
          must be called by the sequence CALL USER(X,F)
          XMIN,XMAX (REAL*4) arrays, upper and lower limits on X
          N (INTEGER*2) number of elements in X
          M (INTEGER*2) number of residuals
          W (REAL*4) working space of dimension MW
          MW (INTEGER*2) dimension of working space MW=5N+3M+N*N+N*M
          MCALL (INTEGER*2) maximum function calls (set to zero for unlimited)
          IPRT (INTEGER*2) print select (for IPRT=0 nothing is printed,
          for IPRT=1 only a one line summary at the end is printed, for
          IPRT=2 X, F, and G are printed at each iteration, and for IPRT=3
          the Jacobian is also printed at each iteration)
          need to open file 6 on the PC)
OUTPUT:   X (REAL*4) solution vector
          F (REAL*4) final residual
          IER (INTEGER*2) error indicator
               IER=0  indicates normal return
               IER=1  indicates N<2
               IER=2  indicates M<N
               IER=3  indicates XMIN>XMAX conflict
               IER=4  indicates insufficient work space (increase MW)
          note that the printer output is going to go directly to the
          printer unless you first spool it with by CALLing
          SPOOL('FILE.EXT',IER); and no, the redirect ">" will not work
NOTE:     This is a FAR more difficult problem than most people even
          imagine!  You can compare it to finding the lowest spot on
          the face of the Earth given the rules of the game "Battle
          Ship" where all you can do is call out coordinates and your
          opponent answers with the elevation.  You may find a rut or a
          ditch somewhere; but it's highly unlikely that you will find
          death valley, let alone some trench in the Pacific. Now extend
          this analogy to many-dimensional space... get the point?  So
          don't get too upset if BRYDN can't find the solution you want.

          I've tried all sorts of algorithms and this seems to work the
          best for general problems.

          If you have 6 equations and 6 unknowns then M=N=6.  If you
          have 12 equations and 6 unknowns then M=12 and N=6.

          If all of your equations are linear, by all means don't use BRYDN,
          you want LLSLC.  If you have only one unknown then use FMIN.

          For double precision use BRYDND.

          Also available are the brute force method (BFNLQ & BFNLQD), the
          conjugate gradient method (CONJG & CONJGD), and Newton's method
          (NTNLQ & NTNLQD). The calling sequence, input, and output are
          all the same as for BRYDN & BRYDND. The only differences is the
          method and the required working space, MW (for BFNLQ, MW=3*N+M;
          for CONJG, MW=8*N+M; and for NTNLQ, MW=6*N+2*M+N*N+N*M).

          This is a little tricky, so I have included an example at the end of
          this section.


NAME:     CUBIC
PURPOSE:  solve a cubic equation
TYPE:     subroutine (far external)
SYNTAX:   CALL SUBROUTINE CUBIC(A1,A2,A3,A4,R1,C1,R2,C2,R3,C3)
INPUT:    A1,A2,A3,A4 (REAL*4) coefficients in cubic equation
OUTPUT:   (R1,C1),(R2,C2),(R3,C3) (REAL*4) roots
NOTE:     the equation has the form A1*X**3+A2*X**2+A3*X+A4=0
          for double precision use CUBICD and don't forget to declare
          A1,A2,A3,A4,R1,C1,R2,C2,R3,C3 all to be REAL*8 and if you
          use constants, don't forget the 1.D0 or whatever.
NOTE:     for double precision use CUBICD


NAME:     FMIN
PURPOSE:  find the minimum of a function within an interval
TYPE:     REAL*4 function (far external)
SYNTAX:   XMIN=FMIN(F,A,B)
INPUT:    A,B (REAL*4) the search interval
          F (a far external REAL*4 function) that you must supply
          that may be called by FX=F(X)
OUTPUT:   location of minimum
NOTE:     for double precision use FMIND


NAME:     FZERO
PURPOSE:  find the zero (root) of a function within an interval
TYPE:     REAL*4 function (far external)
SYNTAX:   X0=FZERO(F,A,B)
INPUT:    A,B (REAL*4) interval to look for root in
          F (REAL*4 function) that YOU MUST SUPPLY and must be
          called by the sequence FX=F(X)
OUTPUT:   location of the root
NOTE:     for double precision use FZEROD


NAME:     GAUSP
PURPOSE:  Gauss-Jordan elimination of simultaneous equations
TYPE:     subroutine (far external)
SYNTAX:   CALL GAUSP(A,B,X,JPIVOT,N,D,C,IER)
INPUT:    A,B (REAL*4) arrays containing the equations to be solved
          having dimensions (N,N) and (N) respectively
          JPIVOT (INTEGER*2) working space of dimension "N"
          N (INTEGER*2) the number of equations
OUTPUT:   X (REAL*4) solution vector of dimension N
          D (REAL*4) ALOG10(ABS(DETERMINANT)) log of the determinant
          C (REAL*4) ALOG10(AMAX1(ABS(PIVOT))/AMIN1(ABS(PIVOT))) a
          measure of the gradedness or condition of the matrix
          IER (INTEGER*2) error indicator (IER=0 is normal)
NOTE:     For double precision use GAUSPD and don't forget to declare
          A,B,X,D,C all to be REAL*8.
          Note that A and B will be destroyed upon return.
          Note that A is in row-order as indicated below

               A(1,1)*X(1)+A(1,2)*X(2)+...A(1,N)*X(N)=B(1)
               A(2,1)*X(1)+A(2,2)*X(2)+...A(2,N)*X(N)=B(2)
               A(3,1)*X(1)+A(3,2)*X(2)+...A(3,N)*X(N)=B(3)
               ..........................................
               A(N,1)*X(1)+A(N,2)*X(2)+...A(N,N)*X(N)=B(N)

          When you really have to be careful about this is with DATA
          statements as FORTRAN fills arrays in column-order.

          Full column pivoting is used.
          Vector emulation is used throughout to give maximum speed.

NOTE:     for double precision use GAUSPD


NAME:     JLDAY
PURPOSE:  convert month/day/year to Julian (annual) day
TYPE:     subroutine (far external)
SYNTAX:   CALL JLDAY(MONTH,IDAY,IYEAR,JDAY)
INPUT:    MONTH,IDAY,IYEAR (INTEGER*2)
OUTPUT:   JDAY (INTEGER*2)
NOTE:     I know that "Julian" day is a misnomer; but that's what
          everyone else calls it.


NAME:     LLSLC
PURPOSE:  linear least-squares with linear constraints
TYPE:     subroutine (far external)
SYNTAX:   CALL LLSLC(X,Y,A,NE,NC,W,NW,IOPT,IER)
INPUT:    X,Y,A,W (REAL*8) see example below
          NE (INTEGER*2) number of variables (X or A)
          NC (INTEGER*2) number of constraints (may be zero)
          NW (INTEGER*2) working space, NW>=(NE+NC)*(NE+NC+3))
          IOPT (INTEGER*2) option, see example below
OUTPUT:   A (REAL*8)
          IER (INTEGER*2) error indicator
               IER=0   normal
               IER=1   dimension error
               IER=2   singular matrix - no inverse
               IER=3   too few constants+constraints (minimum=2)
               IER=4   too few equations (minimum=1)
               IER=5   too few constraints (minimum=0)
               IER=6   too many constraints (maximum=ne-1)
               IER=7   IOPT not defined
               IER=8   matrices not initialized
               IER=9   too few points for least-squares
               IER=10  constraints added before least-squares points
               IER=11  too many constraints were added
               IER=12  too few constraints were specified
               IER=13  working space too small (increase NW)
NOTE:     This is a little tricky, so I have included an example at the end of
          this section. Sorry, there is no REAL*4 equivalent.


NAME:     LNREG
PURPOSE:  linear regression
TYPE:     subroutine (far external)
SYNTAX:   CALL LNREG(X,Y,N,A,B,R,IER)
INPUT:    X,Y (REAL*4) data points
          N (INTEGER*2) number of points
OUTPUT:   A,B (REAL*4) slope, intercept
          R (REAL*4) regression coefficient in %
          IER (INTEGER*2) error indicator (IER=0 is normal)
NOTE:     LNREG only fits a straight line of the form Y=A*X+B to a set
          of points.  If you want something more elaborate use LLSLC.


NAME:     MADD
PURPOSE:  matrix addition
TYPE:     subroutine (far external)
SYNTAX:   CALL MADD(A,B,C,N,M)
INPUT:    A (REAL*4) array of dimension A(N,M)
          B (REAL*4) array of dimension B(N,M)
          N (INTEGER*2) number of rows
          M (INTEGER*2) number of columns
OUTPUT:   C (REAL*4) array of dimension C(N,M)
NOTE:     for double precision use MADDD


NAME:     MCPY
PURPOSE:  matrix copy
TYPE:     subroutine (far external)
SYNTAX:   CALL MCPY(A,B,N,M)
INPUT:    A (REAL*4) array of dimension A(N,M)
          N (INTEGER*2) number of rows
          M (INTEGER*2) number of columns
OUTPUT:   B (REAL*4) array of dimension B(N,M)
NOTE:     for double precision use MCPYD


NAME:     MDAY
PURPOSE:  convert Julian (annual) day to month/day/year
TYPE:     subroutine (far external)
SYNTAX:   CALL MDAY(JDAY,IYEAR,MONTH,IDAY)
INPUT:    JDAY,IYEAR (INTEGER*2)
OUTPUT:   MONTH,IDAY (INTEGER*2)
NOTE:     I know that "Julian" day is a misnomer; but that's what
          everyone else calls it.


NAME:     MINV
PURPOSE:  matrix inversion
TYPE:     subroutine (far external)
SYNTAX:   CALL MINV(A,IPIVOT,JPIVOT,N,D,C,IER)
INPUT:    A (REAL*4) array containing the matrix to be inverted
          IPIVOT,JPIVOT (INTEGER*2) working spaces of dimension "N"
          N (INTEGER*2) rank (the number of equations)
OUTPUT:   A (REAL*4) inverted matrix
          D (REAL*4) ALOG10(ABS(DETERMINANT)) log of the determinant
          C (REAL*4) ALOG10(AMAX1(ABS(PIVOT))/AMIN1(ABS(PIVOT))) a
          measure of the gradedness or condition of the matrix
          IER (INTEGER*2) error indicator (IER=0 is normal)
NOTE:     For double precision use MINVD and don't forget to declare
          A,D,C all to be REAL*8.
          Note that A is inverted "in place."
          Note that A is in row-order (see GAUSP).


NAME:     MPRD
PURPOSE:  matrix product (multiply)
TYPE:     subroutine (far external)
SYNTAX:   CALL MPRD(A,B,C,N,M,L)
INPUT:    A (REAL*4) array of dimension A(N,M)
          B (REAL*4) array of dimension B(M,L)
          N (INTEGER*2) number of rows in A and C
          M (INTEGER*2) number of columns in A and rows in B
          L (INTEGER*2) number of columns in B and C
OUTPUT:   C (REAL*4) array of dimension C(N,L)
NOTE:     for double precision use MPRDD


NAME:     MSUB
PURPOSE:  matrix subtraction
TYPE:     subroutine (far external)
SYNTAX:   CALL MSUB(A,B,C,N,M)
INPUT:    A (REAL*4) array of dimension A(N,M)
          B (REAL*4) array of dimension B(N,M)
          N (INTEGER*2) number of rows
          M (INTEGER*2) number of columns
OUTPUT:   C (REAL*4) array of dimension C(N,M)
NOTE:     for double precision use MSUBD


NAME:     MTRA
PURPOSE:  matrix transpose
TYPE:     subroutine (far external)
SYNTAX:   CALL MTRA(A,B,N,M)
INPUT:    A (REAL*4) array of dimension A(N,M)
          N (INTEGER*2) number of rows in A and columns in B
          M (INTEGER*2) number of columns in A and rows in B
OUTPUT:   B (REAL*4) array of dimension B(M,N)
NOTE:     for double precision use MTRAD


NAME:     MTRS
PURPOSE:  in-place transpose of a square matrix
TYPE:     subroutine (far external)
SYNTAX:   CALL MTRS(A,N)
INPUT:    A (REAL*4) array of dimension A(N,N)
          N (INTEGER*2) number of rows and columns in A
OUTPUT:   A (REAL*4) array of dimension A(N,N)
NOTE:     for double precision use MTRSD


NAME:     PROOT
PURPOSE:  polynomial root finder
TYPE:     subroutine (far external)
SYNTAX:   CALL PROOT(P,R,C,N,IER)
INPUT:    P (REAL*8) array containing the polynomial coefficients
          N (INTEGER*2) number of terms in P
OUTPUT:   R (REAL*8) array containing the real part of the roots
          C (REAL*8) array containing the complex part of the roots
          IER (INTEGER*2) error indicator
               IER=0 indicates no error
               IER=1 indicates input error
               IER=2 indicates failure to converge
NOTE:     Bairstow's method is used.
          The polynomial must be of the form

               P1+P2*X+P3*X**2+P4*X**3+....=0

          Also note that this is a VERY difficult problem for large N,
          so don't be too upset if you get complex roots when you
          thought you should be getting real roots.  That's life in the
          real world.


NAME:     QUAD
PURPOSE:  solve a quadratic
TYPE:     subroutine (far external)
SYNTAX:   CALL SUBROUTINE QUAD(A1,A2,A3,R1,C1,R2,C2)
INPUT:    A1,A2,A3 (REAL*4) coefficients in quadratic equation
OUTPUT:   (R1,C1),(R2,C2) (REAL*4) roots
NOTE:     the equation has the form A1*X**2+A2*X+A3=0
          for double precision use QUADD and don't forget to declare
          A1,A2,A3,R1,C1,R2,C2 all to be REAL*8 and if you use
          constants, don't forget the 1.D0 or whatever.
NOTE:     for double precision use QUADD


NAME:     RAND
PURPOSE:  return a random number
TYPE:     subroutine (far external)
SYNTAX:   CALL RAND(X)
OUTPUT:   X (REAL*4) a random number between 0 and 1
NOTE:     this uses the clock for a seed and the "old IBM" algorithm which
          has come under fire as of late, but it works just fine for my
          purposes


NAME:     RK4
PURPOSE:  4th order Runge-Kutta (solve ordinary differential equations)
TYPE:     subroutine (far external)
SYNTAX:   CALL RK4(USER,X,DX,Y,DY,N,W,V)
INPUT:    X (REAL*4) independent variable
          DX (REAL*4) increment in X (step size)
          Y (REAL*4) array, dependent variables at X
          DY (REAL*4) array, change in Y
          N (INTEGER*2) number of dependent variables (Y)
          W (REAL*4) working space of dimension (N,4)
          V (REAL*4) working space of dimension (N)
          USER (far external) a subroutine that YOU MUST PROVIDE and
          must be called by the sequence CALL USER(X,Y,DY)
OUTPUT:   Y (REAL*4) new values of dependent variable
NOTE:     this isn't too tricky; but I have provided an example that may
          help


NAME:     RK6D
PURPOSE:  6th order Runge-Kutta (solve ordinary differential equations)
TYPE:     subroutine (far external)
SYNTAX:   CALL RK6D(USER,X,DX,Y,DY,N,W,V)
INPUT:    X (REAL*8) independent variable
          DX (REAL*8) increment in X (step size)
          Y (REAL*8) array, dependent variables at X
          DY (REAL*8) array, change in Y
          N (INTEGER*2) number of dependent variables (Y)
          W (REAL*8) working space of dimension (N,8) <-- note that the
          "8" here is not a misprint - it takes 8 steps to get 6th order
          even though it only takes 4 to get 4-order Runge-Kutta
          V (REAL*4) working space of dimension (N)
          USER (far external) a subroutine that YOU MUST PROVIDE and
          must be called by the sequence CALL USER(X,Y,DY)
OUTPUT:   Y (REAL*8) new values of dependent variable
NOTE:     This is just like RK4 except for the "8" above and that you
          need to declare everything REAL*8 instead of REAL*4.


NAME:     SPLNE
PURPOSE:  determine coefficients for a cubic spline
TYPE:     subroutine (far external)
SYNTAX:   CALL SPLNE(XI,YI,N,C)
INPUT:    XI,YI (REAL*4) points to be fit
          N (INTEGER*2) number of points
OUTPUT:   C (REAL*4) coefficient array of dimension (3,N)
NOTE:     you need to call this once, then use SEVAL to evaluate spline


NAME:     SVD
PURPOSE:  singular value decomposition using Householder transformations
TYPE:     subroutine (far external)
SYNTAX:   CALL SVD(A,N,M,S,U,V,W,IER)
INPUT:    A (REAL*8) the matrix (N,M) to be decomposed - row order
          N (INTEGER*2) number of rows in A
          M (INTEGER*2) number of columns in A
          W (REAL*8) working space of dimension MAX0(N,M)
OUTPUT:   S (REAL*8) singular values
          U (REAL*8) left eigenvectors (M,N)
          V (REAL*8) right eigenvectors (N,M)
          IER (INTEGER*2) error indicator (IER=0 is normal)


NAME:     TDIAG
PURPOSE:  tridiagonal matrix solver
TYPE:     subroutine (far external)
SYNTAX:   CALL TDIAG(D,B,X,W,N,IER)
INPUT:    D (REAL*4) matrix of dimension (N,3)
          B (REAL*4) right-hand-side of dimension (N)
          W (REAL*4) working space of dimension (N,3)
          N (INTEGER*2) number of equations
OUTPUT:   X (REAL*4) solution vector of dimension (N)
          IER (INTEGER*2) error indicator (IER=0 is normal)
NOTE:     the tridiagonal system of equations is defined by  [D][X]=[B]
          [D] should look something like the following

                 0  2 -2
                -1  2 -1
                -1  2 -1
                 .  .  .
                -2  2  0

          For double precision use TDIAGD and don't forget to declare
          D,B,X,W all to be of type REAL*8.


NAME:     TREND
PURPOSE:  compute trend of periodic function
TYPE:     subroutine (far external)
SYNTAX:   CALL TREND(X,Y,W,NX,N,C,NC)
INPUT:    X,Y (REAL*4) points to be fit
          W (REAL*4) working space of dimension (NX)
          NX (INTEGER*2) number of points
          N (INTEGER*2) number of known points (N<NX)
          NC (INTEGER*2) number of terms in expansion (try 6)
OUTPUT:   Y (REAL*4) upon return, the N+1,N+2,N+3,...NX values of Y will
          be trended
          C (REAL*4) coefficients in trend series
NOTE:     Chebyshev approximation is used