Programmer 7500

home *** CD-ROM | disk | FTP | other *** search

/ Programmer 7500 / MAX_PROGRAMMERS.iso / CLIPPER / MISC / LIBRY31A.ZIP / LIBRY4.DOC < prev next >

Wrap

Text File | 1987-01-20 | 40.9 KB | 1,126 lines

.pa MATHEMATICAL These mathematical procedures provide a wide variety of functions that are useful for engineering applications. I began developing this library after my first experience with SLIM (name changed to protect the guilty) back in 1976. So I guess I have them to thank. If you've ever used SLIM perhaps you've had a similar experience... Let's say you want to solve what one would think to be a simple problem of simultaneous linear equations, and due to some unfortunate turn of events, you are given SLIM with which to do this. You are confronted with no less than 60 subroutines from which to choose. As it turns out, none will do what you want. Anyway, you narrow it down to the Smith-Jones-Simpson-Wilkerson-Feinstein-Gidley method for partially symmetric matrices with strange sparceness stored in convoluted- compacted form (not to be confused with space-saving form) and the Kleidowitz-Yonkers-Lakeville-Massey-Perkins modification of the Orwell- Gauss-Newton-Schmeltz algorithm for ill-conditioned Hilbert matrices stored in space-saving form (not to be confused with convoluted- compacted form). An obvious choice indeed (any fool with 6 Ph.D.s in math could see that the latter is the preferred method). After numerous attempts to link your program with all 347 of the necessary SLIM libraries you are finally ready to run your program. Of course, it bombs out because you were reading tuesday afternoon's manual and it's now wednesday morning. Tough luck... Actually, SLIM is so bad that even I can't achieve hyperbole. Language fails me. The above could easily be a random selection from one of the six volumes that are supposed to be a "user's guide." Anyway, I have a completely different philosophy of programming. I hope that it works well for you: if there is a clearly preferred method of doing anything, then use it and discard the others. An excellent example of this is Gauss quadrature. There are only two methods of numerical integration that will (at least in theory) converge given an infinite order: trapezoidal rule and Gauss quadrature. Gauss quadrature is so far superior in accuracy to any other method (e.g. Newton-Cotes) it's practically incredible. Therefore, why have 50 other methods from which to choose? An interesting note about this is that, contrary to some rumors 10 subdivisions of 20-point Gauss quadrature (200 point evaluation) is not as accurate as 96-point Gauss quadrature with less than half the function evaluations. Another example is Gauss-Jordan elimination. SLIM has subroutines that provide a "high accuracy solution" (If Andy Rooney understood this, he'd probably ask, "Who would want a low accuracy solution anyway?"). And of course, there's the iterative improvement type solutions. What they don't tell you is that the residual must be computed in greater precision than the matrix is reduced in (say single to double precision) which is OK; but why not just use double precision all the way through? If you are already using double precision and that's not good enough for you... well, you're just out of luck. If your matrix is ill- conditioned (as is typically the case with Vandermondes - that's what you get when you're curve-fitting) about the best thing you can use is Gauss-Jordan elimination with full column pivoting. It is fast enough and I have used vector emulation throughout. I have tried all sorts of things like QR decompositions, Givens rotations, and a whole bunch of others; but it all boils down to the precision of the math coprocessor. I only use one algorithm for solving full matrices. I have tested it using a series of standard Hilbert matrices against SLIM's "high accuracy solution" and found mine to be both faster and more accurate. One last example is the solution of ordinary differential equations. If you read numerical mathematics texts, they always give examples of the error in such-and-such a method as compared to the exact solution. If you look closely at the microscopic print in the margin of the legend under the fly wing that got stuck on the copier, you will probably find the pithy little statement "exact solution determined using 4-th order Runge-Kutta with 0.0000000001 step size." Runge-Kutta is self-starting (which is a bigger pain in the neck than you might think if you use a method that isn't), convenient to use, and works very well. So why not use it, I ask? Of course there are certain (unknown before the fact) cases where other methods are much faster; but I would rather spend less time experimenting with zillions of algorithms and let the machine sweat a little more. I've been through it once; and once is enough. Runge-Kutta is the one for me. A word about automatic stepsize control. It has been my experience that this takes so long at each step that you are better off to use a smaller step from the start. I had two such subroutines that I spent a lot of time on; but I purged them one day in a fit of frustration. One final note: WATCH YOUR VARIABLE TYPES! ESPECIALLY DOUBLE PRECISION! I've seen many a programmer bewildered by a program that crashes because they have passed "0." instead of "0.D0" to a subroutine expecting a double precision value or forgotten to put the "REAL*8" out in front of a function. REAL*8 FUNCTION DOFX(D) It is really best to use IMPLICIT statements like IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z) than to individually type each variable. You're bound to miss one. Another biggy is "O" instead of "0" or vise versa. .pa QUICK LIST OF MATHEMATICAL FUNCTIONS AERF: arc-error function APOLY: area of a polygon BETA: complete beta function BETAR: incomplete beta function BINML: binomial coefficient CHEBY: Chebyshev polynomial of N-th degree DWSN: Dawson's integral ERF: error function ERFC: complementary error function FACT: factorial FGQ0I: numerical integration from zero to infinity FGQ0ID: double precision form of FGI0I FGQ1: numerical integration over definite interval FGQ1D: double precision form of FGQ1 FGQ2: numerical integration over definite interval in 2 dimensions FGQ3: numerical integration over definite interval in 3 dimensions FPG: normal probability distribution in percent FRF: cubic error function FTC95: Student's T-distribution for 95% confidence GAMAL: natural log of the Gamma function for large values GAMMA: Gamma function GXNYM: numerical integration of X**N*Y**M*dXdY over polygonal region PLG0: Legendre polynomial of N-th degree PLG1: Legendre polynomial of N-th degree (first derivitive) PLG2: Legendre polynomial of N-th degree (second derivitive) PLG3: Legendre polynomial of N-th degree (third derivitive) PPLG0: Legendre polynomial of N-th degree in 2 dimensions RPQ: rational polynomial evaluation RPQD: double precision form of RPQ SEVAL: cubic spline evaluation QUICK LIST OF MATHEMATICAL SUBROUTINES BANDG: banded matrix solver BRYDN: modified Broyden's method for nonlinear simultaneous equations CUBIC: solve a cubic equation CUBICD: double precision form of CUBIC FMIN: find the minimum of a function within an interval 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 MDAY: convert Julian (annual) day to month/day/year MINV: matrix inversion MINVD: double precision form of MINV PROOT: polynomial root finder QUAD: solve a quadratic QUADD: double precision form of QUAD 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 SWMLR: stepwise multiple linear regression TDIAG: tridiagonal matrix solver TREND: compute trend of periodic function ZERO: find the zero (root) of a function within an interval .pa NAME: AERF PURPOSE: arc-error function TYPE: REAL*4 function (far external) SYNTAX: X=AERF(E) INPUT: E (REAL*4) the erf of something OUTPUT: X (REAL*4) the arc-erf NAME: APOLY PURPOSE: area of a polygon TYPE: REAL*4 function (far external) SYNTAX: A=APOLY(X,Y,N) INPUT: X,Y (REAL*4) an array of points (X,Y) OUTPUT: A (REAL*4) the area enclosed NOTE: points must be in the order you draw a "connect-the-dots" picture, the connection between the last point and the first point is assumed (e.g. for a triangular region N=3) NAME: BETA PURPOSE: complete beta function TYPE: REAL*4 function (far external) SYNTAX: B=BETA(X,Y) INPUT: X,Y (REAL*4) OUTPUT: B (REAL*4) NAME: BETAR PURPOSE: incomplete beta function TYPE: REAL*4 function (far external) SYNTAX: B=BETAR(X,Y,R) INPUT: X,Y,R (REAL*4) OUTPUT: B (REAL*4) NAME: BINML PURPOSE: binomial coefficient TYPE: REAL*4 function (far external) SYNTAX: B=BINML(N,I) INPUT: N,I (INTEGER*2) OUTPUT: B (REAL*4) NAME: CHEBY PURPOSE: Chebyshev polynomial of N-th degree TYPE: REAL*4 function (far external) SYNTAX: C=CHEBY(N,X) INPUT: X (REAL*4), N (INTEGER*2) OUTPUT: C (REAL*4) NAME: DWSN PURPOSE: Dawson's integral TYPE: REAL*4 function (far external) SYNTAX: D=DWSN(X) INPUT: X (REAL*4) OUTPUT: D (REAL*4) NAME: ERF PURPOSE: error function TYPE: REAL*4 function (far external) SYNTAX: E=ERF(X) INPUT: X (REAL*4) OUTPUT: E (REAL*4) NOTE: this is an unusually fast full machine precision algorithm NAME: ERFC PURPOSE: complementary error function TYPE: REAL*4 function (far external) SYNTAX: E=ERFC(X) INPUT: X (REAL*4) OUTPUT: E (REAL*4) NAME: FACT PURPOSE: factorial TYPE: REAL*4 function (far external) SYNTAX: F=FACT(N) INPUT: N (INTEGER*2) OUTPUT: F (REAL*4) NAME: FGQ0I (note the "0" (zero)) PURPOSE: numerical integration from zero to infinity TYPE: REAL*4 function (far external) SYNTAX: F=FGQ0I(SPARE) INPUT: SPARE (REAL*4) whatever you like (put 0. if no spare) FOFX (REAL*4 function) that YOU MUST SUPPLY and must be called by the sequence F=FOFX(X,SPARE) OUTPUT: F (REAL*4) NOTE: for double precision use FGQ0ID and name your function DOFX 20-point Gauss quadrature is used (96-point for double precision) NAME: FGQ1 PURPOSE: numerical integration over definite interval TYPE: REAL*4 function (far external) SYNTAX: F=FGQ1(A,B,SPARE) INPUT: A,B (REAL*4) interval over which to integrate FOFX SPARE (REAL*4) whatever you like (put 0. if no spare) FOFX (REAL*4 function) that YOU MUST SUPPLY and must be called by the sequence F=FOFX(X,SPARE) OUTPUT: F (REAL*4) NOTE: your function will be integrated from A to B for double precision use FGQ1D and name your function DOFX 20-point Gauss quadrature is used (96-point for double precision) NAME: FGQ2 PURPOSE: numerical integration over definite interval in 2-D TYPE: REAL*4 function (far external) SYNTAX: F=FGQ3(A,B,SPARE) INPUT: A,B (REAL*4) interval over which to integrate FOFX SPARE (REAL*4) whatever you like (put 0. if no spare) FOFX (REAL*4 function) that YOU MUST SUPPLY and must be called by the sequence F=FOFX(X,SPARE) FY1OFX (REAL*4 function) that YOU MUST SUPPLY and must be called by the sequence Y1=FY1OFX(X) FY2OFX (REAL*4 function) that YOU MUST SUPPLY and must be called by the sequence Y2=FY2OFX(X) OUTPUT: F (REAL*4) NOTE: your function will be integrated in X from A to B and Y from FY1OFX(X) to FY2OFX(X), 20*20=400 point Gauss quadrature is used NAME: FGQ3 PURPOSE: numerical integration over definite interval in 3-D TYPE: REAL*4 function (far external) SYNTAX: F=FGQ3(A,B,SPARE) INPUT: A,B (REAL*4) interval over which to integrate FOFX SPARE (REAL*4) whatever you like (put 0. if no spare) FOFX (REAL*4 function) that YOU MUST SUPPLY and must be called by the sequence F=FOFX(X,SPARE) FY1OFX (REAL*4 function) that YOU MUST SUPPLY and must be called by the sequence Y1=FY1OFX(X) FY2OFX (REAL*4 function) that YOU MUST SUPPLY and must be called by the sequence Y2=FY2OFX(X) FZ1OFXY (REAL*4 function) that YOU MUST SUPPLY and must be called by the sequence Z1=FZ1OFXY(X,Y) FZ2OFXY (REAL*4 function) that YOU MUST SUPPLY and must be called by the sequence Z2=FZ2OFXY(X,Y) OUTPUT: F (REAL*4) NOTE: your function will be integrated in X from A to B and Y from FY1OFX(X) to FY2OFX(X) and Z from FZ1OFXY(X,Y) to FZ2OFXY(X,Y), 20*20*20=8000 point Gauss quadrature is used NAME: FPG PURPOSE: normal probability distribution in percent TYPE: REAL*4 function (far external) SYNTAX: F=FPG(XBAR,SIGMA,X,DX) INPUT: XBAR (REAL*4) mean SIGMA (REAL*4) standard deviation X (REAL*4) independent variable DX (REAL*4) increment in X (if you want to know the probability of 0,5,10,15,20,25,...,95,100%, then X=0,5,10,... and DX=5) OUTPUT: F (REAL*4) probability in % NAME: FRF PURPOSE: cubic error function TYPE: REAL*4 function (far external) SYNTAX: F=FRF(X) INPUT: X (REAL*4) OUTPUT: F (REAL*4) NOTE: FRF is similar to ERF and also varies from -1 to 1 as X varies from -infinity to +infinity, it pops up in some problems like the error function (e.g. statistical thermodynamics) NAME: FTC95 PURPOSE: Student's T-distribution for 95% confidence TYPE: REAL*4 function (far external) SYNTAX: F=FTC95(N) INPUT: N (INTEGER*2) number of degrees of freedom OUTPUT: F (REAL*4) NOTE: "Student" is just the guy's hokey pseudonym - it doesn't mean anything NAME: GAMAL PURPOSE: natural log of the Gamma function for large values TYPE: REAL*4 function (far external) SYNTAX: G=GAMAL(X) INPUT: X (REAL*4) OUTPUT: G (REAL*4) NAME: GAMMA PURPOSE: Gamma function TYPE: REAL*4 function (far external) SYNTAX: G=GAMMA(X) INPUT: X (REAL*4) OUTPUT: G (REAL*4) NAME: GXNYM PURPOSE: numerical integration of X**N*Y**M*dXdY over polygonal region TYPE: REAL*4 function (far external) SYNTAX: G=GXNYM(X,Y,N,M,L) INPUT: X,Y (REAL*4) the points defining the boundary of a polygon N (INTEGER*2) the exponent on X (may be -+0) M (INTEGER*2) the exponent on Y (may be -+0) L (INTEGER*2) the number of points OUTPUT: G (REAL*4) NOTE: 8-point Gauss quadrature and Green's Lemma are used, integration will be exact (machine precision) for N+M<16. If N=M=0 then you get the area and you should use APOLY. To get the location of the centroid and the moment of inertia use the following XC=GXNYM(X,Y,1,0,L)/GXNYM(X,Y,0,0,L) YC=GXNYM(X,Y,0,1,L)/GXNYM(X,Y,0,0,L) AI=GXNYM(X,Y,1,1,L) You can get radii of gyration and other useful things like FEM coefficients from GXNYM. NAME: PLG0 PURPOSE: Legendre polynomial of N-th degree TYPE: REAL*4 function (far external) SYNTAX: P=PLG0(N,X) INPUT: N (INTEGER*2) degree X (REAL*4) OUTPUT: P (REAL*4) NAME: PLG1 PURPOSE: Legendre polynomial of N-th degree (first derivitive) TYPE: REAL*4 function (far external) SYNTAX: P=PLG1(N,X) INPUT: N (INTEGER*2) degree X (REAL*4) OUTPUT: P (REAL*4) NAME: PLG2 PURPOSE: Legendre polynomial of N-th degree (second derivitive) TYPE: REAL*4 function (far external) SYNTAX: P=PLG2(N,X) INPUT: N (INTEGER*2) degree X (REAL*4) OUTPUT: P (REAL*4) NAME: PLG3 PURPOSE: Legendre polynomial of N-th degree (third derivitive) TYPE: REAL*4 function (far external) SYNTAX: P=PLG3(N,X) INPUT: N (INTEGER*2) degree X (REAL*4) OUTPUT: P (REAL*4) NAME: PPLG0 PURPOSE: Legendre polynomial of N-th degree (2 dimensions) TYPE: REAL*4 function (far external) SYNTAX: P=PPLG0(N,X,Y) INPUT: N (INTEGER*2) degree X,Y (REAL*4) OUTPUT: P (REAL*4) NAME: RPQ PURPOSE: rational polynomial evaluation TYPE: REAL*4 function (far external) SYNTAX: R=RPQ(P,N,Q,M,X) INPUT: P (REAL*4) coefficients in numerator polynomial N (INTEGER*2) number of terms in P Q (REAL*4) coefficients in denominator polynomial M (INTEGER*2) number of terms in Q X (REAL*4) OUTPUT: R (REAL*4) NOTE: RPQ evaluates a rational polynomial as shown below using Horner's method R=(P1+P2*X+P3*X**2+P4*X**3+...)/(Q1+Q2*X+Q3*X**2+Q4*X**3...) for double precision use RPQD and make sure that you declare R,P,Q,RPQD all to be REAL*8 NAME: SEVAL PURPOSE: cubic spline evaluation TYPE: REAL*4 function (far external) SYNTAX: Y=SEVAL(XI,YI,N,C,X) INPUT: XI,YI (REAL*4) specified points N (INTEGER*2) number of points C (REAL*4) coefficients determined by first calling SPLNE X (REAL*4) point where you want to interpolate OUTPUT: Y (REAL*4) interpolated value NOTE: you need to first call SPLNE, you only need to do this once NAME: BANDG PURPOSE: banded matrix solver TYPE: subroutine (far external) SYNTAX: CALL BANDG(A,B,X,N,M,IER) INPUT: A,B (REAL*4) matrices N (INTEGER*2) number of equations M (INTEGER*2) bandwidth OUTPUT: X (REAL*4) solution vector IER (INTEGER*2) error indicator IER=0 indicates normal return IER=1 indicates M too small (M<3) IER=2 indicates m too big (M>N) IER=3 indicates M not odd IER=4 indicates [A] singular NOTE: The Gauss-Seidel method is used. The system of equations must be of the form [A][X]=[B] [A] should look something like the following (note zeroes)... if N=3 if N=5 if N=7 dimension A(3,N) A(5,N) A(7,N) 0 2 -2 0 0 3 -2 -1 0 0 0 6 -3 -2 -1 -1 2 -1 0 -2 5 -2 -1 0 0 -3 9 -3 -2 -1 -1 2 -1 -1 -2 6 -2 -1 0 -2 -3 11 -3 -2 -1 -1 2 -1 -1 -2 6 -2 -1 -1 -2 -3 12 -3 -2 -1 : : : : : : : : : : : : : : : -1 2 1 -1 -2 6 -2 -1 -1 -2 -3 12 -3 -2 -1 -1 2 1 -1 -2 6 -2 -1 -1 -2 -3 11 -3 -2 0 -1 2 1 -1 -2 5 -2 0 -1 -2 -3 9 -3 0 0 -2 2 0 -1 -2 3 0 0 -1 -2 -3 6 0 0 0 NAME: BRYDN PURPOSE: modified Broyden's method for nonlinear simultaneous equations TYPE: subroutine (far external) SYNTAX: CALL BRYDN(XMIN,XMAX,X,F,N,M,W,MW,PRINT,IER) INPUT: 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>=4N+3M+N*N+N*M PRINT (LOGICAL*2) set to .true. for printer output (also you need to open file 6 on the PC) RESID (a far external subroutine) which YOU MUST SUPPLY that may be called by the following CALL RESID(X,F) and returns "M" different Fs given "N" Xs 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: 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. 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. NAME: FMIN PURPOSE: find the minimum of a function within an interval TYPE: subroutine (far external) SYNTAX: CALL FMIN(A,B,X,F,SPARE) INPUT: A,B (REAL*4) the search interval SPARE (REAL*4) whatever you want to pass to your function FOFX (a far external REAL*4 function) that you must supply that may be called by the following F=FOFX(X,SPARE) OUTPUT: X (REAL*4) location of minimum F (REAL*4) function at minimum 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. 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. Sorry, there is no REAL*4 equivalent. eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee PROGRAM EXMPL C C THIS IS AN EXAMPLE OF HOW TO USE LLSLC C C IN THIS EXAMPLE A POWER SERIES (1,X,X**2,X**3,...) IS USED TO C APPROXIMATE THE FUNCTION F(T)=ERFC(T), THE APPROXIMATION C IS OVER T=0,1 AND THE CONSTRAINTS ARE THAT THE APPROXIMATION C FIT EXACTLY AT THE END POINTS C IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z) REAL*4 ERFC PARAMETER (NE=5,NC=2,NS=NE+NC,NF=NE-NC,NW=NS*(NS+3),M=20) DIMENSION X(NE),A(NE),WORK(NW),T(M),F(M) C OPEN(6,FILE='PRN') C C DEFINE THE FUNCTION F(T)=ERFC(T) C DO 100 I=1,M T(I)=DBLE(I-1)/DBLE(M-1) 100 F(I)=DBLE(ERFC(SNGL(T(I)))) C WRITE(6,1000) NE,NC,NF,NS 1000 FORMAT('1LINEAR LEAST-SQUARES WITH LINEAR CONSTRAINTS',//, &' NUMBER OF TERMS.................... ',I3,/, &' NUMBER OF CONSTRAINTS.............. ',I3,/, &' NUMBER OF FREE CONSTANTS........... ',I3,/, &' NUMBER OF SIMULTANEOUS EQUATIONS... ',I3) C C STEP#1 INITIALIZE (NOTE: IOPT=1) C CALL LLSLC(X,Y,A,NE,NC,WORK,NW,1,IER) IF(IER.NE.0) GO TO 900 C C STEP#2 FEED THE DATA TO LLSLC (NOTE: IOPT=2) C DO 210 I=1,M TI=T(I) FI=F(I) X(1)=1.D0 DO 200 J=2,NE 200 X(J)=X(J-1)*TI Y=FI C CALL LLSLC(X,Y,A,NE,NC,WORK,NW,2,IER) 210 IF(IER.NE.0) GO TO 900 C C STEP#3.1 ADD CONSTRAINT#1 (NOTE: IOPT=3) C NOTE: IF YOU HAVE NO CONSTRAINTS SET NC=0 AND SKIP SECTION 3 C I=1 X(I)=1.D0 TI=T(I) FI=F(I) DO 310 J=2,NE 310 X(J)=X(J-1)*TI Y=FI C CALL LLSLC(X,Y,A,NE,NC,WORK,NW,3,IER) IF(IER.NE.0) GO TO 900 C C STEP#3.2 ADD CONSTRAINT#2 (NOTE: IOPT=3) C I=M TI=T(I) FI=F(I) X(I)=1.D0 DO 320 J=2,NE 320 X(J)=X(J-1)*TI Y=FI C CALL LLSLC(X,Y,A,NE,NC,WORK,NW,3,IER) IF(IER.NE.0) GO TO 900 C C STEP#4 SOLVE FOR COEFFICIENTS (NOTE: IOPT=4) C CALL LLSLC(X,Y,A,NE,NC,WORK,NW,4,IER) IF(IER.NE.0) GO TO 900 C C LIST Q-FACTOR, COEFFICIENTS, AND SIGNIFICANCE LEVELS C WRITE(6,5000) Y 5000 FORMAT(/,' Q-FACTOR=',1PG12.5,' (Q<1 INDICATES STABLE, ', &'Q>1 INDICATES UNSTABLE)',//,' I A S') C DO 500 I=1,NE 500 WRITE(6,5100) I,A(I),X(I) 5100 FORMAT(1X,I3,2(1X,1PG12.5)) C C LIST AGREEMENT C WRITE(6,5200) 5200 FORMAT(/,' AGREEMENT',/, &' I X F Y E') C EMAX=0.D0 EAVE=0.D0 C DO 520 I=1,M X(1)=1.D0 Y=A(1)*X(1) DO 510 J=2,NE X(J)=X(J-1)*T(I) 510 Y=Y+A(J)*X(J) C E=Y-F(I) EAVE=EAVE+DABS(E) IF(DABS(E).GT.DABS(EMAX)) EMAX=E 520 WRITE(6,5300) I,T(I),F(I),Y,E 5300 FORMAT(1X,I3,4(1X,1PG12.5)) C C LIST AVERAGE AND MAXIMUM ERROR C EAVE=EAVE/DBLE(M) WRITE(6,5400) EAVE,EMAX 5400 FORMAT(30X,'AVERAGE ERROR=',1PG12.5,/, &30X,'MAXIMUM ERROR=',1PG12.5) GO TO 999 C 900 WRITE(6,9000) IER 9000 FORMAT(' ERROR CODE:',I5) C 999 CLOSE(6) STOP END eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee 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: 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: 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. NAME: RK4 PURPOSE: 4th order Runge-Kutta (solve ordinary differential equations) TYPE: subroutine (far external) SYNTAX: CALL RK4(X,DX,Y,DY,N,W,V,SPARE) 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) SPARE (REAL*4) whatever you want to pass to DYDX DYDX (far external) a subroutine that YOU MUST PROVIDE and must be called by the sequence CALL DYDX(X,Y,DY,N,SPARE) OUTPUT: Y (REAL*4) new values of dependent variable NOTE: this isn't too tricky; but I have provided an example that may help eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee PROGRAM EXMPL C C THIS IS AN EXAMPLE OF HOW TO USE RK4 C IN THIS EXAMPLE A 3-rd ORDER DIFFERENTIAL EQUATION WILL BE SOLVED C IMPLICIT INTEGER*2 (I-N),REAL*4 (A-H,O-Z) PARAMETER (N=3) DIMENSION Y(N),DY(N),W(N,4),V(N) C C DEFINE THE INITIAL CONDITIONS, THE STEP SIZE, AND THE NUMBER OF STEPS C DATA Y/0.,0.,0./ DATA DX,NSTEP/.1,100/ C OPEN(6,FILE='PRN') C WRITE(6,1000) 1000 FORMAT('1SOLUTION OF 3-RD ORDER EQUATION',//, &' X DDY/DXX DY/DX Y') C DO 100 ISTEP=1,NSTEP CALL RK4(X,DX,Y,DY,N,W,V,SPARE) 100 WRITE(6,1100) X,Y 1100 FORMAT(4(1X,1PG12.5)) C CLOSE(6) STOP END SUBROUTINE DYDX(X,Y,DY,N,SPARE) C C YOU MUST PROVIDE THIS SUBROUTINE! C IMPLICIT INTEGER*2 (I-N),REAL*4 (A-H,O-Z) DIMENSION Y(N),DY(N) C DY(1)=X*COS(X)-SIN(X) DY(2)=Y(1) DY(3)=Y(2) C RETURN END eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee NAME: RK6D PURPOSE: 6th order Runge-Kutta (solve ordinary differential equations) TYPE: subroutine (far external) SYNTAX: CALL RK6D(X,DX,Y,DY,N,W,V,SPARE) 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) SPARE (REAL*8) whatever you want to pass to DYDX DYDX (far external) a subroutine that YOU MUST PROVIDE and must be called by the sequence CALL DYDX(X,Y,DY,N,SPARE) 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: SWMLR PURPOSE: stepwise multiple linear regression TYPE: subroutine (far external) SYNTAX: CALL SWMLR(MA,MF,IER,ND,NF,R,X,Y,A,ATA,ATB,JPIVOT, &F,FTF,FTR,LF,SPARE) INPUT: MA (INTEGER*2) maximum number of terms in regression MF (INTEGER*2) maximum number of functions IER (INTEGER*2) error indicator (IER=0 is normal) ND (INTEGER*2) number of data points (X,Y) X,Y (REAL*8) data points to fit ATA (REAL*8) working space of dimension (MA,MA) ATB (REAL*8) working space of dimension (MA) ATB (INTEGER*2) working space of dimension (MA) F (REAL*8) working space of dimension (MA) FTF (REAL*8) working space of dimension (MA) FTR (REAL*8) working space of dimension (MA) OUTPUT: NF (INTEGER*2) number of functions used R (REAL*8) residual error after approximation A (REAL*8) coefficients LF (INTEGER*2) order in which the functions are used NOTE: This is a little tricky, so an example is provided below. Note that SWMLR will find the set of functions providing the most compact fit. This is different from the fit that has the least error or has the largest F-value and a LOT HARDER TO DETERMINE. eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee PROGRAM EXAMPL C C THIS IS AN EXAMPLE OF HOW TO USE SWMLR C IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z) PARAMETER (MA=8,MF=8,ND=20) REAL*4 ERF CHARACTER NAME*10 DIMENSION X(ND),Y(ND),A(MA),ATA(MA,MA),ATB(MA),JPIVOT(MA),F(MF), &FTF(MF),FTR(MF),LF(MF),R(ND),NAME(MF) DATA NAME/ &'1. ', &'X ', &'X*X ', &'X*X*X ', &'EXP(X) ', &'EXP(-X) ', &'X*EXP(X) ', &'X*EXP(-X) '/ C OPEN(6,FILE='PRN') WRITE(6,1000) 1000 FORMAT('1EXAMPLE OF STEPWISE LINEAR REGRESSION') C C DEFINE THE DATA TO BE FIT (THE ERROR FUNCTION) C DO 100 I=1,ND X(I)=DBLE(I-1)/DBLE(ND-1) 100 Y(I)=DBLE(ERF(SNGL(X(I)))) C CALL SWMLR(MA,MF,IER,ND,NF,R,X,Y,A,ATA,ATB,JPIVOT, &F,FTF,FTR,LF,SPARE) C WRITE(6,2000) 2000 FORMAT(/,' I J A NAME') C DO 200 I=1,NF 200 WRITE(6,2010) I,LF(I),A(I),NAME(LF(I)) 2010 FORMAT(2(2X,I1),2X,1PG12.5,2X,A10) C WRITE(6,2100) 2100 FORMAT(/,' X Y APPROX ERROR') C EAVE=0. EMAX=0. C DO 211 I=1,ND CALL SOFX(X(I),F,SPARE) C APPROX=0.D0 DO 210 J=1,NF 210 APPROX=APPROX+A(J)*F(LF(J)) C E=APPROX-Y(I) EAVE=EAVE+DABS(E) EMAX=DMAX1(EMAX,DABS(E)) 211 WRITE(6,2110) X(I),Y(I),APPROX,E 2110 FORMAT(4(1X,1PG12.5)) C EAVE=EAVE/DBLE(ND) WRITE(6,2200) EAVE,EMAX 2200 FORMAT(/, &' AVERAGE ERROR=',1PG12.5,/, &' MAXIMUM ERROR=',1PG12.5) C CLOSE(6) STOP END SUBROUTINE SOFX(X,F,SPARE) C C YOU MUST PROVIDE THIS SUBROUTINE! C IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z) PARAMETER (MF=8) DIMENSION F(MF) C F(1)=1. F(2)=X F(3)=X*X F(4)=X*X*X F(5)=EXP(X) F(6)=EXP(-X) F(7)=X*EXP(X) F(8)=X*EXP(-X) C RETURN END eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee 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 NAME: ZERO PURPOSE: find the zero (root) of a function within an interval TYPE: subroutine (far external) SYNTAX: CALL ZERO(A,B,X,IER,SPARE) INPUT: A,B (REAL*4) interval to look for root in SPARE (REAL*4) whatever you want to pass to your function FOFX (REAL*4 function) that YOU MUST SUPPLY and must be called by the sequence F=FOFX(X,SPARE) OUTPUT: X (REAL*4) the location of the root within (A<=X<=B) IER (INTEGER*2) error indicator (IER=0 is normal)