home *** CD-ROM | disk | FTP | other *** search
/ Phoenix CD 2.0 / Phoenix_CD.cdr / 01e / libry31a.zip / LIBRY4.DOC < prev    next >
Text File  |  1987-01-20  |  42KB  |  1,126 lines
  1. .pa
  2.                               MATHEMATICAL
  3.  
  4. These  mathematical  procedures provide a wide variety of functions that
  5. are useful  for  engineering  applications.   I  began  developing  this
  6. library after my first experience with SLIM (name changed to protect the
  7. guilty) back in 1976.  So I guess I have them to thank.  If you've  ever
  8. used SLIM perhaps you've had a similar experience...
  9.  
  10. Let's  say you want to solve what one would think to be a simple problem
  11. of simultaneous linear equations, and due to some  unfortunate  turn  of
  12. events,  you  are  given SLIM with which to do this.  You are confronted
  13. with no less than 60 subroutines from which to choose.  As it turns out,
  14. none  will  do  what  you  want.   Anyway,  you  narrow it  down to  the
  15. Smith-Jones-Simpson-Wilkerson-Feinstein-Gidley  method   for   partially
  16. symmetric   matrices  with  strange  sparceness  stored  in  convoluted-
  17. compacted form (not to be  confused  with  space-saving  form)  and  the
  18. Kleidowitz-Yonkers-Lakeville-Massey-Perkins  modification of the Orwell-
  19. Gauss-Newton-Schmeltz algorithm  for  ill-conditioned  Hilbert  matrices
  20. stored  in  space-saving  form  (not  to  be  confused  with convoluted-
  21. compacted form).  An obvious choice indeed (any fool with  6  Ph.D.s  in
  22. math could see that the latter is the preferred method).  After numerous
  23. attempts to link your  program  with  all  347  of  the  necessary  SLIM
  24. libraries  you  are  finally  ready  to run your program.  Of course, it
  25. bombs out because you were reading tuesday afternoon's manual  and  it's
  26. now wednesday morning.  Tough luck...
  27.  
  28. Actually,  SLIM is so bad that even I can't achieve hyperbole.  Language
  29. fails me.  The above could easily be a random selection from one of  the
  30. six volumes that are supposed to be a "user's guide."
  31.  
  32. Anyway, I have a completely different philosophy of programming.  I hope
  33. that it works well for you:  if there is a clearly preferred  method  of
  34. doing anything, then use it and discard the others.
  35.  
  36. An  excellent  example  of this is Gauss quadrature.  There are only two
  37. methods of numerical integration that will (at least in theory) converge
  38. given  an infinite order:  trapezoidal rule and Gauss quadrature.  Gauss
  39. quadrature is so far superior in accuracy  to  any  other  method  (e.g.
  40. Newton-Cotes) it's practically incredible.  Therefore, why have 50 other
  41. methods from which to choose?  An interesting note about this  is  that,
  42. contrary  to  some  rumors  10 subdivisions of 20-point Gauss quadrature
  43. (200 point evaluation) is not as accurate as 96-point  Gauss  quadrature
  44. with  less  than half the function evaluations.
  45.  
  46. Another example is Gauss-Jordan elimination.  SLIM has subroutines  that
  47. provide a "high accuracy solution" (If Andy Rooney understood this, he'd
  48. probably ask, "Who would want a low accuracy solution anyway?").  And of
  49. course,  there's  the  iterative  improvement type solutions.  What they
  50. don't tell you  is  that  the  residual  must  be  computed  in  greater
  51. precision than the matrix is reduced in (say single to double precision)
  52. which is OK;  but why not just use double precision all the way through?
  53. If you are already using double precision and that's not good enough for
  54. you...  well,  you're  just  out  of  luck.   If  your  matrix  is  ill-
  55. conditioned  (as  is  typically the case with Vandermondes - that's what
  56. you  get  when you're curve-fitting) about the best thing you can use is
  57. Gauss-Jordan elimination with full column pivoting.  It is  fast  enough
  58. and  I have used vector emulation throughout.  I have tried all sorts of
  59. things like QR decompositions, Givens rotations, and a  whole  bunch  of
  60. others;  but it all boils down to the precision of the math coprocessor.
  61. I  only  use  one algorithm for solving full matrices.  I have tested it
  62. using a  series  of  standard  Hilbert  matrices  against  SLIM's  "high
  63. accuracy solution" and found mine to be both faster and more accurate.
  64.  
  65. One last example is the solution of ordinary differential equations.  If
  66. you  read  numerical mathematics texts, they always give examples of the
  67. error in such-and-such a method as compared to the exact  solution.   If
  68. you  look  closely  at the microscopic print in the margin of the legend
  69. under the fly wing that got stuck on the copier, you will probably  find
  70. the  pithy  little statement "exact solution determined using 4-th order
  71. Runge-Kutta with 0.0000000001 step size." Runge-Kutta  is  self-starting
  72. (which  is  a  bigger pain in the neck than you might think if you use a
  73. method that isn't), convenient to use, and works very well.  So why  not
  74. use  it,  I  ask?  Of course there are certain (unknown before the fact)
  75. cases where other methods are much faster;  but  I  would  rather  spend
  76. less  time experimenting with zillions of algorithms and let the machine
  77. sweat a little more.  I've been through it once;  and  once  is  enough.
  78. Runge-Kutta is the one for me.  A word about automatic stepsize control.
  79. It has been my experience that this takes so long at each step that  you
  80. are  better  off  to  use a smaller step from the start.  I had two such
  81. subroutines that I spent a lot of time on;  but I purged them one day in
  82. a fit of frustration.
  83.  
  84. One  final  note:   WATCH  YOUR  VARIABLE  TYPES!    ESPECIALLY   DOUBLE
  85. PRECISION!   I've  seen  many  a programmer bewildered by a program that
  86. crashes because they have passed "0." instead of "0.D0" to a  subroutine
  87. expecting  a double precision value or forgotten to put the "REAL*8" out
  88. in front of a function.
  89.  
  90.           REAL*8 FUNCTION DOFX(D)
  91.  
  92. It is really best to use IMPLICIT statements like
  93.  
  94.           IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z)
  95.  
  96. than to individually type each variable.   You're  bound  to  miss  one.
  97.  
  98. Another biggy is "O" instead of "0" or vise versa.
  99. .pa
  100.                   QUICK LIST OF MATHEMATICAL FUNCTIONS
  101.  
  102. AERF:   arc-error function
  103. APOLY:  area of a polygon
  104. BETA:   complete beta function
  105. BETAR:  incomplete beta function
  106. BINML:  binomial coefficient
  107. CHEBY:  Chebyshev polynomial of N-th degree
  108. DWSN:   Dawson's integral
  109. ERF:    error function
  110. ERFC:   complementary error function
  111. FACT:   factorial
  112. FGQ0I:  numerical integration from zero to infinity
  113. FGQ0ID: double precision form of FGI0I
  114. FGQ1:   numerical integration over definite interval
  115. FGQ1D:  double precision form of FGQ1
  116. FGQ2:   numerical integration over definite interval in 2 dimensions
  117. FGQ3:   numerical integration over definite interval in 3 dimensions
  118. FPG:    normal probability distribution in percent
  119. FRF:    cubic error function
  120. FTC95:  Student's T-distribution for 95% confidence
  121. GAMAL:  natural log of the Gamma function for large values
  122. GAMMA:  Gamma function
  123. GXNYM:  numerical integration of X**N*Y**M*dXdY over polygonal region
  124. PLG0:   Legendre polynomial of N-th degree
  125. PLG1:   Legendre polynomial of N-th degree (first derivitive)
  126. PLG2:   Legendre polynomial of N-th degree (second derivitive)
  127. PLG3:   Legendre polynomial of N-th degree (third derivitive)
  128. PPLG0:  Legendre polynomial of N-th degree in 2 dimensions
  129. RPQ:    rational polynomial evaluation
  130. RPQD:   double precision form of RPQ
  131. SEVAL:  cubic spline evaluation
  132.  
  133.  
  134.                  QUICK LIST OF MATHEMATICAL SUBROUTINES
  135.  
  136. BANDG:  banded matrix solver
  137. BRYDN:  modified Broyden's method for nonlinear simultaneous equations
  138. CUBIC:  solve a cubic equation
  139. CUBICD: double precision form of CUBIC
  140. FMIN:   find the minimum of a function within an interval
  141. GAUSP:  Gauss-Jordan elimination of simultaneous equations
  142. GAUSPD: double precision form of GAUSP
  143. JLDAY:  convert month/day/year to Julian (annual) day
  144. LLSLC:  linear least-squares with linear constraints
  145. LNREG:  linear regression
  146. MDAY:   convert Julian (annual) day to month/day/year
  147. MINV:   matrix inversion
  148. MINVD:  double precision form of MINV
  149. PROOT:  polynomial root finder
  150. QUAD:   solve a quadratic
  151. QUADD:  double precision form of QUAD
  152. RK4:    4th order Runge-Kutta (solution of ordinary differential eqns.)
  153. RK6D:   6th order Runge-Kutta (solution of ordinary differential eqns.)
  154. SPLNE:  determine coefficients for a cubic spline
  155. SVD:    singular value decomposition using Householder transformations
  156. SWMLR:  stepwise multiple linear regression
  157. TDIAG:  tridiagonal matrix solver
  158. TREND:  compute trend of periodic function
  159. ZERO:   find the zero (root) of a function within an interval
  160. .pa
  161. NAME:     AERF
  162. PURPOSE:  arc-error function
  163. TYPE:     REAL*4 function (far external)
  164. SYNTAX:   X=AERF(E)
  165. INPUT:    E (REAL*4) the erf of something
  166. OUTPUT:   X (REAL*4) the arc-erf
  167.  
  168.  
  169. NAME:     APOLY
  170. PURPOSE:  area of a polygon
  171. TYPE:     REAL*4 function (far external)
  172. SYNTAX:   A=APOLY(X,Y,N)
  173. INPUT:    X,Y (REAL*4) an array of points (X,Y)
  174. OUTPUT:   A (REAL*4) the area enclosed
  175. NOTE:     points must be in the order you draw a "connect-the-dots"
  176.           picture, the connection between the last point and the first
  177.           point is assumed (e.g. for a triangular region N=3)
  178.  
  179.  
  180. NAME:     BETA
  181. PURPOSE:  complete beta function
  182. TYPE:     REAL*4 function (far external)
  183. SYNTAX:   B=BETA(X,Y)
  184. INPUT:    X,Y (REAL*4)
  185. OUTPUT:   B (REAL*4)
  186.  
  187.  
  188. NAME:     BETAR
  189. PURPOSE:  incomplete beta function
  190. TYPE:     REAL*4 function (far external)
  191. SYNTAX:   B=BETAR(X,Y,R)
  192. INPUT:    X,Y,R (REAL*4)
  193. OUTPUT:   B (REAL*4)
  194.  
  195.  
  196. NAME:     BINML
  197. PURPOSE:  binomial coefficient
  198. TYPE:     REAL*4 function (far external)
  199. SYNTAX:   B=BINML(N,I)
  200. INPUT:    N,I (INTEGER*2)
  201. OUTPUT:   B (REAL*4)
  202.  
  203.  
  204. NAME:     CHEBY
  205. PURPOSE:  Chebyshev polynomial of N-th degree
  206. TYPE:     REAL*4 function (far external)
  207. SYNTAX:   C=CHEBY(N,X)
  208. INPUT:    X (REAL*4), N (INTEGER*2)
  209. OUTPUT:   C (REAL*4)
  210.  
  211.  
  212. NAME:     DWSN
  213. PURPOSE:  Dawson's integral
  214. TYPE:     REAL*4 function (far external)
  215. SYNTAX:   D=DWSN(X)
  216. INPUT:    X (REAL*4)
  217. OUTPUT:   D (REAL*4)
  218.  
  219.  
  220. NAME:     ERF
  221. PURPOSE:  error function
  222. TYPE:     REAL*4 function (far external)
  223. SYNTAX:   E=ERF(X)
  224. INPUT:    X (REAL*4)
  225. OUTPUT:   E (REAL*4)
  226. NOTE:     this is an unusually fast full machine precision algorithm
  227.  
  228.  
  229. NAME:     ERFC
  230. PURPOSE:  complementary error function
  231. TYPE:     REAL*4 function (far external)
  232. SYNTAX:   E=ERFC(X)
  233. INPUT:    X (REAL*4)
  234. OUTPUT:   E (REAL*4)
  235.  
  236.  
  237. NAME:     FACT
  238. PURPOSE:  factorial
  239. TYPE:     REAL*4 function (far external)
  240. SYNTAX:   F=FACT(N)
  241. INPUT:    N (INTEGER*2)
  242. OUTPUT:   F (REAL*4)
  243.  
  244.  
  245. NAME:     FGQ0I (note the "0" (zero))
  246. PURPOSE:  numerical integration from zero to infinity
  247. TYPE:     REAL*4 function (far external)
  248. SYNTAX:   F=FGQ0I(SPARE)
  249. INPUT:    SPARE (REAL*4) whatever you like (put 0. if no spare)
  250.           FOFX (REAL*4 function) that YOU MUST SUPPLY and must be
  251.           called by the sequence F=FOFX(X,SPARE)
  252. OUTPUT:   F (REAL*4)
  253. NOTE:     for double precision use FGQ0ID and name your function DOFX
  254.           20-point Gauss quadrature is used (96-point for double
  255.           precision)
  256.  
  257.  
  258. NAME:     FGQ1
  259. PURPOSE:  numerical integration over definite interval
  260. TYPE:     REAL*4 function (far external)
  261. SYNTAX:   F=FGQ1(A,B,SPARE)
  262. INPUT:    A,B (REAL*4) interval over which to integrate FOFX
  263.           SPARE (REAL*4) whatever you like (put 0. if no spare)
  264.           FOFX (REAL*4 function) that YOU MUST SUPPLY and must be
  265.           called by the sequence F=FOFX(X,SPARE)
  266. OUTPUT:   F (REAL*4)
  267. NOTE:     your function will be integrated from A to B
  268.           for double precision use FGQ1D and name your function DOFX
  269.           20-point Gauss quadrature is used (96-point for double
  270.           precision)
  271.  
  272.  
  273. NAME:     FGQ2
  274. PURPOSE:  numerical integration over definite interval in 2-D
  275. TYPE:     REAL*4 function (far external)
  276. SYNTAX:   F=FGQ3(A,B,SPARE)
  277. INPUT:    A,B (REAL*4) interval over which to integrate FOFX
  278.           SPARE (REAL*4) whatever you like (put 0. if no spare)
  279.           FOFX (REAL*4 function) that YOU MUST SUPPLY and must be
  280.           called by the sequence F=FOFX(X,SPARE)
  281.           FY1OFX (REAL*4 function) that YOU MUST SUPPLY and must be
  282.           called by the sequence Y1=FY1OFX(X)
  283.           FY2OFX (REAL*4 function) that YOU MUST SUPPLY and must be
  284.           called by the sequence Y2=FY2OFX(X)
  285. OUTPUT:   F (REAL*4)
  286. NOTE:     your function will be integrated in X from A to B and Y
  287.           from FY1OFX(X) to FY2OFX(X), 20*20=400 point Gauss quadrature
  288.           is used
  289.  
  290.  
  291. NAME:     FGQ3
  292. PURPOSE:  numerical integration over definite interval in 3-D
  293. TYPE:     REAL*4 function (far external)
  294. SYNTAX:   F=FGQ3(A,B,SPARE)
  295. INPUT:    A,B (REAL*4) interval over which to integrate FOFX
  296.           SPARE (REAL*4) whatever you like (put 0. if no spare)
  297.           FOFX (REAL*4 function) that YOU MUST SUPPLY and must be
  298.           called by the sequence F=FOFX(X,SPARE)
  299.           FY1OFX (REAL*4 function) that YOU MUST SUPPLY and must be
  300.           called by the sequence Y1=FY1OFX(X)
  301.           FY2OFX (REAL*4 function) that YOU MUST SUPPLY and must be
  302.           called by the sequence Y2=FY2OFX(X)
  303.           FZ1OFXY (REAL*4 function) that YOU MUST SUPPLY and must be
  304.           called by the sequence Z1=FZ1OFXY(X,Y)
  305.           FZ2OFXY (REAL*4 function) that YOU MUST SUPPLY and must be
  306.           called by the sequence Z2=FZ2OFXY(X,Y)
  307. OUTPUT:   F (REAL*4)
  308. NOTE:     your function will be integrated in X from A to B and Y
  309.           from FY1OFX(X) to FY2OFX(X) and Z from FZ1OFXY(X,Y) to
  310.           FZ2OFXY(X,Y), 20*20*20=8000 point Gauss quadrature
  311.           is used
  312.  
  313.  
  314. NAME:     FPG
  315. PURPOSE:  normal probability distribution in percent
  316. TYPE:     REAL*4 function (far external)
  317. SYNTAX:   F=FPG(XBAR,SIGMA,X,DX)
  318. INPUT:    XBAR (REAL*4) mean
  319.           SIGMA (REAL*4) standard deviation
  320.           X (REAL*4) independent variable
  321.           DX (REAL*4) increment in X (if you want to know the
  322.           probability of 0,5,10,15,20,25,...,95,100%, then X=0,5,10,...
  323.           and DX=5)
  324. OUTPUT:   F (REAL*4) probability in %
  325.  
  326.  
  327. NAME:     FRF
  328. PURPOSE:  cubic error function
  329. TYPE:     REAL*4 function (far external)
  330. SYNTAX:   F=FRF(X)
  331. INPUT:    X (REAL*4)
  332. OUTPUT:   F (REAL*4)
  333. NOTE:     FRF is similar to ERF and also varies from -1 to 1 as X varies
  334.           from -infinity to +infinity, it pops up in some problems like
  335.           the error function (e.g. statistical thermodynamics)
  336.  
  337.  
  338. NAME:     FTC95
  339. PURPOSE:  Student's T-distribution for 95% confidence
  340. TYPE:     REAL*4 function (far external)
  341. SYNTAX:   F=FTC95(N)
  342. INPUT:    N (INTEGER*2) number of degrees of freedom
  343. OUTPUT:   F (REAL*4)
  344. NOTE:     "Student" is just the guy's hokey pseudonym - it doesn't mean
  345.           anything
  346.  
  347.  
  348. NAME:     GAMAL
  349. PURPOSE:  natural log of the Gamma function for large values
  350. TYPE:     REAL*4 function (far external)
  351. SYNTAX:   G=GAMAL(X)
  352. INPUT:    X (REAL*4)
  353. OUTPUT:   G (REAL*4)
  354.  
  355.  
  356. NAME:     GAMMA
  357. PURPOSE:  Gamma function
  358. TYPE:     REAL*4 function (far external)
  359. SYNTAX:   G=GAMMA(X)
  360. INPUT:    X (REAL*4)
  361. OUTPUT:   G (REAL*4)
  362.  
  363.  
  364. NAME:     GXNYM
  365. PURPOSE:  numerical integration of X**N*Y**M*dXdY over polygonal region
  366. TYPE:     REAL*4 function (far external)
  367. SYNTAX:   G=GXNYM(X,Y,N,M,L)
  368. INPUT:    X,Y (REAL*4) the points defining the boundary of a polygon
  369.           N (INTEGER*2) the exponent on X (may be -+0)
  370.           M (INTEGER*2) the exponent on Y (may be -+0)
  371.           L (INTEGER*2) the number of points
  372. OUTPUT:   G (REAL*4)
  373. NOTE:     8-point Gauss quadrature and Green's Lemma are used,
  374.           integration will be exact (machine precision) for N+M<16.
  375.           If N=M=0 then you get the area and you should use APOLY.
  376.           To get the location of the centroid and the moment of inertia
  377.           use the following
  378.  
  379.           XC=GXNYM(X,Y,1,0,L)/GXNYM(X,Y,0,0,L)
  380.           YC=GXNYM(X,Y,0,1,L)/GXNYM(X,Y,0,0,L)
  381.           AI=GXNYM(X,Y,1,1,L)
  382.  
  383.           You can get radii of gyration and other useful things like FEM
  384.           coefficients from GXNYM.
  385.  
  386.  
  387. NAME:     PLG0
  388. PURPOSE:  Legendre polynomial of N-th degree
  389. TYPE:     REAL*4 function (far external)
  390. SYNTAX:   P=PLG0(N,X)
  391. INPUT:    N (INTEGER*2) degree
  392.           X (REAL*4)
  393. OUTPUT:   P (REAL*4)
  394.  
  395.  
  396. NAME:     PLG1
  397. PURPOSE:  Legendre polynomial of N-th degree (first derivitive)
  398. TYPE:     REAL*4 function (far external)
  399. SYNTAX:   P=PLG1(N,X)
  400. INPUT:    N (INTEGER*2) degree
  401.           X (REAL*4)
  402. OUTPUT:   P (REAL*4)
  403.  
  404.  
  405. NAME:     PLG2
  406. PURPOSE:  Legendre polynomial of N-th degree (second derivitive)
  407. TYPE:     REAL*4 function (far external)
  408. SYNTAX:   P=PLG2(N,X)
  409. INPUT:    N (INTEGER*2) degree
  410.           X (REAL*4)
  411. OUTPUT:   P (REAL*4)
  412.  
  413.  
  414. NAME:     PLG3
  415. PURPOSE:  Legendre polynomial of N-th degree (third derivitive)
  416. TYPE:     REAL*4 function (far external)
  417. SYNTAX:   P=PLG3(N,X)
  418. INPUT:    N (INTEGER*2) degree
  419.           X (REAL*4)
  420. OUTPUT:   P (REAL*4)
  421.  
  422.  
  423. NAME:     PPLG0
  424. PURPOSE:  Legendre polynomial of N-th degree (2 dimensions)
  425. TYPE:     REAL*4 function (far external)
  426. SYNTAX:   P=PPLG0(N,X,Y)
  427. INPUT:    N (INTEGER*2) degree
  428.           X,Y (REAL*4)
  429. OUTPUT:   P (REAL*4)
  430.  
  431.  
  432. NAME:     RPQ
  433. PURPOSE:  rational polynomial evaluation
  434. TYPE:     REAL*4 function (far external)
  435. SYNTAX:   R=RPQ(P,N,Q,M,X)
  436. INPUT:    P (REAL*4) coefficients in numerator polynomial
  437.           N (INTEGER*2) number of terms in P
  438.           Q (REAL*4) coefficients in denominator polynomial
  439.           M (INTEGER*2) number of terms in Q
  440.           X (REAL*4)
  441. OUTPUT:   R (REAL*4)
  442. NOTE:     RPQ evaluates a rational polynomial as shown below using
  443.           Horner's method
  444.  
  445.           R=(P1+P2*X+P3*X**2+P4*X**3+...)/(Q1+Q2*X+Q3*X**2+Q4*X**3...)
  446.  
  447.           for double precision use RPQD and make sure that you declare
  448.           R,P,Q,RPQD all to be REAL*8
  449.  
  450.  
  451. NAME:     SEVAL
  452. PURPOSE:  cubic spline evaluation
  453. TYPE:     REAL*4 function (far external)
  454. SYNTAX:   Y=SEVAL(XI,YI,N,C,X)
  455. INPUT:    XI,YI (REAL*4) specified points
  456.           N (INTEGER*2) number of points
  457.           C (REAL*4) coefficients determined by first calling SPLNE
  458.           X (REAL*4) point where you want to interpolate
  459. OUTPUT:   Y (REAL*4) interpolated value
  460. NOTE:     you need to first call SPLNE, you only need to do this once
  461.  
  462.  
  463. NAME:     BANDG
  464. PURPOSE:  banded matrix solver
  465. TYPE:     subroutine (far external)
  466. SYNTAX:   CALL BANDG(A,B,X,N,M,IER)
  467. INPUT:    A,B (REAL*4) matrices
  468.           N (INTEGER*2) number of equations
  469.           M (INTEGER*2) bandwidth
  470. OUTPUT:   X (REAL*4) solution vector
  471.           IER (INTEGER*2) error indicator
  472.                IER=0  indicates normal return
  473.                IER=1  indicates M too small (M<3)
  474.                IER=2  indicates m too big (M>N)
  475.                IER=3  indicates M not odd
  476.                IER=4  indicates [A] singular
  477. NOTE:     The Gauss-Seidel method is used.
  478.           The system of equations must be of the form
  479.  
  480.                [A][X]=[B]
  481.  
  482.           [A] should look something like the following (note zeroes)...
  483.  
  484.             if N=3           if N=5               if N=7
  485. dimension   A(3,N)           A(5,N)               A(7,N)
  486.            0  2 -2       0  0  3 -2 -1      0  0  0  6 -3 -2 -1
  487.           -1  2 -1       0 -2  5 -2 -1      0  0 -3  9 -3 -2 -1
  488.           -1  2 -1      -1 -2  6 -2 -1      0 -2 -3 11 -3 -2 -1
  489.           -1  2 -1      -1 -2  6 -2 -1     -1 -2 -3 12 -3 -2 -1
  490.            :  :  :       :  :  :  :  :      :  :  :  :  :  :  :
  491.           -1  2  1      -1 -2  6 -2 -1     -1 -2 -3 12 -3 -2 -1
  492.           -1  2  1      -1 -2  6 -2 -1     -1 -2 -3 11 -3 -2  0
  493.           -1  2  1      -1 -2  5 -2  0     -1 -2 -3  9 -3  0  0
  494.           -2  2  0      -1 -2  3  0  0     -1 -2 -3  6  0  0  0
  495.  
  496.  
  497. NAME:     BRYDN
  498. PURPOSE:  modified Broyden's method for nonlinear simultaneous equations
  499. TYPE:     subroutine (far external)
  500. SYNTAX:   CALL BRYDN(XMIN,XMAX,X,F,N,M,W,MW,PRINT,IER)
  501. INPUT:    XMIN,XMAX (REAL*4) arrays, upper and lower limits on X
  502.           N (INTEGER*2) number of elements in X
  503.           M (INTEGER*2) number of residuals
  504.           W (REAL*4) working space of dimension MW
  505.           MW (INTEGER*2) dimension of working space MW>=4N+3M+N*N+N*M
  506.           PRINT (LOGICAL*2) set to .true. for printer output (also you
  507.           need to open file 6 on the PC)
  508.           RESID (a far external subroutine) which YOU MUST SUPPLY that
  509.           may be called by the following
  510.  
  511.                CALL RESID(X,F)
  512.  
  513.           and returns "M" different Fs given "N" Xs
  514. OUTPUT:   X (REAL*4) solution vector
  515.           F (REAL*4) final residual
  516.           IER (INTEGER*2) error indicator
  517.                IER=0  indicates normal return
  518.                IER=1  indicates N<2
  519.                IER=2  indicates M<N
  520.                IER=3  indicates XMIN>XMAX conflict
  521.                IER=4  indicates insufficient work space (increase MW)
  522. NOTE:     This is a FAR more difficult problem than most people even
  523.           imagine!  You can compare it to finding the lowest spot on
  524.           the face of the Earth given the rules of the game "Battle
  525.           Ship" where all you can do is call out coordinates and your
  526.           opponent answers with the elevation.  You may find a rut or a
  527.           ditch somewhere; but it's highly unlikely that you will find
  528.           death valley, let alone some trench in the Pacific. Now extend
  529.           this analogy to many-dimensional space... get the point?  So
  530.           don't get too upset if BRYDN can't find the solution you want.
  531.           I've tried all sorts of algorithms and this seems to work the
  532.           best for general problems.
  533.  
  534.           If you have 6 equations and 6 unknowns then M=N=6.  If you
  535.           have 12 equations and 6 unknowns then M=12 and N=6.
  536.  
  537.           If all of your equations are linear, by all means don't use
  538.           BRYDN, you want LLSLC.
  539.  
  540.           If you have only one unknown then use FMIN.
  541.  
  542.  
  543. NAME:     CUBIC
  544. PURPOSE:  solve a cubic equation
  545. TYPE:     subroutine (far external)
  546. SYNTAX:   CALL SUBROUTINE CUBIC(A1,A2,A3,A4,R1,C1,R2,C2,R3,C3)
  547. INPUT:    A1,A2,A3,A4 (REAL*4) coefficients in cubic equation
  548. OUTPUT:   (R1,C1),(R2,C2),(R3,C3) (REAL*4) roots
  549. NOTE:     the equation has the form A1*X**3+A2*X**2+A3*X+A4=0
  550.           for double precision use CUBICD and don't forget to declare
  551.           A1,A2,A3,A4,R1,C1,R2,C2,R3,C3 all to be REAL*8 and if you
  552.           use constants, don't forget the 1.D0 or whatever.
  553.  
  554.  
  555. NAME:     FMIN
  556. PURPOSE:  find the minimum of a function within an interval
  557. TYPE:     subroutine (far external)
  558. SYNTAX:   CALL FMIN(A,B,X,F,SPARE)
  559. INPUT:    A,B (REAL*4) the search interval
  560.           SPARE (REAL*4) whatever you want to pass to your function
  561.           FOFX (a far external REAL*4 function) that you must supply
  562.           that may be called by the following
  563.  
  564.                F=FOFX(X,SPARE)
  565.  
  566. OUTPUT:   X (REAL*4) location of minimum
  567.           F (REAL*4) function at minimum
  568.  
  569.  
  570. NAME:     GAUSP
  571. PURPOSE:  Gauss-Jordan elimination of simultaneous equations
  572. TYPE:     subroutine (far external)
  573. SYNTAX:   CALL GAUSP(A,B,X,JPIVOT,N,D,C,IER)
  574. INPUT:    A,B (REAL*4) arrays containing the equations to be solved
  575.           having dimensions (N,N) and (N) respectively
  576.           JPIVOT (INTEGER*2) working space of dimension "N"
  577.           N (INTEGER*2) the number of equations
  578. OUTPUT:   X (REAL*4) solution vector of dimension N
  579.           D (REAL*4) ALOG10(ABS(DETERMINANT)) log of the determinant
  580.           C (REAL*4) ALOG10(AMAX1(ABS(PIVOT))/AMIN1(ABS(PIVOT))) a
  581.           measure of the gradedness or condition of the matrix
  582.           IER (INTEGER*2) error indicator (IER=0 is normal)
  583. NOTE:     For double precision use GAUSPD and don't forget to declare
  584.           A,B,X,D,C all to be REAL*8.
  585.           Note that A and B will be destroyed upon return.
  586.           Note that A is in row-order as indicated below
  587.  
  588.                A(1,1)*X(1)+A(1,2)*X(2)+...A(1,N)*X(N)=B(1)
  589.                A(2,1)*X(1)+A(2,2)*X(2)+...A(2,N)*X(N)=B(2)
  590.                A(3,1)*X(1)+A(3,2)*X(2)+...A(3,N)*X(N)=B(3)
  591.                ..........................................
  592.                A(N,1)*X(1)+A(N,2)*X(2)+...A(N,N)*X(N)=B(N)
  593.  
  594.           When you really have to be careful about this is with DATA
  595.           statements as FORTRAN fills arrays in column-order.
  596.  
  597.           Full column pivoting is used.
  598.           Vector emulation is used throughout to give maximum speed.
  599.  
  600.  
  601. NAME:     JLDAY
  602. PURPOSE:  convert month/day/year to Julian (annual) day
  603. TYPE:     subroutine (far external)
  604. SYNTAX:   CALL JLDAY(MONTH,IDAY,IYEAR,JDAY)
  605. INPUT:    MONTH,IDAY,IYEAR (INTEGER*2)
  606. OUTPUT:   JDAY (INTEGER*2)
  607. NOTE:     I know that "Julian" day is a misnomer; but that's what
  608.           everyone else calls it.
  609.  
  610.  
  611. NAME:     LLSLC
  612. PURPOSE:  linear least-squares with linear constraints
  613. TYPE:     subroutine (far external)
  614. SYNTAX:   CALL LLSLC(X,Y,A,NE,NC,W,NW,IOPT,IER)
  615. INPUT:    X,Y,A,W (REAL*8) see example below
  616.           NE (INTEGER*2) number of variables (X or A)
  617.           NC (INTEGER*2) number of constraints (may be zero)
  618.           NW (INTEGER*2) working space, NW>=(NE+NC)*(NE+NC+3))
  619.           IOPT (INTEGER*2) option, see example below
  620. OUTPUT:   A (REAL*8)
  621.           IER (INTEGER*2) error indicator
  622.                IER=0   normal
  623.                IER=1   dimension error
  624.                IER=2   singular matrix - no inverse
  625.                IER=3   too few constants+constraints (minimum=2)
  626.                IER=4   too few equations (minimum=1)
  627.                IER=5   too few constraints (minimum=0)
  628.                IER=6   too many constraints (maximum=ne-1)
  629.                IER=7   IOPT not defined
  630.                IER=8   matrices not initialized
  631.                IER=9   too few points for least-squares
  632.                IER=10  constraints added before least-squares points
  633.                IER=11  too many constraints were added
  634.                IER=12  too few constraints were specified
  635.                IER=13  working space too small (increase NW)
  636. NOTE:     This is a little tricky, so I have included an example.
  637.           Sorry, there is no REAL*4 equivalent.
  638.  
  639.  
  640. eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
  641.       PROGRAM EXMPL
  642. C
  643. C  THIS IS AN EXAMPLE OF HOW TO USE LLSLC
  644. C
  645. C  IN THIS EXAMPLE A POWER SERIES (1,X,X**2,X**3,...) IS USED TO
  646. C  APPROXIMATE THE FUNCTION F(T)=ERFC(T),  THE APPROXIMATION
  647. C  IS OVER T=0,1 AND THE CONSTRAINTS ARE THAT THE APPROXIMATION
  648. C  FIT EXACTLY AT THE END POINTS
  649. C
  650.       IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z)
  651.       REAL*4 ERFC
  652.       PARAMETER (NE=5,NC=2,NS=NE+NC,NF=NE-NC,NW=NS*(NS+3),M=20)
  653.       DIMENSION X(NE),A(NE),WORK(NW),T(M),F(M)
  654. C
  655.       OPEN(6,FILE='PRN')
  656. C
  657. C  DEFINE THE FUNCTION F(T)=ERFC(T)
  658. C
  659.       DO 100 I=1,M
  660.       T(I)=DBLE(I-1)/DBLE(M-1)
  661.   100 F(I)=DBLE(ERFC(SNGL(T(I))))
  662. C
  663.       WRITE(6,1000) NE,NC,NF,NS
  664.  1000 FORMAT('1LINEAR LEAST-SQUARES WITH LINEAR CONSTRAINTS',//,
  665.      &' NUMBER OF TERMS.................... ',I3,/,
  666.      &' NUMBER OF CONSTRAINTS.............. ',I3,/,
  667.      &' NUMBER OF FREE CONSTANTS........... ',I3,/,
  668.      &' NUMBER OF SIMULTANEOUS EQUATIONS... ',I3)
  669. C
  670. C  STEP#1  INITIALIZE  (NOTE: IOPT=1)
  671. C
  672.       CALL LLSLC(X,Y,A,NE,NC,WORK,NW,1,IER)
  673.       IF(IER.NE.0) GO TO 900
  674. C
  675. C  STEP#2  FEED THE DATA TO LLSLC  (NOTE: IOPT=2)
  676. C
  677.       DO 210 I=1,M
  678.       TI=T(I)
  679.       FI=F(I)
  680.       X(1)=1.D0
  681.       DO 200 J=2,NE
  682.   200 X(J)=X(J-1)*TI
  683.       Y=FI
  684. C
  685.       CALL LLSLC(X,Y,A,NE,NC,WORK,NW,2,IER)
  686.   210 IF(IER.NE.0) GO TO 900
  687. C
  688. C  STEP#3.1  ADD CONSTRAINT#1  (NOTE: IOPT=3)
  689. C  NOTE: IF YOU HAVE NO CONSTRAINTS SET NC=0 AND SKIP SECTION 3
  690. C
  691.       I=1
  692.       X(I)=1.D0
  693.       TI=T(I)
  694.       FI=F(I)
  695.       DO 310 J=2,NE
  696.   310 X(J)=X(J-1)*TI
  697.       Y=FI
  698. C
  699.       CALL LLSLC(X,Y,A,NE,NC,WORK,NW,3,IER)
  700.       IF(IER.NE.0) GO TO 900
  701. C
  702. C  STEP#3.2  ADD CONSTRAINT#2  (NOTE: IOPT=3)
  703. C
  704.       I=M
  705.       TI=T(I)
  706.       FI=F(I)
  707.       X(I)=1.D0
  708.       DO 320 J=2,NE
  709.   320 X(J)=X(J-1)*TI
  710.       Y=FI
  711. C
  712.       CALL LLSLC(X,Y,A,NE,NC,WORK,NW,3,IER)
  713.       IF(IER.NE.0) GO TO 900
  714. C
  715. C  STEP#4  SOLVE FOR COEFFICIENTS  (NOTE: IOPT=4)
  716. C
  717.       CALL LLSLC(X,Y,A,NE,NC,WORK,NW,4,IER)
  718.       IF(IER.NE.0) GO TO 900
  719. C
  720. C  LIST Q-FACTOR, COEFFICIENTS, AND SIGNIFICANCE LEVELS
  721. C
  722.       WRITE(6,5000) Y
  723.  5000 FORMAT(/,' Q-FACTOR=',1PG12.5,'  (Q<1 INDICATES STABLE, ',
  724.      &'Q>1 INDICATES UNSTABLE)',//,'   I       A            S')
  725. C
  726.       DO 500 I=1,NE
  727.   500 WRITE(6,5100) I,A(I),X(I)
  728.  5100 FORMAT(1X,I3,2(1X,1PG12.5))
  729. C
  730. C  LIST AGREEMENT
  731. C
  732.       WRITE(6,5200)
  733.  5200 FORMAT(/,' AGREEMENT',/,
  734.      &'   I       X            F            Y            E')
  735. C
  736.       EMAX=0.D0
  737.       EAVE=0.D0
  738. C
  739.       DO 520 I=1,M
  740.       X(1)=1.D0
  741.       Y=A(1)*X(1)
  742.       DO 510 J=2,NE
  743.       X(J)=X(J-1)*T(I)
  744.   510 Y=Y+A(J)*X(J)
  745. C
  746.       E=Y-F(I)
  747.       EAVE=EAVE+DABS(E)
  748.       IF(DABS(E).GT.DABS(EMAX)) EMAX=E
  749.   520 WRITE(6,5300) I,T(I),F(I),Y,E
  750.  5300 FORMAT(1X,I3,4(1X,1PG12.5))
  751. C
  752. C  LIST AVERAGE AND MAXIMUM ERROR
  753. C
  754.       EAVE=EAVE/DBLE(M)
  755.       WRITE(6,5400) EAVE,EMAX
  756.  5400 FORMAT(30X,'AVERAGE ERROR=',1PG12.5,/,
  757.      &30X,'MAXIMUM ERROR=',1PG12.5)
  758.       GO TO 999
  759. C
  760.   900 WRITE(6,9000) IER
  761.  9000 FORMAT(' ERROR CODE:',I5)
  762. C
  763.   999 CLOSE(6)
  764.       STOP
  765.       END
  766. eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
  767.  
  768.  
  769. NAME:     LNREG
  770. PURPOSE:  linear regression
  771. TYPE:     subroutine (far external)
  772. SYNTAX:   CALL LNREG(X,Y,N,A,B,R,IER)
  773. INPUT:    X,Y (REAL*4) data points
  774.           N (INTEGER*2) number of points
  775. OUTPUT:   A,B (REAL*4) slope, intercept
  776.           R (REAL*4) regression coefficient in %
  777.           IER (INTEGER*2) error indicator (IER=0 is normal)
  778. NOTE:     LNREG only fits a straight line of the form Y=A*X+B to a set
  779.           of points.  If you want something more elaborate use LLSLC.
  780.  
  781.  
  782. NAME:     MDAY
  783. PURPOSE:  convert Julian (annual) day to month/day/year
  784. TYPE:     subroutine (far external)
  785. SYNTAX:   CALL MDAY(JDAY,IYEAR,MONTH,IDAY)
  786. INPUT:    JDAY,IYEAR (INTEGER*2)
  787. OUTPUT:   MONTH,IDAY (INTEGER*2)
  788. NOTE:     I know that "Julian" day is a misnomer; but that's what
  789.           everyone else calls it.
  790.  
  791.  
  792. NAME:     MINV
  793. PURPOSE:  matrix inversion
  794. TYPE:     subroutine (far external)
  795. SYNTAX:   CALL MINV(A,IPIVOT,JPIVOT,N,D,C,IER)
  796. INPUT:    A (REAL*4) array containing the matrix to be inverted
  797.           IPIVOT,JPIVOT (INTEGER*2) working spaces of dimension "N"
  798.           N (INTEGER*2) rank (the number of equations)
  799. OUTPUT:   A (REAL*4) inverted matrix
  800.           D (REAL*4) ALOG10(ABS(DETERMINANT)) log of the determinant
  801.           C (REAL*4) ALOG10(AMAX1(ABS(PIVOT))/AMIN1(ABS(PIVOT))) a
  802.           measure of the gradedness or condition of the matrix
  803.           IER (INTEGER*2) error indicator (IER=0 is normal)
  804. NOTE:     For double precision use MINVD and don't forget to declare
  805.           A,D,C all to be REAL*8.
  806.           Note that A is inverted "in place."
  807.           Note that A is in row-order (see GAUSP).
  808.  
  809.  
  810. NAME:     PROOT
  811. PURPOSE:  polynomial root finder
  812. TYPE:     subroutine (far external)
  813. SYNTAX:   CALL PROOT(P,R,C,N,IER)
  814. INPUT:    P (REAL*8) array containing the polynomial coefficients
  815.           N (INTEGER*2) number of terms in P
  816. OUTPUT:   R (REAL*8) array containing the real part of the roots
  817.           C (REAL*8) array containing the complex part of the roots
  818.           IER (INTEGER*2) error indicator
  819.                IER=0 indicates no error
  820.                IER=1 indicates input error
  821.                IER=2 indicates failure to converge
  822. NOTE:     Bairstow's method is used.
  823.           The polynomial must be of the form
  824.  
  825.                P1+P2*X+P3*X**2+P4*X**3+....=0
  826.  
  827.           Also note that this is a VERY difficult problem for large N,
  828.           so don't be too upset if you get complex roots when you
  829.           thought you should be getting real roots.  That's life in the
  830.           real world.
  831.  
  832.  
  833. NAME:     QUAD
  834. PURPOSE:  solve a quadratic
  835. TYPE:     subroutine (far external)
  836. SYNTAX:   CALL SUBROUTINE QUAD(A1,A2,A3,R1,C1,R2,C2)
  837. INPUT:    A1,A2,A3 (REAL*4) coefficients in quadratic equation
  838. OUTPUT:   (R1,C1),(R2,C2) (REAL*4) roots
  839. NOTE:     the equation has the form A1*X**2+A2*X+A3=0
  840.           for double precision use QUADD and don't forget to declare
  841.           A1,A2,A3,R1,C1,R2,C2 all to be REAL*8 and if you use
  842.           constants, don't forget the 1.D0 or whatever.
  843.  
  844.  
  845. NAME:     RK4
  846. PURPOSE:  4th order Runge-Kutta (solve ordinary differential equations)
  847. TYPE:     subroutine (far external)
  848. SYNTAX:   CALL RK4(X,DX,Y,DY,N,W,V,SPARE)
  849. INPUT:    X (REAL*4) independent variable
  850.           DX (REAL*4) increment in X (step size)
  851.           Y (REAL*4) array, dependent variables at X
  852.           DY (REAL*4) array, change in Y
  853.           N (INTEGER*2) number of dependent variables (Y)
  854.           W (REAL*4) working space of dimension (N,4)
  855.           V (REAL*4) working space of dimension (N)
  856.           SPARE (REAL*4) whatever you want to pass to DYDX
  857.           DYDX (far external) a subroutine that YOU MUST PROVIDE and
  858.           must be called by the sequence
  859.  
  860.                CALL DYDX(X,Y,DY,N,SPARE)
  861.  
  862. OUTPUT:   Y (REAL*4) new values of dependent variable
  863. NOTE:     this isn't too tricky; but I have provided an example that may
  864.           help
  865.  
  866.  
  867. eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
  868.       PROGRAM EXMPL
  869. C
  870. C  THIS IS AN EXAMPLE OF HOW TO USE RK4
  871. C  IN THIS EXAMPLE A 3-rd ORDER DIFFERENTIAL EQUATION WILL BE SOLVED
  872. C
  873.       IMPLICIT INTEGER*2 (I-N),REAL*4 (A-H,O-Z)
  874.       PARAMETER (N=3)
  875.       DIMENSION Y(N),DY(N),W(N,4),V(N)
  876. C
  877. C  DEFINE THE INITIAL CONDITIONS, THE STEP SIZE, AND THE NUMBER OF STEPS
  878. C
  879.       DATA Y/0.,0.,0./
  880.       DATA DX,NSTEP/.1,100/
  881. C
  882.       OPEN(6,FILE='PRN')
  883. C
  884.       WRITE(6,1000)
  885.  1000 FORMAT('1SOLUTION OF 3-RD ORDER EQUATION',//,
  886.      &'       X         DDY/DXX       DY/DX          Y')
  887. C
  888.       DO 100 ISTEP=1,NSTEP
  889.       CALL RK4(X,DX,Y,DY,N,W,V,SPARE)
  890.   100 WRITE(6,1100) X,Y
  891.  1100 FORMAT(4(1X,1PG12.5))
  892. C
  893.       CLOSE(6)
  894.       STOP
  895.       END
  896.       SUBROUTINE DYDX(X,Y,DY,N,SPARE)
  897. C
  898. C  YOU MUST PROVIDE THIS SUBROUTINE!
  899. C
  900.       IMPLICIT INTEGER*2 (I-N),REAL*4 (A-H,O-Z)
  901.       DIMENSION Y(N),DY(N)
  902. C
  903.       DY(1)=X*COS(X)-SIN(X)
  904.       DY(2)=Y(1)
  905.       DY(3)=Y(2)
  906. C
  907.       RETURN
  908.       END
  909. eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
  910.  
  911.  
  912. NAME:     RK6D
  913. PURPOSE:  6th order Runge-Kutta (solve ordinary differential equations)
  914. TYPE:     subroutine (far external)
  915. SYNTAX:   CALL RK6D(X,DX,Y,DY,N,W,V,SPARE)
  916. INPUT:    X (REAL*8) independent variable
  917.           DX (REAL*8) increment in X (step size)
  918.           Y (REAL*8) array, dependent variables at X
  919.           DY (REAL*8) array, change in Y
  920.           N (INTEGER*2) number of dependent variables (Y)
  921.           W (REAL*8) working space of dimension (N,8) <-- note that the
  922.           "8" here is not a misprint - it takes 8 steps to get 6th order
  923.           even though it only takes 4 to get 4-order Runge-Kutta
  924.           V (REAL*4) working space of dimension (N)
  925.           SPARE (REAL*8) whatever you want to pass to DYDX
  926.           DYDX (far external) a subroutine that YOU MUST PROVIDE and
  927.           must be called by the sequence
  928.  
  929.                CALL DYDX(X,Y,DY,N,SPARE)
  930.  
  931. OUTPUT:   Y (REAL*8) new values of dependent variable
  932. NOTE:     This is just like RK4 except for the "8" above and that you
  933.           need to declare everything REAL*8 instead of REAL*4.
  934.  
  935.  
  936. NAME:     SPLNE
  937. PURPOSE:  determine coefficients for a cubic spline
  938. TYPE:     subroutine (far external)
  939. SYNTAX:   CALL SPLNE(XI,YI,N,C)
  940. INPUT:    XI,YI (REAL*4) points to be fit
  941.           N (INTEGER*2) number of points
  942. OUTPUT:   C (REAL*4) coefficient array of dimension (3,N)
  943. NOTE:     you need to call this once, then use SEVAL to evaluate spline
  944.  
  945.  
  946. NAME:     SVD
  947. PURPOSE:  singular value decomposition using Householder transformations
  948. TYPE:     subroutine (far external)
  949. SYNTAX:   CALL SVD(A,N,M,S,U,V,W,IER)
  950. INPUT:    A (REAL*8) the matrix (N,M) to be decomposed - row order
  951.           N (INTEGER*2) number of rows in A
  952.           M (INTEGER*2) number of columns in A
  953.           W (REAL*8) working space of dimension MAX0(N,M)
  954. OUTPUT:   S (REAL*8) singular values
  955.           U (REAL*8) left eigenvectors (M,N)
  956.           V (REAL*8) right eigenvectors (N,M)
  957.           IER (INTEGER*2) error indicator (IER=0 is normal)
  958.  
  959.  
  960. NAME:     SWMLR
  961. PURPOSE:  stepwise multiple linear regression
  962. TYPE:     subroutine (far external)
  963. SYNTAX:   CALL SWMLR(MA,MF,IER,ND,NF,R,X,Y,A,ATA,ATB,JPIVOT,
  964.          &F,FTF,FTR,LF,SPARE)
  965. INPUT:    MA (INTEGER*2) maximum number of terms in regression
  966.           MF (INTEGER*2) maximum number of functions
  967.           IER (INTEGER*2) error indicator (IER=0 is normal)
  968.           ND (INTEGER*2) number of data points (X,Y)
  969.           X,Y (REAL*8) data points to fit
  970.           ATA (REAL*8) working space of dimension (MA,MA)
  971.           ATB (REAL*8) working space of dimension (MA)
  972.           ATB (INTEGER*2) working space of dimension (MA)
  973.           F (REAL*8) working space of dimension (MA)
  974.           FTF (REAL*8) working space of dimension (MA)
  975.           FTR (REAL*8) working space of dimension (MA)
  976. OUTPUT:   NF (INTEGER*2) number of functions used
  977.           R (REAL*8) residual error after approximation
  978.           A (REAL*8) coefficients
  979.           LF (INTEGER*2) order in which the functions are used
  980. NOTE:     This is a little tricky, so an example is provided below.
  981.           Note that SWMLR will find the set of functions providing the
  982.           most compact fit.  This is different from the fit that has
  983.           the least error or has the largest F-value and a LOT HARDER
  984.           TO DETERMINE.
  985.  
  986.  
  987. eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
  988.       PROGRAM EXAMPL
  989. C
  990. C  THIS IS AN EXAMPLE OF HOW TO USE SWMLR
  991. C
  992.       IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z)
  993.       PARAMETER (MA=8,MF=8,ND=20)
  994.       REAL*4 ERF
  995.       CHARACTER NAME*10
  996.       DIMENSION X(ND),Y(ND),A(MA),ATA(MA,MA),ATB(MA),JPIVOT(MA),F(MF),
  997.      &FTF(MF),FTR(MF),LF(MF),R(ND),NAME(MF)
  998.       DATA NAME/
  999.      &'1.        ',
  1000.      &'X         ',
  1001.      &'X*X       ',
  1002.      &'X*X*X     ',
  1003.      &'EXP(X)    ',
  1004.      &'EXP(-X)   ',
  1005.      &'X*EXP(X)  ',
  1006.      &'X*EXP(-X) '/
  1007. C
  1008.       OPEN(6,FILE='PRN')
  1009.       WRITE(6,1000)
  1010.  1000 FORMAT('1EXAMPLE OF STEPWISE LINEAR REGRESSION')
  1011. C
  1012. C  DEFINE THE DATA TO BE FIT (THE ERROR FUNCTION)
  1013. C
  1014.       DO 100 I=1,ND
  1015.       X(I)=DBLE(I-1)/DBLE(ND-1)
  1016.   100 Y(I)=DBLE(ERF(SNGL(X(I))))
  1017. C
  1018.       CALL SWMLR(MA,MF,IER,ND,NF,R,X,Y,A,ATA,ATB,JPIVOT,
  1019.      &F,FTF,FTR,LF,SPARE)
  1020. C
  1021.       WRITE(6,2000)
  1022.  2000 FORMAT(/,'  I  J       A        NAME')
  1023. C
  1024.       DO 200 I=1,NF
  1025.   200 WRITE(6,2010) I,LF(I),A(I),NAME(LF(I))
  1026.  2010 FORMAT(2(2X,I1),2X,1PG12.5,2X,A10)
  1027. C
  1028.       WRITE(6,2100)
  1029.  2100 FORMAT(/,'      X            Y         APPROX        ERROR')
  1030. C
  1031.       EAVE=0.
  1032.       EMAX=0.
  1033. C
  1034.       DO 211 I=1,ND
  1035.       CALL SOFX(X(I),F,SPARE)
  1036. C
  1037.       APPROX=0.D0
  1038.       DO 210 J=1,NF
  1039.   210 APPROX=APPROX+A(J)*F(LF(J))
  1040. C
  1041.       E=APPROX-Y(I)
  1042.       EAVE=EAVE+DABS(E)
  1043.       EMAX=DMAX1(EMAX,DABS(E))
  1044.   211 WRITE(6,2110) X(I),Y(I),APPROX,E
  1045.  2110 FORMAT(4(1X,1PG12.5))
  1046. C
  1047.       EAVE=EAVE/DBLE(ND)
  1048.       WRITE(6,2200) EAVE,EMAX
  1049.  2200 FORMAT(/,
  1050.      &' AVERAGE ERROR=',1PG12.5,/,
  1051.      &' MAXIMUM ERROR=',1PG12.5)
  1052. C
  1053.       CLOSE(6)
  1054.       STOP
  1055.       END
  1056.       SUBROUTINE SOFX(X,F,SPARE)
  1057. C
  1058. C  YOU MUST PROVIDE THIS SUBROUTINE!
  1059. C
  1060.       IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z)
  1061.       PARAMETER (MF=8)
  1062.       DIMENSION F(MF)
  1063. C
  1064.       F(1)=1.
  1065.       F(2)=X
  1066.       F(3)=X*X
  1067.       F(4)=X*X*X
  1068.       F(5)=EXP(X)
  1069.       F(6)=EXP(-X)
  1070.       F(7)=X*EXP(X)
  1071.       F(8)=X*EXP(-X)
  1072. C
  1073.       RETURN
  1074.       END
  1075. eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
  1076.  
  1077.  
  1078. NAME:     TDIAG
  1079. PURPOSE:  tridiagonal matrix solver
  1080. TYPE:     subroutine (far external)
  1081. SYNTAX:   CALL TDIAG(D,B,X,W,N,IER)
  1082. INPUT:    D (REAL*4) matrix of dimension (N,3)
  1083.           B (REAL*4) right-hand-side of dimension (N)
  1084.           W (REAL*4) working space of dimension (N,3)
  1085.           N (INTEGER*2) number of equations
  1086. OUTPUT:   X (REAL*4) solution vector of dimension (N)
  1087.           IER (INTEGER*2) error indicator (IER=0 is normal)
  1088. NOTE:     the tridiagonal system of equations is defined by  [D][X]=[B]
  1089.           [D] should look something like the following
  1090.  
  1091.                  0  2 -2
  1092.                 -1  2 -1
  1093.                 -1  2 -1
  1094.                  .  .  .
  1095.                 -2  2  0
  1096.  
  1097.           For double precision use TDIAGD and don't forget to declare
  1098.           D,B,X,W all to be of type REAL*8.
  1099.  
  1100.  
  1101. NAME:     TREND
  1102. PURPOSE:  compute trend of periodic function
  1103. TYPE:     subroutine (far external)
  1104. SYNTAX:   CALL TREND(X,Y,W,NX,N,C,NC)
  1105. INPUT:    X,Y (REAL*4) points to be fit
  1106.           W (REAL*4) working space of dimension (NX)
  1107.           NX (INTEGER*2) number of points
  1108.           N (INTEGER*2) number of known points (N<NX)
  1109.           NC (INTEGER*2) number of terms in expansion (try 6)
  1110. OUTPUT:   Y (REAL*4) upon return, the N+1,N+2,N+3,...NX values of Y will
  1111.           be trended
  1112.           C (REAL*4) coefficients in trend series
  1113. NOTE:     Chebyshev approximation is used
  1114.  
  1115.  
  1116. NAME:     ZERO
  1117. PURPOSE:  find the zero (root) of a function within an interval
  1118. TYPE:     subroutine (far external)
  1119. SYNTAX:   CALL ZERO(A,B,X,IER,SPARE)
  1120. INPUT:    A,B (REAL*4) interval to look for root in
  1121.           SPARE (REAL*4) whatever you want to pass to your function
  1122.           FOFX (REAL*4 function) that YOU MUST SUPPLY and must be
  1123.           called by the sequence F=FOFX(X,SPARE)
  1124. OUTPUT:   X (REAL*4) the location of the root within (A<=X<=B)
  1125.           IER (INTEGER*2) error indicator (IER=0 is normal)
  1126.