# for the C, double precision version say 'send dloess from a'
# to unbundle, sh this file (in an empty directory)
echo README 1>&2
-Fortran Routines for Locally-Weighted Regression 31 Aug 1990
-
-William S. Cleveland
-Eric Grosse
-
-Locally-weighted regression, or loess, is a procedure for estimating a
-regression surface by a multivariate smoothing procedure: fitting a
-linear or quadratic function of the independent variables in a moving
-fashion that is analogous to how a moving average is computed for a
-time series. Compared to classical approaches - fitting global
-parametric functions - loess substantially increases the domain of
-surfaces that can be estimated without distortion. Also, a pleasant
-fact about loess is that analogues of the statistical procedures used
-in parametric function fitting - for example, ANOVA and t intervals -
-involve statistics whose distributions are well approximated by
-familiar distributions.
-
-The file sample.f is an example main program which reads one line
- n p alpha totaldeg fcell
-and then n lines of p+1 numbers (x,y). The paramater alpha controls
-the amount of smoothing; a typical starting value might be .5;
-totaldeg is the degree polynomials to use, either 1 or 2; fcell is
-a stopping criterion, typically set to alpha/5 or alpha/10.
-For more detail on calling sequences of the Fortran routines, see
-the file "internal".
-
-Both single precision and double precision versions are available;
-double is preferred when you use quadratic fitting on 32-bit machines
-like IEEE standard, IBM, and VAX.
-
-For the version distributed by electronic mail from netlib, subroutines
-from linpack have been omitted. If those are not already on your system,
- mail netlib@netlib.bell-labs.com
- send r1mach snrm2 ssvdc sqrdc sdot sqrsl isamax from linpack core
-if you are using the single precision version; for double precision,
- send d1mach dnrm2 dsvdc dqrdc ddot dqrsl idamax from linpack core.
-When installing, don't forget to uncomment the appropriate DATA statements
-in r1mach or d1mach, as described by the comments in those functions.
-
-Bug reports are appreciated. Send electronic mail to
- ehg@netlib.bell-labs.com
-including the words "this is not spam" in the Subject line
-or send paper mail to
- Eric Grosse
- Bell Labs 2T-502
- Murray Hill NJ 07974
-Remember that this is experimental software distributed free of charge
-and comes with no warranty! Exercise professional caution.
-
-Happy Smoothing!
//GO.SYSIN DD README
echo internal 1>&2
-***************************************************************
-* LOESS smoothing scattered data in one or more variables *
-* documentation of Fortran routines *
-* Cleveland, Devlin, Grosse, Shyu *
-***************************************************************
-
-1. The typical program would call lowesd, set tolerances in iv,v if
- desired, then call lowesb and lowese.
-2. To save the k-d tree, call lowesd, lowesb and then losave; subsequent
- programs would call lohead, set liv and lv, then call loread and lowese.
-3. For statistics, get diagL and then call lowesa or get the full hat
- matrix and call lowesc. Robustness iterations can take advantage of
- lowesw and lowesp.
-
-lowesd(106,iv,liv,lv,v,d,n,f,tdeg,nvmax,setLf) setup workspace
-lowesf(x,y,w,iv,liv,lv,v,m,z,L,hat,s) slow smooth at z
-lowesb(x,y,w,diagL,infl,iv,liv,lv,v) build k-d tree
-lowesr(y,iv,liv,lv,v) rebuild with new data values
- (does not change y)
-lowese(iv,liv,lv,v,m,z, s) evaluate smooth at z
-lowesl(iv,liv,lv,v,m,z, L) explicit hat matrix,
- which maps from y to z
-lofort(iunit,iv,liv,lv,v) save k-d tree as Fortran
-losave(iunit,iv,liv,lv,v) save k-d tree in file
-lohead(iunit,d,vc,nc,nv) read d,vc,nc,nv from file
- liv = 50+(vc+3)*nc determine space
- lv = 50+(2*d+1)*nv+nc requirements
-loread(iunit,d,vc,nc,nv,iv,liv,lv,v) finish reading k-d tree,
- ready for lowese
-lowesa(trL,n,d,tau,nsing, del1,del2) approximate delta
-lowesc(n,L,LL, trL,del1,del2) exact delta
-lowesp(n,y,yhat,w,rw, pi,ytilde) pseudo-values
-lowesw(res,n, rw,pi) robustness weights
-
-=== arguments ===
-d number of independent variables [integer] (called "p" elsewhere)
-del1,del2 delta1, delta2
-diagL diagonal of hat matrix, only set if infl=.true. (n)
-f fraction of points to use in local smooth (called "alpha" elsewhere)
-fc don't refine cells with less than fc*n points; ordinarily=.05
-hat is hat matrix desired? [integer]
- 0 = none
- 1 = diagonal only
- 2 = full matrix
-infl is diagonal of hat matrix desired? [logical]
-iunit Fortran unit number for i/o
-iv workspace (liv)
-L hat matrix (m,n) [real]
- in lowesf, only computed if hat nonzero; if hat=1 only size (n)
-LL workspace (n,n)
-liv 50+(2^d+4)*nvmax+2*n
- if setLf, add nf*nvmax
-lv 50+(3*d+3)*nvmax+n+(tau0+2)*nf
- if setLf, add (d+1)*nf*nvmax
-m number of points to smooth at; ordinarily=n
-n number of observations
-nf min(n,floor(n*f))
-nsing if 0, print warning in lowesa when trL<tau; typically nsing=iv(30)
-nvmax limit on number of vertices for kd-tree; e.g. max(200,n)
-pi workspace (n) [integer]
-res residual yhat-y (n)
-rw robustness weights (n)
-s smoothed values at z (m)
-setLf in lowesb, save matrix factorizations [logical]
- (needed for lowesr and lowesl)
-tau dimension of local model = iv(DIM);
- =d+1 for linear, (d+2)(d+1)/2 for quadratic
- reduced if dropping squares
-tau0 =d+1 for linear, (d+2)(d+1)/2 for quadratic
-tdeg polynomials to fit; 0=constants, 1=linear, 2=quadratics
-trL trace L = sum diagL
-v workspace (lv)
-w weights (n) local regression: min sum wi * (f(xi)-yi)^2
-x sample locations (n,d)
-y observations (n)
-yhat smoothed y (n)
-ytilde pseudo y (n)
-z locations where smooth is desired (m,d)
-
-If using the double precision version, [real] above should be understood
-as Fortran "double precision".
-
-The first argument to lowesd is a version number, updated when calling
-sequences change.
-
-If you peek inside the fortran, you will quickly notice that it
-was machine generated; the typeset original (in the language "pine")
-is much easier to read.
-
-=== iv indices ===
-1 INFO return code (not currently used)
-2 D number of independent variables
-3 N number of observations
-4 VC 2^d (number of vertices of a cell)
-5 NC number of k-d cells
-6 NV number of k-d vertices
-7 A1 starting index in iv of a
-8 C1 starting index in iv of c
-9 HI1 starting index in iv of hi
-10 LO1 starting index in iv of lo
-11 V1 starting index in v of vertices
-12 XI1 starting index in v of cut values
-13 VV1 starting index in v of vertex values
-14 NVMAX maximum allowed value of nv
-15 WORK1 starting index in v of workspace
-16 WORK2 starting index in v of workspace
-17 NCMAX maximum allowed value of nc
-18 WORK3 starting index in v of workspace
-19 NF floor(n*f) (number of points used as neighborhood)
-20 KERNEL 1=tricube, 2=unif
-21 KIND 1=k-d,cubic blend, (not implemented:2=quadtree,3=triangulation)
-22 PI1 starting index in iv of tree permutation
-23 VH starting index in iv of vhit
-24 VV2 starting index in v of work vval used in trL computation
-25 LQ starting index in iv of Lq
-26 WORK4 starting index in v of workspace
-27 PSI1 starting index in iv of workspace permutation
-28 SEQ sequence number, to check if routines called out of order
- takes on values:
- 171 after lowesd
- 172 after lowesf
- 173 after lowesb
-29 DIM dimension of local regression
- 1 constant
- d+1 linear (default)
- (d+2)(d+1)/2 quadratic
- Modified by ehg127 if cdeg<tdeg.
-30 SING number of times singular tolerance was met in l2fit, l2tr
-31 PRINT verbose output?
-32 DEG total degree (of polynomial for local model)
-33 NDIST dd = variables 1:dd enter into distance calculation
-34 LF starting index in v of Lf
-35..40 reserved for future use
-41..49 CDEG componentwise degree
-iv(A1) a coordinate of cut; 0 for leaf (nc)
-iv(C1) c pointers to corners (index into vertex array v) (vc,nc)
-iv(HI1) hi right subcell (nc)
-iv(LO1) lo left subcell (nc)
- Leaf cell j encloses points x(pi(i),), lo(j)<=i<=hi(j).
- Also, iv(C1),...,iv(PI1-1) is used as workspace (t) by l2fit
-------------------------eval only needs workspace up to here
-iv(PI1) pi permutation of 1:n for listing points in cells
-iv(VH) vhit cell whose subdivision creates vertex (nv)
- 0 if vertex is corner of original bounding box.
-iv(LQ) Lq active point indices for block of Lf (nvmax*nf)
-iv(PSI1) psi workspace permutation of 1:n for sorting distances
-
-=== v indices ===
-1 F fraction of n to be used as neighborhood. See also iv(19).
-2 FCELL no refinement if #points <= fcell * n
- default .05
-3 FDIAM no refinement if diameter is fdiam * overall bounding box
- default 0; Warning: reset to 0 by ehg142 when nsteps>0.
-4 RCOND reciprocal condition number
-... 49 reserved for future use
-iv(V1) v vertices (nv,d)
-iv(VV1) vval vertex values (0:d,nv)
-iv(XI1) xi cut values (nc)
-------------------------eval only needs workspace up to here
-iv(WORK1) workspace (n) l2fit:dist
-iv(WORK2) workspace (nf) l2fit:eta
-iv(WORK3) workspace (dim,nf) l2fit:X
-iv(VV2) vval2 workspace ((d+1)*nv) pseudo-vval for trL computation
-iv(LF) Lf hat matrix (data to vertex) ((d+1)*nvmax*nf)
-iv(WORK4) workspace (nf) l2fit:w
-
-Internal routine names have been hidden as follows:
-ehg106 select q-th smallest by partial sorting
-ehg124 rbuild
-ehg125 cpvert
-ehg126 bbox
-ehg127 l2fit,l2tr computational kernel
-ehg128 eval
-ehg129 spread
-ehg131 lowesb after workspace expansion
-ehg133 lowese after workspace expansion
-ehg134 abort by calling S Recover function
-ehg136 l2fit with hat matrix L
-ehg137 vleaf
-ehg138 descend
-ehg139 l2tr
-ehg140(w,i,j) w(i)=j used when w is declared real, but should store an int
-ehg141 delta1,2 from trL
-ehg142 robust iteration
-ehg144 now called lowesc
-ehg152 like ehg142, but for lowesf
-ehg167 kernel for losave
-ehg168 kernel for loread
-ehg169 compute derived k-d tree information
-ehg170 generate Fortran
-ehg176,ehg177,ehg178,ehg179,ehg180,ehg181 loeval for delta
-ehg182 ehgdie(i)
-ehg183 warning(message,i,n,inc)
-ehg184 warning(message,x,n,inc)
-ehg190 now called lowesa, with slight change in calling sequence
-ehg191 lowesl after workspace expansion
-ehg192 lowesr after workspace expansion
-ehg196(tau,d,f,trl) trL approximation
-ehg197 for deg=1,2
-m9rwt now called lowesw
-pseudo now called lowesp
//GO.SYSIN DD internal
echo changes 1>&2
-CHANGES PLANNED SOMEDAY
-1) more vertices in k-d tree for dimension > 2, to get continuity.
-2) triangulation based method.
-----------------------
-
-19 Nov 1987 workspace not big enough for degree=2
-
-22 Jan 1988 switched from depth first to breadth first tree build
-
-14 Mar 1988 lostt.3 extra space needed if (method mod 1000 = 0),
- not the documented (method/1000=0)
-
-28 Apr 1988 l2tr.g vval2 needed to be initialized to 0
-
- galaxy smooth needs double precision on vax
-
-26 May 1988 bbox.g add 10% margin to allow limited extrapolation
-
-6 June 1988 loess/lostt.f trL wasn't set if method/1000==0
-
-10 June 1988 losave, loread
-
- v(RCOND) 1 / max condition number
-
-12 June 1988 lofort
-
-21 June 1988 additional workspace for explicit L
-
-27 June 1988 workspace checking in lowesf was slightly pessimistic
-
-30 June 1988 Changed default fdiam to 0.
- Added warning messages for memory limits and pseudoinverse.
-
-4 Aug 1988 bbox.g changed margin from 10% to 0.5%.
-
-24 Aug 1988 loser documentation should have specified workspace
- of size ...+m*n, not ...+m**2.
-
-Sep 1988
- loess-based approximations of delta1,2.
- pseudo-values, so statistics are available with robustness iterations.
- reorganize error messages to better fit into S.
- sample driver program.
- somewhat shorter code generated by ehg170.
-
-20 Dec 1988
- workspace in loser
-
-27 Jan 1989
- workspace checking in lostt was a bit pessimistic.
-
-3 Feb 1989
- l2fit, l2tr: error message should contain sqrt(rho)
-
-18 Dec 1989
- ehg141, ehg179-ehg181: new delta approximations
-
-24 Jan 1990
- master copy moved from Sun3/180 to SGI 4D/240S
- (no intentional changes)
-
-25 Jan 1990
- (many routines touched; ehg127 added) cleaned up computational
- kernel, added provision for only first dd<=d variables to enter
- the distance calculation ("conditionally parametric variables"),
- added independent bounds on total and componentwise degree for
- local polynomial model, made extrapolation warning message print
- a bit more detail.
-
-14 Mar 1990
- added setLf argument to lowesd; added lowesr, lowesl for resmoothing.
-
--------------------------------------------------------
-Converting to the new version of loess
-5 April 1990
-
-Over the past few months, a number of changes have been made to the
-loess package, to provide more control over the local model, to allow
-conditionally parametric variables, and to return exact statistical
-quantities for the blending method. Unlike earlier internal
-algorithmic improvements, this round of changes added some extra
-arguments in the Fortran calling sequences. The purpose of this note
-is to assist in converting programs that called the old version.
-
-An explicit argument setLf has been added to lowesd(), since it affects
-the partitioning of the workspace. To help protect against inadvertent
-version mismatches, the version number that lowesd() checks has also
-been changed. The componentwise degree and the specification of
-conditionally nonparametric variables can be changed from the default
-by modifying iv(CDEG) and iv(NDIST).
-
-The influence matrix L for blending is now explicitly available by
-calling a new subroutine lowesl(), but this loses the speed
-advantage of blending. A faster, sometimes equivalent method is
-to use the influence matrix that carries data values to coefficients
-at the vertices of the k-d tree. This information is saved in iv(iv(Lq))
-and v(iv(Lf)), for the afficionado.
-
-The new subroutine lowesr() takes advantage of Lq and Lf to allow rapid
-resmoothing for applications when only y, not x, is subject to change.
--------------------------------------------------------
-
-7 May 1990
- new delta approximations.
- added prior weights to input format for sample driver.
-
-29 May 1990
- loess,lostt,loser,pseudo moved from Fortran to S.
-
-11 Jul 1990
- column equilibration, so pseudoinverse is needed less often.
-
-27 May 1991
-lowesd version 105; increased nvmax,ncmax to max(200,n).
-l2fit added ihat=1 (diagL only).
-ehg133,lowese removed unused arguments dist,eta.
-ehg190,ehg141 changed name to lowesa, slight change to calling sequence.
-ehg144 changed name to lowesc
-m9rwt changed name to lowesw
-pseudo changed name to lowesp
-
-22 Jul 1991 IMPORTANT BUG FIX!
-ehg131 vval2 should be dimensioned 0:d, not 0:8
-
-26 Jul 1991
-lowesd change calling sequence to provide tighter memory allocation
-diff old/internal new/internal
-< lowesd(105,iv,liv,lv,v,d,n,f,tdeg,setLf) setup workspace
-> lowesd(106,iv,liv,lv,v,d,n,f,tdeg,nvmax,setLf) setup workspace
-< liv 50+(2^d+6)*max(200,n)
-< if setLf, add nf*max(200,n)
-< lv 50+(3*d+4)*max(200,n)+(tau+2)*nf
-< if setLf, add (d+1)*nf*max(200,n)
-> liv 50+(2^d+6)*nvmax
-> if setLf, add nf*nvmax
-> lv 50+(3*d+4)*nvmax+(tau+2)*nf
-> if setLf, add (d+1)*nf*nvmax
-> nvmax limit on number of vertices for kd-tree; e.g. max(200,n)
-
-20 Sep 1991
-sample.f brought in sync with recent loess changes.
-
-24 Dec 1991
-l2fit.f fixed comment in single precision version
-
-10 Jan 1992
-ehg197.f new formula for approximating trL, valid for small f
-
-15 May 1992
-netlib/a/dloess now includes C drivers (written by William Shyu,
- adapted from code used inside the S system)
-
-22 Jun 1992
-ehg191.f Loop 11 ran too far, picking up one more value than necessary.
- The value was not used, so the loess computation itself is unaffected,
- but on some systems the old code could conceivably cause a reference
- to an invalid memory address and abort with a segmentation fault
- message.
-
-23 Jun 1992
-S.h #include <math.h>, since loessc.c calls floor() and pow().
-
-25 Mar 1996
-update email address
-
-4 Jul 2006
- bug fix ehg127() when k=1 and od>1
-
//GO.SYSIN DD changes
echo depend.ps 1>&2
-%!
-/Courier-Bold findfont 10 scalefont setfont
-%draw a box
-%x y width height box
-/box { newpath
- /height exch def
- /width exch def
- /y exch def
- /x exch def
- x width 2 div sub
- y height 2 div sub moveto
- width 0 rlineto
- 0 height rlineto
- width neg 0 rlineto
- closepath } def
-
-%draw a circle
-%x y radius circle
-/circle { newpath 0 360 arc } def
-
-%draw an ellipse
-%x y width height ellipse
-/ellipse { gsave
- /height exch def
- /width exch def
- 1 height width div scale
- width height div mul
- width 2 div
- circle stroke
- grestore } def
-
-%draw a centered label
-%x y str
-/label {
- /str exch def
- /y exch def
- /x exch def
- str stringwidth
- pop /width exch def
- x width 2 div sub
- y 10 3 div sub moveto str show
- } def
-
-%draw a line
-%x1 y1 x2 y2 drawline
-/drawline { 4 -2 roll moveto lineto stroke } def
-
-277 684 42 14 box stroke
-277 684 (lowesd) label
-349 630 42 14 box stroke
-349 630 (lowesf) label
-205 630 42 14 box stroke
-205 630 (lowesb) label
-155 565 42 14 box stroke
-155 565 (lowesr) label
-146 427 42 14 box stroke
-146 427 (lowese) label
-277 576 42 14 box stroke
-277 576 (lowesl) label
-203 464 42 14 box stroke
-203 464 (lofort) label
-81 576 42 14 box stroke
-81 576 (losave) label
-81 522 42 14 box stroke
-81 522 (lohead) label
-81 468 42 14 box stroke
-81 468 (loread) label
-405 540 42 14 box stroke
-405 540 (lowesa) label
-342 539 42 14 box stroke
-342 539 (lowesc) label
-92 461 134 434 drawline
-124.266363 435.502104 134.000000 434.000000 drawline
-134.000000 434.000000 128.592424 442.231532 drawline
-81 515 81 475 drawline
-77.000000 484.000000 81.000000 475.000000 drawline
-81.000000 475.000000 85.000000 484.000000 drawline
-81 569 81 529 drawline
-77.000000 538.000000 81.000000 529.000000 drawline
-81.000000 529.000000 85.000000 538.000000 drawline
-289 569 329 546 drawline
-319.203959 547.018615 329.000000 546.000000 drawline
-329.000000 546.000000 323.191728 553.953865 drawline
-154 558 146 434 drawline
-142.587739 443.238857 146.000000 434.000000 drawline
-146.000000 434.000000 150.571142 442.723799 drawline
-188 623 97 583 drawline
-103.629564 590.283466 97.000000 583.000000 drawline
-97.000000 583.000000 106.848776 582.959760 drawline
-204 623 203 471 drawline
-199.059296 480.026120 203.000000 471.000000 drawline
-203.000000 471.000000 207.059123 479.973490 drawline
-214 623 267 583 drawline
-257.406670 585.228906 267.000000 583.000000 drawline
-267.000000 583.000000 262.225925 591.614419 drawline
-199 623 160 572 drawline
-162.289620 581.579021 160.000000 572.000000 drawline
-160.000000 572.000000 168.644482 576.719420 drawline
-220 623 389 547 drawline
-379.151237 547.043173 389.000000 547.000000 drawline
-389.000000 547.000000 382.432359 554.339352 drawline
-202 623 148 434 drawline
-146.626394 443.752600 148.000000 434.000000 drawline
-148.000000 434.000000 154.318586 441.554831 drawline
-348 623 342 546 drawline
-338.711268 555.283547 342.000000 546.000000 drawline
-342.000000 546.000000 346.687091 554.662054 drawline
-353 623 400 547 drawline
-391.864262 552.550655 400.000000 547.000000 drawline
-400.000000 547.000000 398.668290 556.758409 drawline
-267 677 214 637 drawline
-218.774075 645.614419 214.000000 637.000000 drawline
-214.000000 637.000000 223.593330 639.228906 drawline
-286 677 339 637 drawline
-329.406670 639.228906 339.000000 637.000000 drawline
-339.000000 637.000000 334.225925 645.614419 drawline
-showpage
//GO.SYSIN DD depend.ps
echo sample.in 1>&2
-88 2 .6 2 .05
-12. 0.907 3.741 1
-12. 0.761 2.295 1
-12. 1.108 1.498 1
-12. 1.016 2.881 1
-12. 1.189 0.76 1
-9. 1.001 3.12 1
-9. 1.231 0.638 1
-9. 1.123 1.17 1
-12. 1.042 2.358 1
-12. 1.215 0.606 1
-12. 0.930 3.669 1
-12. 1.152 1.00 1
-15. 1.138 0.981 1
-18. 0.601 1.192 1
-7.5 0.696 0.926 1
-12. 0.686 1.59 1
-12. 1.072 1.806 1
-15. 1.074 1.962 1
-15. 0.934 4.028 1
-9. 0.808 3.148 1
-9. 1.071 1.836 1
-7.5 1.009 2.845 1
-7.5 1.142 1.013 1
-18. 1.229 0.414 1
-18. 1.175 0.812 1
-15. 0.568 0.374 1
-15. 0.977 3.623 1
-7.5 0.767 1.869 1
-7.5 1.006 2.836 1
-9. 0.893 3.567 1
-15. 1.152 0.866 1
-15. 0.693 1.369 1
-15. 1.232 0.542 1
-15. 1.036 2.739 1
-15. 1.125 1.20 1
-9. 1.081 1.719 1
-9. 0.868 3.423 1
-7.5 0.762 1.634 1
-7.5 1.144 1.021 1
-7.5 1.045 2.157 1
-18. 0.797 3.361 1
-18. 1.115 1.39 1
-18. 1.070 1.947 1
-18. 1.219 0.962 1
-9. 0.637 0.571 1
-9. 0.733 2.219 1
-9. 0.715 1.419 1
-9. 0.872 3.519 1
-7.5 0.765 1.732 1
-7.5 0.878 3.206 1
-7.5 0.811 2.471 1
-15. 0.676 1.777 1
-18. 1.045 2.571 1
-18. 0.968 3.952 1
-15. 0.846 3.931 1
-15. 0.684 1.587 1
-7.5 0.729 1.397 1
-7.5 0.911 3.536 1
-7.5 0.808 2.202 1
-7.5 1.168 0.756 1
-7.5 0.749 1.62 1
-7.5 0.892 3.656 1
-7.5 1.002 2.964 1
-18. 0.812 3.76 1
-18. 1.230 0.672 1
-18. 0.804 3.677 1
-12. 0.813 3.517 1
-12. 1.002 3.29 1
-9. 0.696 1.139 1
-9. 1.199 0.727 1
-9. 1.030 2.581 1
-15. 0.602 0.923 1
-15. 0.694 1.527 1
-15. 0.816 3.388 1
-15. 1.037 2.085 1
-15. 1.181 0.966 1
-7.5 0.899 3.488 1
-7.5 1.227 0.754 1
-9. 1.180 0.797 1
-7.5 0.795 2.064 1
-18. 0.990 3.732 1
-18. 1.201 0.586 1
-7.5 0.629 0.561 1
-9. 0.608 0.563 1
-12. 0.584 0.678 1
-15. 0.562 0.37 1
-18. 0.535 0.53 1
-18. 0.655 1.90 1
//GO.SYSIN DD sample.in
echo sample.f 1>&2
- integer d,i,tdeg,k,liv,lv,n,nvmax,tau0,nf,livm,ivm,nm
- parameter (nm=1000)
- parameter (livm=10000)
- parameter (ivm=20000)
- integer iwork(livm)
- real alpha,fcell
- real s(nm),w(nm),work(ivm),xx(nm,8),yy(nm)
- integer ifloor
- external ifloor
- read(5,*)n,d,alpha,tdeg,fcell
- nf=min(n,ifloor(n*alpha))
- nvmax=max(200,n)
- tau0=d+1
- if(tdeg.eq.2) tau0=((d+2)*(d+1))/2
- liv=50+(2**d+4)*nvmax+2*n
- lv =50+(3*d+3)*nvmax+n+(tau0+2)*nf
- if(nm.lt.n .or. livm.lt.liv .or. ivm.lt.lv)then
- print *, 'sample program needs more workspace'
- else
- call getxy(n,d,xx,yy,w)
- call lowesd(106,iwork,liv,lv,work,d,n,alpha,
- $ tdeg,nvmax,.false.)
- work(2)=fcell
- call lowesb(xx,yy,w,xx,.false.,iwork,liv,lv,work)
- call lowese(iwork,liv,lv,work,n,xx,s)
- do 4 i=1,n
- print *, s(i)
- 4 continue
- end if
- end
//GO.SYSIN DD sample.f
echo getxy.f 1>&2
- subroutine getxy(n,d,x,y,w)
- integer n, d, i, j
- real x(n,d), y(n), w(n)
- do 100 i=1,n
- read(5,*,end=110) (x(i,j),j=1,d),y(i),w(i)
-100 continue
- return
-110 print *, 'unexpected eof i=',i,' n=',n
- stop
- end
//GO.SYSIN DD getxy.f
echo bbox.f 1>&2
- subroutine ehg126(d,n,vc,x,v,nvmax)
- integer d,execnt,i,j,k,n,nv,nvmax,vc
- real machin,alpha,beta,mu,t
- real v(nvmax,d),x(n,d)
- real r1mach
- external r1mach
- save machin,execnt
- data execnt /0/
-c MachInf -> machin
- execnt=execnt+1
- if(execnt.eq.1)then
- machin=r1mach(2)
- end if
-c fill in vertices for bounding box of $x$
-c lower left, upper right
- do 3 k=1,d
- alpha=machin
- beta=-machin
- do 4 i=1,n
- t=x(i,k)
- alpha=min(alpha,t)
- beta=max(beta,t)
- 4 continue
-c expand the box a little
- mu=0.005*max(beta-alpha,1.e-10*max(abs(alpha),abs(beta))+1.e-30
- $)
- alpha=alpha-mu
- beta=beta+mu
- v(1,k)=alpha
- v(vc,k)=beta
- 3 continue
-c remaining vertices
- do 5 i=2,vc-1
- j=i-1
- do 6 k=1,d
- v(i,k)=v(1+mod(j,2)*(vc-1),k)
- j=float(j)/2.
- 6 continue
- 5 continue
- return
- end
//GO.SYSIN DD bbox.f
echo cpvert.f 1>&2
- subroutine ehg125(p,nv,v,vhit,nvmax,d,k,t,r,s,f,l,u)
- logical i1,i2,match
- integer d,execnt,h,i,i3,j,k,m,mm,nv,nvmax,p,r,s
- integer f(r,0:1,s),l(r,0:1,s),u(r,0:1,s),vhit(nvmax)
- real t
- real v(nvmax,d)
- external ehg182
- save execnt
- data execnt /0/
- execnt=execnt+1
- h=nv
- do 3 i=1,r
- do 4 j=1,s
- h=h+1
- do 5 i3=1,d
- v(h,i3)=v(f(i,0,j),i3)
- 5 continue
- v(h,k)=t
-c check for redundant vertex
- match=.false.
- m=1
-c top of while loop
- 6 if(.not.match)then
- i1=(m.le.nv)
- else
- i1=.false.
- end if
- if(.not.(i1))goto 7
- match=(v(m,1).eq.v(h,1))
- mm=2
-c top of while loop
- 8 if(match)then
- i2=(mm.le.d)
- else
- i2=.false.
- end if
- if(.not.(i2))goto 9
- match=(v(m,mm).eq.v(h,mm))
- mm=mm+1
- goto 8
-c bottom of while loop
- 9 m=m+1
- goto 6
-c bottom of while loop
- 7 m=m-1
- if(match)then
- h=h-1
- else
- m=h
- if(vhit(1).ge.0)then
- vhit(m)=p
- end if
- end if
- l(i,0,j)=f(i,0,j)
- l(i,1,j)=m
- u(i,0,j)=m
- u(i,1,j)=f(i,1,j)
- 4 continue
- 3 continue
- nv=h
- if(.not.(nv.le.nvmax))then
- call ehg182(180)
- end if
- return
- end
//GO.SYSIN DD cpvert.f
echo descend.f 1>&2
- integer function ehg138(i,z,a,xi,lo,hi,ncmax)
- logical i1
- integer d,execnt,i,j,nc,ncmax
- integer a(ncmax),hi(ncmax),lo(ncmax)
- real xi(ncmax),z(8)
- save execnt
- data execnt /0/
- execnt=execnt+1
-c descend tree until leaf or ambiguous
- j=i
-c top of while loop
- 3 if(a(j).ne.0)then
- i1=(z(a(j)).ne.xi(j))
- else
- i1=.false.
- end if
- if(.not.(i1))goto 4
- if(z(a(j)).le.xi(j))then
- j=lo(j)
- else
- j=hi(j)
- end if
- goto 3
-c bottom of while loop
- 4 ehg138=j
- return
- end
//GO.SYSIN DD descend.f
echo ehg131.f 1>&2
- subroutine ehg131(x,y,rw,trl,diagl,kernel,k,n,d,nc,ncmax,vc,nv,nvm
- $ax,nf,f,a,c,hi,lo,pi,psi,v,vhit,vval,xi,dist,eta,b,ntol,fd,w,vval2
- $,rcond,sing,dd,tdeg,cdeg,lq,lf,setlf)
- logical setlf
- integer identi,d,dd,execnt,i1,i2,j,k,kernel,n,nc,ncmax,nf,ntol,nv,
- $nvmax,sing,tdeg,vc
- integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),cdeg(8),hi(ncmax),lo(ncm
- $ax),pi(n),psi(n),vhit(nvmax)
- real f,fd,rcond,trl
- real lf(0:d,nvmax,nf),b(*),delta(8),diagl(n),dist(n),eta(nf),rw(n)
- $,v(nvmax,d),vval(0:d,nvmax),vval2(0:d,nvmax),w(nf),x(n,d),xi(ncmax
- $),y(n)
- external ehg126,ehg182,ehg139,ehg124
- real snrm2
- external snrm2
- save execnt
- data execnt /0/
-c Identity -> identi
-c X -> b
- execnt=execnt+1
- if(.not.(d.le.8))then
- call ehg182(101)
- end if
-c build $k$-d tree
- call ehg126(d,n,vc,x,v,nvmax)
- nv=vc
- nc=1
- do 3 j=1,vc
- c(j,nc)=j
- vhit(j)=0
- 3 continue
- do 4 i1=1,d
- delta(i1)=v(vc,i1)-v(1,i1)
- 4 continue
- fd=fd*snrm2(d,delta,1)
- do 5 identi=1,n
- pi(identi)=identi
- 5 continue
- call ehg124(1,n,d,n,nv,nc,ncmax,vc,x,pi,a,xi,lo,hi,c,v,vhit,nvmax,
- $ntol,fd,dd)
-c smooth
- if(trl.ne.0)then
- do 6 i2=1,nv
- do 7 i1=0,d
- vval2(i1,i2)=0
- 7 continue
- 6 continue
- end if
- call ehg139(v,nvmax,nv,n,d,nf,f,x,pi,psi,y,rw,trl,kernel,k,dist,di
- $st,eta,b,d,w,diagl,vval2,nc,vc,a,xi,lo,hi,c,vhit,rcond,sing,dd,tde
- $g,cdeg,lq,lf,setlf,vval)
- return
- end
//GO.SYSIN DD ehg131.f
echo ehg133.f 1>&2
- subroutine ehg133(n,d,vc,nvmax,nc,ncmax,a,c,hi,lo,v,vval,xi,m,z,s)
- integer d,execnt,i,i1,m,nc,ncmax,nv,nvmax,vc
- integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax)
- real delta(8),s(m),v(nvmax,d),vval(0:d,nvmax),xi(ncmax),z(m,d)
- real ehg128
- external ehg128
- save execnt
- data execnt /0/
- execnt=execnt+1
- do 3 i=1,m
- do 4 i1=1,d
- delta(i1)=z(i,i1)
- 4 continue
- s(i)=ehg128(delta,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vval)
- 3 continue
- return
- end
//GO.SYSIN DD ehg133.f
echo eval.f 1>&2
- real function ehg128(z,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vval)
- logical i10,i2,i3,i4,i5,i6,i7,i8,i9
- integer d,execnt,i,i1,i11,i12,ig,ii,j,lg,ll,m,nc,ncmax,nt,nv,nvmax
- $,ur,vc
- integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax),t(20)
- real ge,gn,gs,gw,gpe,gpn,gps,gpw,h,phi0,phi1,psi0,psi1,s,sew,sns,v
- $0,v1,xibar
- real g(0:8,256),g0(0:8),g1(0:8),v(nvmax,d),vval(0:d,nvmax),xi(ncma
- $x),z(d)
- external ehg182,ehg184
- save execnt
- data execnt /0/
- execnt=execnt+1
-c locate enclosing cell
- nt=1
- t(nt)=1
- j=1
-c top of while loop
- 3 if(.not.(a(j).ne.0))goto 4
- nt=nt+1
-c bug fix 2006-07-18 (thanks, btyner@gmail.com)
- if(z(a(j)).le.xi(j))then
- i1=lo(j)
- else
- i1=hi(j)
- end if
- t(nt)=i1
- if(.not.(nt.lt.20))then
- call ehg182(181)
- end if
- j=t(nt)
- goto 3
-c bottom of while loop
- 4 continue
-c tensor
- do 5 i12=1,vc
- do 6 i11=0,d
- g(i11,i12)=vval(i11,c(i12,j))
- 6 continue
- 5 continue
- lg=vc
- ll=c(1,j)
- ur=c(vc,j)
- do 7 i=d,1,-1
- h=(z(i)-v(ll,i))/(v(ur,i)-v(ll,i))
- if(h.lt.-.001)then
- call ehg184('eval ',z,d,1)
- call ehg184('lowerlimit ',v(ll,1),d,nvmax)
- else
- if(1.001.lt.h)then
- call ehg184('eval ',z,d,1)
- call ehg184('upperlimit ',v(ur,1),d,nvmax)
- end if
- end if
- if(-.001.le.h)then
- i2=(h.le.1.001)
- else
- i2=.false.
- end if
- if(.not.i2)then
- call ehg182(122)
- end if
- lg=float(lg)/2.
- do 8 ig=1,lg
-c Hermite basis
- phi0=(1-h)**2*(1+2*h)
- phi1=h**2*(3-2*h)
- psi0=h*(1-h)**2
- psi1=h**2*(h-1)
- g(0,ig)=phi0*g(0,ig)+phi1*g(0,ig+lg)+(psi0*g(i,ig)+psi1*g(i,
- $ig+lg))*(v(ur,i)-v(ll,i))
- do 9 ii=1,i-1
- g(ii,ig)=phi0*g(ii,ig)+phi1*g(ii,ig+lg)
- 9 continue
- 8 continue
- 7 continue
- s=g(0,1)
-c blending
- if(d.eq.2)then
-c ----- North -----
- v0=v(ll,1)
- v1=v(ur,1)
- do 10 i11=0,d
- g0(i11)=vval(i11,c(3,j))
- 10 continue
- do 11 i11=0,d
- g1(i11)=vval(i11,c(4,j))
- 11 continue
- xibar=v(ur,2)
- m=nt-1
-c top of while loop
- 12 if(m.eq.0)then
- i4=.true.
- else
- if(a(t(m)).eq.2)then
- i3=(xi(t(m)).eq.xibar)
- else
- i3=.false.
- end if
- i4=i3
- end if
- if(.not.(.not.i4))goto 13
- m=m-1
-c voidp junk
- goto 12
-c bottom of while loop
- 13 if(m.ge.1)then
- m=hi(t(m))
-c top of while loop
- 14 if(.not.(a(m).ne.0))goto 15
- if(z(a(m)).le.xi(m))then
- m=lo(m)
- else
- m=hi(m)
- end if
- goto 14
-c bottom of while loop
- 15 if(v0.lt.v(c(1,m),1))then
- v0=v(c(1,m),1)
- do 16 i11=0,d
- g0(i11)=vval(i11,c(1,m))
- 16 continue
- end if
- if(v(c(2,m),1).lt.v1)then
- v1=v(c(2,m),1)
- do 17 i11=0,d
- g1(i11)=vval(i11,c(2,m))
- 17 continue
- end if
- end if
- h=(z(1)-v0)/(v1-v0)
-c Hermite basis
- phi0=(1-h)**2*(1+2*h)
- phi1=h**2*(3-2*h)
- psi0=h*(1-h)**2
- psi1=h**2*(h-1)
- gn=phi0*g0(0)+phi1*g1(0)+(psi0*g0(1)+psi1*g1(1))*(v1-v0)
- gpn=phi0*g0(2)+phi1*g1(2)
-c ----- South -----
- v0=v(ll,1)
- v1=v(ur,1)
- do 18 i11=0,d
- g0(i11)=vval(i11,c(1,j))
- 18 continue
- do 19 i11=0,d
- g1(i11)=vval(i11,c(2,j))
- 19 continue
- xibar=v(ll,2)
- m=nt-1
-c top of while loop
- 20 if(m.eq.0)then
- i6=.true.
- else
- if(a(t(m)).eq.2)then
- i5=(xi(t(m)).eq.xibar)
- else
- i5=.false.
- end if
- i6=i5
- end if
- if(.not.(.not.i6))goto 21
- m=m-1
-c voidp junk
- goto 20
-c bottom of while loop
- 21 if(m.ge.1)then
- m=lo(t(m))
-c top of while loop
- 22 if(.not.(a(m).ne.0))goto 23
- if(z(a(m)).le.xi(m))then
- m=lo(m)
- else
- m=hi(m)
- end if
- goto 22
-c bottom of while loop
- 23 if(v0.lt.v(c(3,m),1))then
- v0=v(c(3,m),1)
- do 24 i11=0,d
- g0(i11)=vval(i11,c(3,m))
- 24 continue
- end if
- if(v(c(4,m),1).lt.v1)then
- v1=v(c(4,m),1)
- do 25 i11=0,d
- g1(i11)=vval(i11,c(4,m))
- 25 continue
- end if
- end if
- h=(z(1)-v0)/(v1-v0)
-c Hermite basis
- phi0=(1-h)**2*(1+2*h)
- phi1=h**2*(3-2*h)
- psi0=h*(1-h)**2
- psi1=h**2*(h-1)
- gs=phi0*g0(0)+phi1*g1(0)+(psi0*g0(1)+psi1*g1(1))*(v1-v0)
- gps=phi0*g0(2)+phi1*g1(2)
-c ----- East -----
- v0=v(ll,2)
- v1=v(ur,2)
- do 26 i11=0,d
- g0(i11)=vval(i11,c(2,j))
- 26 continue
- do 27 i11=0,d
- g1(i11)=vval(i11,c(4,j))
- 27 continue
- xibar=v(ur,1)
- m=nt-1
-c top of while loop
- 28 if(m.eq.0)then
- i8=.true.
- else
- if(a(t(m)).eq.1)then
- i7=(xi(t(m)).eq.xibar)
- else
- i7=.false.
- end if
- i8=i7
- end if
- if(.not.(.not.i8))goto 29
- m=m-1
-c voidp junk
- goto 28
-c bottom of while loop
- 29 if(m.ge.1)then
- m=hi(t(m))
-c top of while loop
- 30 if(.not.(a(m).ne.0))goto 31
- if(z(a(m)).le.xi(m))then
- m=lo(m)
- else
- m=hi(m)
- end if
- goto 30
-c bottom of while loop
- 31 if(v0.lt.v(c(1,m),2))then
- v0=v(c(1,m),2)
- do 32 i11=0,d
- g0(i11)=vval(i11,c(1,m))
- 32 continue
- end if
- if(v(c(3,m),2).lt.v1)then
- v1=v(c(3,m),2)
- do 33 i11=0,d
- g1(i11)=vval(i11,c(3,m))
- 33 continue
- end if
- end if
- h=(z(2)-v0)/(v1-v0)
-c Hermite basis
- phi0=(1-h)**2*(1+2*h)
- phi1=h**2*(3-2*h)
- psi0=h*(1-h)**2
- psi1=h**2*(h-1)
- ge=phi0*g0(0)+phi1*g1(0)+(psi0*g0(2)+psi1*g1(2))*(v1-v0)
- gpe=phi0*g0(1)+phi1*g1(1)
-c ----- West -----
- v0=v(ll,2)
- v1=v(ur,2)
- do 34 i11=0,d
- g0(i11)=vval(i11,c(1,j))
- 34 continue
- do 35 i11=0,d
- g1(i11)=vval(i11,c(3,j))
- 35 continue
- xibar=v(ll,1)
- m=nt-1
-c top of while loop
- 36 if(m.eq.0)then
- i10=.true.
- else
- if(a(t(m)).eq.1)then
- i9=(xi(t(m)).eq.xibar)
- else
- i9=.false.
- end if
- i10=i9
- end if
- if(.not.(.not.i10))goto 37
- m=m-1
-c voidp junk
- goto 36
-c bottom of while loop
- 37 if(m.ge.1)then
- m=lo(t(m))
-c top of while loop
- 38 if(.not.(a(m).ne.0))goto 39
- if(z(a(m)).le.xi(m))then
- m=lo(m)
- else
- m=hi(m)
- end if
- goto 38
-c bottom of while loop
- 39 if(v0.lt.v(c(2,m),2))then
- v0=v(c(2,m),2)
- do 40 i11=0,d
- g0(i11)=vval(i11,c(2,m))
- 40 continue
- end if
- if(v(c(4,m),2).lt.v1)then
- v1=v(c(4,m),2)
- do 41 i11=0,d
- g1(i11)=vval(i11,c(4,m))
- 41 continue
- end if
- end if
- h=(z(2)-v0)/(v1-v0)
-c Hermite basis
- phi0=(1-h)**2*(1+2*h)
- phi1=h**2*(3-2*h)
- psi0=h*(1-h)**2
- psi1=h**2*(h-1)
- gw=phi0*g0(0)+phi1*g1(0)+(psi0*g0(2)+psi1*g1(2))*(v1-v0)
- gpw=phi0*g0(1)+phi1*g1(1)
-c NS
- h=(z(2)-v(ll,2))/(v(ur,2)-v(ll,2))
-c Hermite basis
- phi0=(1-h)**2*(1+2*h)
- phi1=h**2*(3-2*h)
- psi0=h*(1-h)**2
- psi1=h**2*(h-1)
- sns=phi0*gs+phi1*gn+(psi0*gps+psi1*gpn)*(v(ur,2)-v(ll,2))
-c EW
- h=(z(1)-v(ll,1))/(v(ur,1)-v(ll,1))
-c Hermite basis
- phi0=(1-h)**2*(1+2*h)
- phi1=h**2*(3-2*h)
- psi0=h*(1-h)**2
- psi1=h**2*(h-1)
- sew=phi0*gw+phi1*ge+(psi0*gpw+psi1*gpe)*(v(ur,1)-v(ll,1))
- s=(sns+sew)-s
- end if
- ehg128=s
- return
- end
//GO.SYSIN DD eval.f
echo f77.f 1>&2
- integer function ifloor(x)
- real x
- ifloor=x
- if(ifloor.gt.x) ifloor=ifloor-1
- end
- real function sign(a1,a2)
- real a1, a2
- sign=abs(a1)
- if(a2.ge.0) sign=-sign
- end
//GO.SYSIN DD f77.f
echo l2fit.f 1>&2
- subroutine ehg136(u,lm,m,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,o
- $d,o,ihat,w,rcond,sing,dd,tdeg,cdeg,s)
- integer identi,d,dd,execnt,i,i1,ihat,info,j,k,kernel,l,lm,m,n,nf,o
- $d,sing,tdeg
- integer cdeg(8),psi(n)
- real f,i2,rcond,scale,tol
- real o(m,n),sigma(15),e(15,15),g(15,15),b(nf,k),dist(n),eta(nf),ga
- $mma(15),q(8),qraux(15),rw(n),s(0:od,m),u(lm,d),w(nf),work(15),x(n,
- $d),y(n)
- external ehg127,ehg182,sqrsl
- real sdot
- external sdot
- save execnt
- data execnt /0/
-c V -> g
-c U -> e
-c Identity -> identi
-c L -> o
-c X -> b
- execnt=execnt+1
- if(.not.(k.le.nf-1))then
- call ehg182(104)
- end if
- if(.not.(k.le.15))then
- call ehg182(105)
- end if
- do 3 identi=1,n
- psi(identi)=identi
- 3 continue
- do 4 l=1,m
- do 5 i1=1,d
- q(i1)=u(l,i1)
- 5 continue
- call ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w,rcon
- $d,sing,sigma,e,g,gamma,qraux,work,tol,dd,tdeg,cdeg,s(0,l))
- if(ihat.eq.1)then
-c $L sub {l,l} =
-c V sub {1,:} SIGMA sup {+} U sup T
-c (Q sup T W e sub i )$
- if(.not.(m.eq.n))then
- call ehg182(123)
- end if
-c find $i$ such that $l = psi sub i$
- i=1
-c top of while loop
- 6 if(.not.(l.ne.psi(i)))goto 7
- i=i+1
- if(.not.(i.lt.nf))then
- call ehg182(123)
- end if
- goto 6
-c bottom of while loop
- 7 do 8 i1=1,nf
- eta(i1)=0
- 8 continue
- eta(i)=w(i)
-c $eta = Q sup T W e sub i$
- call sqrsl(b,nf,nf,k,qraux,eta,eta,eta,eta,eta,eta,1000,info
- $)
-c $gamma = U sup T eta sub {1:k}$
- do 9 i1=1,k
- gamma(i1)=0
- 9 continue
- do 10 j=1,k
- i2=eta(j)
- do 11 i1=1,k
- gamma(i1)=gamma(i1)+i2*e(j,i1)
- 11 continue
- 10 continue
-c $gamma = SIGMA sup {+} gamma$
- do 12 j=1,k
- if(tol.lt.sigma(j))then
- gamma(j)=gamma(j)/sigma(j)
- else
- gamma(j)=0.
- end if
- 12 continue
-c voidp junk
-c voidp junk
- o(l,1)=sdot(k,g(1,1),15,gamma,1)
- else
- if(ihat.eq.2)then
-c $L sub {l,:} =
-c V sub {1,:} SIGMA sup {+}
-c ( U sup T Q sup T ) W $
- do 13 i1=1,n
- o(l,i1)=0
- 13 continue
- do 14 j=1,k
- do 15 i1=1,nf
- eta(i1)=0
- 15 continue
- do 16 i1=1,k
- eta(i1)=e(i1,j)
- 16 continue
- call sqrsl(b,nf,nf,k,qraux,eta,eta,work,work,work,work
- $,10000,info)
- if(tol.lt.sigma(j))then
- scale=1./sigma(j)
- else
- scale=0.
- end if
- do 17 i1=1,nf
- eta(i1)=eta(i1)*(scale*w(i1))
- 17 continue
- do 18 i=1,nf
- o(l,psi(i))=o(l,psi(i))+g(1,j)*eta(i)
- 18 continue
- 14 continue
- end if
- end if
- 4 continue
- return
- end
//GO.SYSIN DD l2fit.f
echo l2tr.f 1>&2
- subroutine ehg139(v,nvmax,nv,n,d,nf,f,x,pi,psi,y,rw,trl,kernel,k,d
- $ist,phi,eta,b,od,w,diagl,vval2,ncmax,vc,a,xi,lo,hi,c,vhit,rcond,si
- $ng,dd,tdeg,cdeg,lq,lf,setlf,s)
- logical setlf
- integer identi,d,dd,execnt,i,i2,i3,i5,i6,ii,ileaf,info,j,k,kernel,
- $l,n,nc,ncmax,nf,nleaf,nv,nvmax,od,sing,tdeg,vc
- integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),cdeg(8),hi(ncmax),leaf(2
- $56),lo(ncmax),phi(n),pi(n),psi(n),vhit(nvmax)
- real f,i1,i4,i7,rcond,scale,term,tol,trl
- real lf(0:d,nvmax,nf),sigma(15),u(15,15),e(15,15),b(nf,k),diagl(n)
- $,dist(n),eta(nf),gamma(15),q(8),qraux(15),rw(n),s(0:od,nv),v(nvmax
- $,d),vval2(0:d,nv),w(nf),work(15),x(n,d),xi(ncmax),y(n),z(8)
- external ehg127,ehg182,sqrsl,ehg137
- real ehg128
- external ehg128
- real sdot
- external sdot
- save execnt
- data execnt /0/
-c V -> e
-c Identity -> identi
-c X -> b
- execnt=execnt+1
-c l2fit with trace(L)
- if(.not.(k.le.nf-1))then
- call ehg182(104)
- end if
- if(.not.(k.le.15))then
- call ehg182(105)
- end if
- if(trl.ne.0)then
- do 3 i5=1,n
- diagl(i5)=0
- 3 continue
- do 4 i6=1,nv
- do 5 i5=0,d
- vval2(i5,i6)=0
- 5 continue
- 4 continue
- end if
- do 6 identi=1,n
- psi(identi)=identi
- 6 continue
- do 7 l=1,nv
- do 8 i5=1,d
- q(i5)=v(l,i5)
- 8 continue
- call ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w,rcon
- $d,sing,sigma,u,e,gamma,qraux,work,tol,dd,tdeg,cdeg,s(0,l))
- if(trl.ne.0)then
-c invert $psi$
- do 9 i5=1,n
- phi(i5)=0
- 9 continue
- do 10 i=1,nf
- phi(psi(i))=i
- 10 continue
- do 11 i5=1,d
- z(i5)=v(l,i5)
- 11 continue
- call ehg137(z,vhit(l),leaf,nleaf,d,nv,nvmax,ncmax,vc,a,xi,lo
- $,hi,c,v)
- do 12 ileaf=1,nleaf
- do 13 ii=lo(leaf(ileaf)),hi(leaf(ileaf))
- i=phi(pi(ii))
- if(i.ne.0)then
- if(.not.(psi(i).eq.pi(ii)))then
- call ehg182(194)
- end if
- do 14 i5=1,nf
- eta(i5)=0
- 14 continue
- eta(i)=w(i)
-c $eta = Q sup T W e sub i$
- call sqrsl(b,nf,nf,k,qraux,eta,work,eta,eta,work,wo
- $rk,1000,info)
- do 15 j=1,k
- if(tol.lt.sigma(j))then
- i4=sdot(k,u(1,j),1,eta,1)/sigma(j)
- else
- i4=0.
- end if
- gamma(j)=i4
- 15 continue
- do 16 j=1,d+1
-c bug fix 2006-07-15 for k=1, od>1. (thanks btyner@gmail.com)
- if(j.le.k)then
- vval2(j-1,l)=sdot(k,e(j,1),15,gamma,1)
- else
- vval2(j-1,l)=0
- end if
- 16 continue
- do 17 i5=1,d
- z(i5)=x(pi(ii),i5)
- 17 continue
- term=ehg128(z,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vval2
- $)
- diagl(pi(ii))=diagl(pi(ii))+term
- do 18 i5=0,d
- vval2(i5,l)=0
- 18 continue
- end if
- 13 continue
- 12 continue
- end if
- if(setlf)then
-c $Lf sub {:,l,:} = V SIGMA sup {+} U sup T Q sup T W$
- if(.not.(k.ge.d+1))then
- call ehg182(196)
- end if
- do 19 i5=1,nf
- lq(l,i5)=psi(i5)
- 19 continue
- do 20 i6=1,nf
- do 21 i5=0,d
- lf(i5,l,i6)=0
- 21 continue
- 20 continue
- do 22 j=1,k
- do 23 i5=1,nf
- eta(i5)=0
- 23 continue
- do 24 i5=1,k
- eta(i5)=u(i5,j)
- 24 continue
- call sqrsl(b,nf,nf,k,qraux,eta,eta,work,work,work,work,10
- $000,info)
- if(tol.lt.sigma(j))then
- scale=1./sigma(j)
- else
- scale=0.
- end if
- do 25 i5=1,nf
- eta(i5)=eta(i5)*(scale*w(i5))
- 25 continue
- do 26 i=1,nf
- i7=eta(i)
- do 27 i5=0,d
- if(i5.lt.k)then
- lf(i5,l,i)=lf(i5,l,i)+e(1+i5,j)*i7
- else
- lf(i5,l,i)=0
- end if
- 27 continue
- 26 continue
- 22 continue
- end if
- 7 continue
- if(trl.ne.0)then
- if(n.le.0)then
- trl=0.
- else
- i3=n
- i1=diagl(i3)
- do 28 i2=i3-1,1,-1
- i1=diagl(i2)+i1
- 28 continue
- trl=i1
- end if
- end if
- return
- end
//GO.SYSIN DD l2tr.f
echo lowesb.f 1>&2
- subroutine lowesb(xx,yy,ww,diagl,infl,iv,liv,lv,wv)
- logical infl,setlf
- integer execnt
- integer iv(*)
- real trl
- real diagl(*),wv(*),ww(*),xx(*),yy(*)
- external ehg131,ehg182,ehg183
- integer ifloor
- external ifloor
- save execnt
- data execnt /0/
- execnt=execnt+1
- if(.not.(iv(28).ne.173))then
- call ehg182(174)
- end if
- if(iv(28).ne.172)then
- if(.not.(iv(28).eq.171))then
- call ehg182(171)
- end if
- end if
- iv(28)=173
- if(infl)then
- trl=1.
- else
- trl=0.
- end if
- setlf=(iv(27).ne.iv(25))
- call ehg131(xx,yy,ww,trl,diagl,iv(20),iv(29),iv(3),iv(2),iv(5),iv(
- $17),iv(4),iv(6),iv(14),iv(19),wv(1),iv(iv(7)),iv(iv(8)),iv(iv(9)),
- $iv(iv(10)),iv(iv(22)),iv(iv(27)),wv(iv(11)),iv(iv(23)),wv(iv(13)),
- $wv(iv(12)),wv(iv(15)),wv(iv(16)),wv(iv(18)),ifloor(iv(3)*wv(2)),wv
- $(3),wv(iv(26)),wv(iv(24)),wv(4),iv(30),iv(33),iv(32),iv(41),iv(iv(
- $25)),wv(iv(34)),setlf)
- if(iv(14).lt.iv(6)+float(iv(4))/2.)then
- call ehg183('Warning. k-d tree limited by memory; nvmax=',iv(14
- $),1,1)
- else
- if(iv(17).lt.iv(5)+2)then
- call ehg183('Warning. k-d tree limited by memory. ncmax=',iv
- $(17),1,1)
- end if
- end if
- return
- end
//GO.SYSIN DD lowesb.f
echo lowesd.f 1>&2
- subroutine lowesd(versio,iv,liv,lv,v,d,n,f,ideg,nvmax,setlf)
- logical setlf
- integer bound,d,execnt,i,i1,i2,ideg,j,liv,lv,n,ncmax,nf,nvmax,vc,v
- $ersio
- integer iv(liv)
- real f
- real v(lv)
- external ehg182
- integer ifloor
- external ifloor
- save execnt
- data execnt /0/
-c version -> versio
- execnt=execnt+1
- if(.not.(versio.eq.106))then
- call ehg182(100)
- end if
- iv(28)=171
- iv(2)=d
- iv(3)=n
- vc=2**d
- iv(4)=vc
- if(.not.(0.lt.f))then
- call ehg182(120)
- end if
- nf=min(n,ifloor(n*f))
- iv(19)=nf
- iv(20)=1
- if(ideg.eq.0)then
- i1=1
- else
- if(ideg.eq.1)then
- i1=d+1
- else
- if(ideg.eq.2)then
- i1=float((d+2)*(d+1))/2.
- end if
- end if
- end if
- iv(29)=i1
- iv(21)=1
- iv(14)=nvmax
- ncmax=nvmax
- iv(17)=ncmax
- iv(30)=0
- iv(32)=ideg
- if(.not.(ideg.ge.0))then
- call ehg182(195)
- end if
- if(.not.(ideg.le.2))then
- call ehg182(195)
- end if
- iv(33)=d
- do 3 i2=41,49
- iv(i2)=ideg
- 3 continue
- iv(7)=50
- iv(8)=iv(7)+ncmax
- iv(9)=iv(8)+vc*ncmax
- iv(10)=iv(9)+ncmax
- iv(22)=iv(10)+ncmax
-c initialize permutation
- j=iv(22)-1
- do 4 i=1,n
- iv(j+i)=i
- 4 continue
- iv(23)=iv(22)+n
- iv(25)=iv(23)+nvmax
- if(setlf)then
- iv(27)=iv(25)+nvmax*nf
- else
- iv(27)=iv(25)
- end if
- bound=iv(27)+n
- if(.not.(bound-1.le.liv))then
- call ehg182(102)
- end if
- iv(11)=50
- iv(13)=iv(11)+nvmax*d
- iv(12)=iv(13)+(d+1)*nvmax
- iv(15)=iv(12)+ncmax
- iv(16)=iv(15)+n
- iv(18)=iv(16)+nf
- iv(24)=iv(18)+iv(29)*nf
- iv(34)=iv(24)+(d+1)*nvmax
- if(setlf)then
- iv(26)=iv(34)+(d+1)*nvmax*nf
- else
- iv(26)=iv(34)
- end if
- bound=iv(26)+nf
- if(.not.(bound-1.le.lv))then
- call ehg182(103)
- end if
- v(1)=f
- v(2)=0.05
- v(3)=0.
- v(4)=1.
- return
- end
//GO.SYSIN DD lowesd.f
echo lowese.f 1>&2
- subroutine lowese(iv,liv,lv,wv,m,z,s)
- integer execnt,m
- integer iv(*)
- real s(m),wv(*),z(m,1)
- external ehg133,ehg182
- save execnt
- data execnt /0/
- execnt=execnt+1
- if(.not.(iv(28).ne.172))then
- call ehg182(172)
- end if
- if(.not.(iv(28).eq.173))then
- call ehg182(173)
- end if
- call ehg133(iv(3),iv(2),iv(4),iv(14),iv(5),iv(17),iv(iv(7)),iv(iv(
- $8)),iv(iv(9)),iv(iv(10)),wv(iv(11)),wv(iv(13)),wv(iv(12)),m,z,s)
- return
- end
//GO.SYSIN DD lowese.f
echo lowesf.f 1>&2
- subroutine lowesf(xx,yy,ww,iv,liv,lv,wv,m,z,l,ihat,s)
- logical i1
- integer execnt,ihat,m,n
- integer iv(*)
- real l(m,*),s(m),wv(*),ww(*),xx(*),yy(*),z(m,1)
- external ehg182,ehg136
- save execnt
- data execnt /0/
- execnt=execnt+1
- if(171.le.iv(28))then
- i1=(iv(28).le.174)
- else
- i1=.false.
- end if
- if(.not.i1)then
- call ehg182(171)
- end if
- iv(28)=172
- if(.not.(iv(14).ge.iv(19)))then
- call ehg182(186)
- end if
- call ehg136(z,m,m,iv(3),iv(2),iv(19),wv(1),xx,iv(iv(22)),yy,ww,iv(
- $20),iv(29),wv(iv(15)),wv(iv(16)),wv(iv(18)),0,l,ihat,wv(iv(26)),wv
- $(4),iv(30),iv(33),iv(32),iv(41),s)
- return
- end
//GO.SYSIN DD lowesf.f
echo rbuild.f 1>&2
- subroutine ehg124(ll,uu,d,n,nv,nc,ncmax,vc,x,pi,a,xi,lo,hi,c,v,vhi
- $t,nvmax,fc,fd,dd)
- logical i1,i2,i3,leaf
- integer d,dd,execnt,fc,i4,inorm2,k,l,ll,m,n,nc,ncmax,nv,nvmax,p,u,
- $uu,vc,lower,upper,check,offset
- integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax),pi(n),vhit(nvmax)
- real diam,fd
- real diag(8),sigma(8),v(nvmax,d),x(n,d),xi(ncmax)
- external ehg125,ehg106,ehg129
- integer isamax
- external isamax
- save execnt
- data execnt /0/
- execnt=execnt+1
- p=1
- l=ll
- u=uu
- lo(p)=l
- hi(p)=u
-c top of while loop
- 3 if(.not.(p.le.nc))goto 4
- do 5 i4=1,dd
- diag(i4)=v(c(vc,p),i4)-v(c(1,p),i4)
- 5 continue
- diam=0
- do 6 inorm2=1,dd
- diam=diam+diag(inorm2)**2
- 6 continue
- diam=sqrt(diam)
- if((u-l)+1.le.fc)then
- i1=.true.
- else
- i1=(diam.le.fd)
- end if
- if(i1)then
- leaf=.true.
- else
- if(ncmax.lt.nc+2)then
- i2=.true.
- else
- i2=(nvmax.lt.nv+float(vc)/2.)
- end if
- leaf=i2
- end if
- if(.not.leaf)then
- call ehg129(l,u,dd,x,pi,n,sigma)
- k=isamax(dd,sigma,1)
- m=float(l+u)/2.
- call ehg106(l,u,m,1,x(1,k),pi,n)
-
-c bug fix from btyner@gmail.com 2006-07-20
- offset = 0
- 7 if(((m+offset).ge.u).or.((m+offset).lt.l))then
- goto 8
- else
- if(offset.lt.0)then
- lower = l
- check = m + offset
- upper = check
- else
- lower = m + offset + 1
- check = lower
- upper = u
- end if
- call ehg106(lower,upper,check,1,x(1,k),pi,n)
- if(x(pi(m + offset),k).eq.x(pi(m+offset+1),k))then
- offset = (-offset)
- if(offset.ge.0)then
- offset = offset + 1
- end if
- goto 7
- else
- m = m + offset
- goto 8
- end if
- end if
-
- 8 if(v(c(1,p),k).eq.x(pi(m),k))then
- leaf=.true.
- else
- leaf=(v(c(vc,p),k).eq.x(pi(m),k))
- end if
- end if
- if(leaf)then
- a(p)=0
- else
- a(p)=k
- xi(p)=x(pi(m),k)
-c left son
- nc=nc+1
- lo(p)=nc
- lo(nc)=l
- hi(nc)=m
-c right son
- nc=nc+1
- hi(p)=nc
- lo(nc)=m+1
- hi(nc)=u
- call ehg125(p,nv,v,vhit,nvmax,d,k,xi(p),2**(k-1),2**(d-k),c(
- $1,p),c(1,lo(p)),c(1,hi(p)))
- end if
- p=p+1
- l=lo(p)
- u=hi(p)
- goto 3
-c bottom of while loop
- 4 return
- end
//GO.SYSIN DD rbuild.f
echo spread.f 1>&2
- subroutine ehg129(l,u,d,x,pi,n,sigma)
- integer d,execnt,i,k,l,n,u
- integer pi(n)
- real machin,alpha,beta,t
- real sigma(d),x(n,d)
- real r1mach
- external r1mach
- save machin,execnt
- data execnt /0/
-c MachInf -> machin
- execnt=execnt+1
- if(execnt.eq.1)then
- machin=r1mach(2)
- end if
- do 3 k=1,d
- alpha=machin
- beta=-machin
- do 4 i=l,u
- t=x(pi(i),k)
- alpha=min(alpha,x(pi(i),k))
- beta=max(beta,t)
- 4 continue
- sigma(k)=beta-alpha
- 3 continue
- return
- end
//GO.SYSIN DD spread.f
echo vleaf.f 1>&2
- subroutine ehg137(z,kappa,leaf,nleaf,d,nv,nvmax,ncmax,vc,a,xi,lo,h
- $i,c,v)
- integer d,execnt,nc,ncmax,nleaf,p,stackt
- integer a(ncmax),hi(ncmax),leaf(256),lo(ncmax),pstack(20)
- real xi(ncmax),z(d)
- external ehg182
- save execnt
- data execnt /0/
-c stacktop -> stackt
- execnt=execnt+1
-c find leaf cells affected by $z$
- stackt=0
- p=1
- nleaf=0
-c top of while loop
- 3 if(.not.(0.lt.p))goto 4
- if(a(p).eq.0)then
-c leaf
- nleaf=nleaf+1
- leaf(nleaf)=p
-c Pop
- if(stackt.ge.1)then
- p=pstack(stackt)
- else
- p=0
- end if
- stackt=max(0,stackt-1)
- else
- if(z(a(p)).eq.xi(p))then
-c Push
- stackt=stackt+1
- if(.not.(stackt.le.20))then
- call ehg182(187)
- end if
- pstack(stackt)=hi(p)
- p=lo(p)
- else
- if(z(a(p)).le.xi(p))then
- p=lo(p)
- else
- p=hi(p)
- end if
- end if
- end if
- goto 3
-c bottom of while loop
- 4 if(.not.(nleaf.le.256))then
- call ehg182(185)
- end if
- return
- end
//GO.SYSIN DD vleaf.f
echo ehg182.f 1>&2
- subroutine ehg182(i)
- integer i
- if(i.eq.100) print *,'wrong version number in lowesd. Probably ty
- +po in caller.'
- if(i.eq.101) print *,'d>dMAX in ehg131. Need to recompile with in
- +creased dimensions.'
- if(i.eq.102) print *,'liv too small. (Discovered by lowesd)'
- if(i.eq.103) print *,'lv too small. (Discovered by lowesd)'
- if(i.eq.104) print *,'alpha too small. fewer data values than deg
- +rees of freedom.'
- if(i.eq.105) print *,'k>d2MAX in ehg136. Need to recompile with i
- +ncreased dimensions.'
- if(i.eq.106) print *,'lwork too small'
- if(i.eq.107) print *,'invalid value for kernel'
- if(i.eq.108) print *,'invalid value for ideg'
- if(i.eq.109) print *,'lowstt only applies when kernel=1.'
- if(i.eq.110) print *,'not enough extra workspace for robustness ca
- +lculation'
- if(i.eq.120) print *,'zero-width neighborhood. make alpha bigger'
- if(i.eq.121) print *,'all data on boundary of neighborhood. make a
- +lpha bigger'
- if(i.eq.122) print *,'extrapolation not allowed with blending'
- if(i.eq.123) print *,'ihat=1 (diag L) in l2fit only makes sense if
- + z=x (eval=data).'
- if(i.eq.171) print *,'lowesd must be called first.'
- if(i.eq.172) print *,'lowesf must not come between lowesb and lowe
- +se, lowesr, or lowesl.'
- if(i.eq.173) print *,'lowesb must come before lowese, lowesr, or l
- +owesl.'
- if(i.eq.174) print *,'lowesb need not be called twice.'
- if(i.eq.180) print *,'nv>nvmax in cpvert.'
- if(i.eq.181) print *,'nt>20 in eval.'
- if(i.eq.182) print *,'svddc failed in l2fit.'
- if(i.eq.183) print *,'didnt find edge in vleaf.'
- if(i.eq.184) print *,'zero-width cell found in vleaf.'
- if(i.eq.185) print *,'trouble descending to leaf in vleaf.'
- if(i.eq.186) print *,'insufficient workspace for lowesf.'
- if(i.eq.187) print *,'insufficient stack space'
- if(i.eq.188) print *,'lv too small for computing explicit L'
- if(i.eq.191) print *,'computed trace L was negative; something is
- +wrong!'
- if(i.eq.192) print *,'computed delta was negative; something is wr
- +ong!'
- if(i.eq.193) print *,'workspace in loread appears to be corrupted'
- if(i.eq.194) print *,'trouble in l2fit/l2tr'
- if(i.eq.195) print *,'only constant, linear, or quadratic local mo
- +dels allowed'
- if(i.eq.196) print *,'degree must be at least 1 for vertex influen
- +ce matrix'
- if(i.eq.999) print *,'not yet implemented'
- print *,'Assert failed, error code ',i
- stop
- end
- subroutine ehg183(s,i,n,inc)
- character*(*) s
- integer n, inc, i(inc,n), j
- print *,s,(i(1,j),j=1,n)
- end
- subroutine ehg184(s,x,n,inc)
- character*(*) s
- integer n, inc, j
- real x(inc,n)
- print *,s,(x(1,j),j=1,n)
- end
//GO.SYSIN DD ehg182.f
echo ehg106.f 1>&2
- subroutine ehg106(il,ir,k,nk,p,pi,n)
- integer execnt,i,ii,il,ir,j,k,l,n,nk,r
- integer pi(n)
- real t
- real p(nk,n)
- save execnt
- data execnt /0/
- execnt=execnt+1
-c find the $k$-th smallest of $n$ elements
-c Floyd+Rivest, CACM Mar '75, Algorithm 489
- l=il
- r=ir
-c top of while loop
- 3 if(.not.(l.lt.r))goto 4
-c to avoid recursion, sophisticated partition deleted
-c partition $x sub {l..r}$ about $t$
- t=p(1,pi(k))
- i=l
- j=r
- ii=pi(l)
- pi(l)=pi(k)
- pi(k)=ii
- if(t.lt.p(1,pi(r)))then
- ii=pi(l)
- pi(l)=pi(r)
- pi(r)=ii
- end if
-c top of while loop
- 5 if(.not.(i.lt.j))goto 6
- ii=pi(i)
- pi(i)=pi(j)
- pi(j)=ii
- i=i+1
- j=j-1
-c top of while loop
- 7 if(.not.(p(1,pi(i)).lt.t))goto 8
- i=i+1
- goto 7
-c bottom of while loop
- 8 continue
-c top of while loop
- 9 if(.not.(t.lt.p(1,pi(j))))goto 10
- j=j-1
- goto 9
-c bottom of while loop
- 10 goto 5
-c bottom of while loop
- 6 if(p(1,pi(l)).eq.t)then
- ii=pi(l)
- pi(l)=pi(j)
- pi(j)=ii
- else
- j=j+1
- ii=pi(r)
- pi(r)=pi(j)
- pi(j)=ii
- end if
- if(j.le.k)then
- l=j+1
- end if
- if(k.le.j)then
- r=j-1
- end if
- goto 3
-c bottom of while loop
- 4 return
- end
//GO.SYSIN DD ehg106.f
echo ehg140.f 1>&2
- subroutine ehg140(iw,i,j)
- integer execnt,i,j
- integer iw(i)
- save execnt
- data execnt /0/
- execnt=execnt+1
- iw(i)=j
- return
- end
//GO.SYSIN DD ehg140.f
echo ehg127.f 1>&2
- subroutine ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w,r
- $cond,sing,sigma,u,e,gamma,qraux,work,tol,dd,tdeg,cdeg,s)
- integer column,d,dd,execnt,i,i3,i9,info,inorm2,j,jj,jpvt,k,kernel,
- $n,nf,od,sing,tdeg
- integer cdeg(8),psi(n)
- real machep,f,i1,i10,i2,i4,i5,i6,i7,i8,rcond,rho,scal,tol
- real g(15),sigma(15),u(15,15),e(15,15),b(nf,k),colnor(15),dist(n),
- $eta(nf),gamma(15),q(d),qraux(15),rw(n),s(0:od),w(nf),work(15),x(n,
- $d),y(n)
- external ehg106,ehg182,ehg184,sqrdc,sqrsl,ssvdc
- integer isamax
- external isamax
- real r1mach
- external r1mach
- real sdot
- external sdot
- save machep,execnt
- data execnt /0/
-c colnorm -> colnor
-c E -> g
-c MachEps -> machep
-c V -> e
-c X -> b
- execnt=execnt+1
- if(execnt.eq.1)then
- machep=r1mach(4)
- end if
-c sort by distance
- do 3 i3=1,n
- dist(i3)=0
- 3 continue
- do 4 j=1,dd
- i4=q(j)
- do 5 i3=1,n
- dist(i3)=dist(i3)+(x(i3,j)-i4)**2
- 5 continue
- 4 continue
- call ehg106(1,n,nf,1,dist,psi,n)
- rho=dist(psi(nf))*max(1.,f)
- if(.not.(0.lt.rho))then
- call ehg182(120)
- end if
-c compute neighborhood weights
- if(kernel.eq.2)then
- do 6 i=1,nf
- if(dist(psi(i)).lt.rho)then
- i1=sqrt(rw(psi(i)))
- else
- i1=0
- end if
- w(i)=i1
- 6 continue
- else
- do 7 i3=1,nf
- w(i3)=sqrt(dist(psi(i3))/rho)
- 7 continue
- do 8 i3=1,nf
- w(i3)=sqrt(rw(psi(i3))*(1-w(i3)**3)**3)
- 8 continue
- end if
- if(abs(w(isamax(nf,w,1))).eq.0)then
- call ehg184('at ',q,dd,1)
- call ehg184('radius ',rho,1,1)
- if(.not..false.)then
- call ehg182(121)
- end if
- end if
-c fill design matrix
- column=1
- do 9 i3=1,nf
- b(i3,column)=w(i3)
- 9 continue
- if(tdeg.ge.1)then
- do 10 j=1,d
- if(cdeg(j).ge.1)then
- column=column+1
- i5=q(j)
- do 11 i3=1,nf
- b(i3,column)=w(i3)*(x(psi(i3),j)-i5)
- 11 continue
- end if
- 10 continue
- end if
- if(tdeg.ge.2)then
- do 12 j=1,d
- if(cdeg(j).ge.1)then
- if(cdeg(j).ge.2)then
- column=column+1
- i6=q(j)
- do 13 i3=1,nf
- b(i3,column)=w(i3)*(x(psi(i3),j)-i6)**2
- 13 continue
- end if
- do 14 jj=j+1,d
- if(cdeg(jj).ge.1)then
- column=column+1
- i7=q(j)
- i8=q(jj)
- do 15 i3=1,nf
- b(i3,column)=w(i3)*(x(psi(i3),j)-i7)*(x(psi(i3),
- $jj)-i8)
- 15 continue
- end if
- 14 continue
- end if
- 12 continue
- k=column
- end if
- do 16 i3=1,nf
- eta(i3)=w(i3)*y(psi(i3))
- 16 continue
-c equilibrate columns
- do 17 j=1,k
- scal=0
- do 18 inorm2=1,nf
- scal=scal+b(inorm2,j)**2
- 18 continue
- scal=sqrt(scal)
- if(0.lt.scal)then
- do 19 i3=1,nf
- b(i3,j)=b(i3,j)/scal
- 19 continue
- colnor(j)=scal
- else
- colnor(j)=1
- end if
- 17 continue
-c singular value decomposition
- call sqrdc(b,nf,nf,k,qraux,jpvt,work,0)
- call sqrsl(b,nf,nf,k,qraux,eta,work,eta,eta,work,work,1000,info)
- do 20 i9=1,k
- do 21 i3=1,k
- u(i3,i9)=0
- 21 continue
- 20 continue
- do 22 i=1,k
- do 23 j=i,k
- u(i,j)=b(i,j)
- 23 continue
- 22 continue
- call ssvdc(u,15,k,k,sigma,g,u,15,e,15,work,21,info)
- if(.not.(info.eq.0))then
- call ehg182(182)
- end if
- tol=sigma(1)*(100*machep)
- rcond=min(rcond,sigma(k)/sigma(1))
- if(sigma(k).le.tol)then
- sing=sing+1
- if(sing.eq.1)then
- call ehg184('Warning. pseudoinverse used at',q,d,1)
- call ehg184('neighborhood radius',sqrt(rho),1,1)
- call ehg184('reciprocal condition number ',rcond,1,1)
- else
- if(sing.eq.2)then
- call ehg184('There are other near singularities as well.'
- $,rho,1,1)
- end if
- end if
- end if
-c compensate for equilibration
- do 24 j=1,k
- i10=colnor(j)
- do 25 i3=1,k
- e(j,i3)=e(j,i3)/i10
- 25 continue
- 24 continue
-c solve least squares problem
- do 26 j=1,k
- if(tol.lt.sigma(j))then
- i2=sdot(k,u(1,j),1,eta,1)/sigma(j)
- else
- i2=0.
- end if
- gamma(j)=i2
- 26 continue
- do 27 j=0,od
-c bug fix 2006-07-04 for k=1, od>1. (thanks btyner@gmail.com)
- if(j.lt.k)then
- s(j)=sdot(k,e(j+1,1),15,gamma,1)
- else
- s(j)=0.
- end if
- 27 continue
- return
- end
//GO.SYSIN DD ehg127.f
echo losave.f 1>&2
- subroutine losave(iunit,iv,liv,lv,v)
- integer execnt,iunit,liv,lv
- integer iv(liv)
- real v(lv)
- external ehg167
- save execnt
- data execnt /0/
- execnt=execnt+1
- call ehg167(iunit,iv(2),iv(4),iv(5),iv(6),iv(14),v(iv(11)),iv(iv(7
- $)),v(iv(12)),v(iv(13)))
- return
- end
//GO.SYSIN DD losave.f
echo ehg167.f 1>&2
- subroutine ehg167(iunit,d,vc,nc,nv,nvmax,v,a,xi,vval)
- integer iunit,d,vc,nc,nv,a(nc),magic,i,j
- real v(nvmax,d),xi(nc),vval(0:d,nv)
- write(iunit,*)d,nc,nv
- do 10 i=1,d
-10 write(iunit,*)v(1,i),v(vc,i)
- j = 0
- do 20 i=1,nc
- if(a(i).ne.0)then
- write(iunit,*)a(i),xi(i)
- else
- write(iunit,*)a(i),j
- end if
-20 continue
- do 30 i=1,nv
-30 write(iunit,*)(vval(j,i),j=0,d)
- end
//GO.SYSIN DD ehg167.f
echo lohead.f 1>&2
- subroutine lohead(iunit,d,vc,nc,nv)
- integer iunit,d,vc,nc,nv
- read(iunit,*)d,nc,nv
- vc = 2**d
- end
//GO.SYSIN DD lohead.f
echo loread.f 1>&2
- subroutine loread(iunit,d,vc,nc,nv,iv,liv,lv,v)
- integer bound,d,execnt,iunit,liv,lv,nc,nv,vc
- integer iv(liv)
- real v(lv)
- external ehg168,ehg169,ehg182
- save execnt
- data execnt /0/
- execnt=execnt+1
- iv(28)=173
- iv(2)=d
- iv(4)=vc
- iv(14)=nv
- iv(17)=nc
- iv(7)=50
- iv(8)=iv(7)+nc
- iv(9)=iv(8)+vc*nc
- iv(10)=iv(9)+nc
- bound=iv(10)+nc
- if(.not.(bound-1.le.liv))then
- call ehg182(102)
- end if
- iv(11)=50
- iv(13)=iv(11)+nv*d
- iv(12)=iv(13)+(d+1)*nv
- bound=iv(12)+nc
- if(.not.(bound-1.le.lv))then
- call ehg182(103)
- end if
- call ehg168(iunit,d,vc,nc,nv,nv,v(iv(11)),iv(iv(7)),v(iv(12)),v(iv
- $(13)))
- call ehg169(d,vc,nc,nc,nv,nv,v(iv(11)),iv(iv(7)),v(iv(12)),iv(iv(8
- $)),iv(iv(9)),iv(iv(10)))
- return
- end
//GO.SYSIN DD loread.f
echo ehg168.f 1>&2
- subroutine ehg168(iunit,d,vc,nc,nv,nvmax,v,a,xi,vval)
- integer iunit,d,vc,nc,nv,a(nc),magic,i,j
- real v(nvmax,d),xi(nc),vval(0:d,nv)
- do 10 i=1,d
-10 read(iunit,*)v(1,i),v(vc,i)
- do 20 i=1,nc
-20 read(iunit,*)a(i),xi(i)
- do 30 i=1,nv
-30 read(iunit,*)(vval(j,i),j=0,d)
- end
//GO.SYSIN DD ehg168.f
echo ehg169.f 1>&2
- subroutine ehg169(d,vc,nc,ncmax,nv,nvmax,v,a,xi,c,hi,lo)
- integer d,execnt,i,j,k,mc,mv,nc,ncmax,nv,nvmax,p,vc
- integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax),novhit(1)
- real v(nvmax,d),xi(ncmax)
- external ehg125,ehg182
- integer ifloor
- external ifloor
- save execnt
- data execnt /0/
- execnt=execnt+1
-c as in bbox
-c remaining vertices
- do 3 i=2,vc-1
- j=i-1
- do 4 k=1,d
- v(i,k)=v(1+mod(j,2)*(vc-1),k)
- j=ifloor(float(j)/2.)
- 4 continue
- 3 continue
-c as in ehg131
- mc=1
- mv=vc
- novhit(1)=-1
- do 5 j=1,vc
- c(j,mc)=j
- 5 continue
-c as in rbuild
- p=1
-c top of while loop
- 6 if(.not.(p.le.nc))goto 7
- if(a(p).ne.0)then
- k=a(p)
-c left son
- mc=mc+1
- lo(p)=mc
-c right son
- mc=mc+1
- hi(p)=mc
- call ehg125(p,mv,v,novhit,nvmax,d,k,xi(p),2**(k-1),2**(d-k),
- $c(1,p),c(1,lo(p)),c(1,hi(p)))
- end if
- p=p+1
- goto 6
-c bottom of while loop
- 7 if(.not.(mc.eq.nc))then
- call ehg182(193)
- end if
- if(.not.(mv.eq.nv))then
- call ehg182(193)
- end if
- return
- end
//GO.SYSIN DD ehg169.f
echo ehg170.f 1>&2
- subroutine ehg170(k,d,vc,nv,nvmax,nc,ncmax,a,c,hi,lo,v,vval,xi)
- integer d,execnt,i,j,nc,ncmax,nv,nvmax,vc
- integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax)
- real v(nvmax,d),vval(0:d,nvmax),xi(ncmax)
- save execnt
- data execnt /0/
- execnt=execnt+1
- write(k,*)' real function loeval(z)'
- write(k,50)d
- write(k,*)' integer d,vc,nv,nc'
- write(k,51)nc,vc,nc
- write(k,52)nc,nc
- write(k,53)nv,d
- write(k,54)d,nv
- write(k,55)nc
- write(k,56)
- write(k,57)d,vc,nv,nc
-50 format(' real z(',i2,')')
-51 format(' integer a(',i5,'), c(',i3,',',i5,')')
-52 format(' integer hi(',i5,'), lo(',i5,')')
-53 format(' real v(',i5,',',i2,')')
-54 format(' real vval(0:',i2,',',i5,')')
-55 format(' real xi(',i5,')')
-56 format(' real ehg128')
-57 format(' data d,vc,nv,nc /',i2,',',i3,',',i5,',',i5,'/')
- do 3 i=1,nc
- write(k,58)i,a(i)
-58 format(' data a(',i5,') /',i5,'/')
- if(a(i).ne.0)then
- write(k,59)i,i,i,hi(i),lo(i),xi(i)
-59 format(' data hi(',i5,'),lo(',i5,'),xi(',i5,') /',
- $ i5,',',i5,',',1pe15.6,'/')
- end if
- do 4 j=1,vc
- write(k,60)j,i,c(j,i)
-60 format(' data c(',i3,',',i5,') /',i5,'/')
- 4 continue
- 3 continue
- do 5 i=1,nv
- write(k,61)i,vval(0,i)
-61 format(' data vval(0,',i5,') /',1pe15.6,'/')
- do 6 j=1,d
- write(k,62)i,j,v(i,j)
-62 format(' data v(',i5,',',i2,') /',1pe15.6,'/')
- write(k,63)j,i,vval(j,i)
-63 format(' data vval(',i2,',',i5,') /',1pe15.6,'/')
- 6 continue
- 5 continue
- write(k,*)' loeval=ehg128(z,d,nc,vc,a,xi,lo,hi,c,v,nv,vval)'
- write(k,*)' end'
- return
- end
//GO.SYSIN DD ehg170.f
echo lofort.f 1>&2
- subroutine lofort(iunit,iv,liv,lv,wv)
- integer execnt,iunit
- integer iv(*)
- real wv(*)
- external ehg170
- save execnt
- data execnt /0/
- execnt=execnt+1
- call ehg170(iunit,iv(2),iv(4),iv(6),iv(14),iv(5),iv(17),iv(iv(7)),
- $iv(iv(8)),iv(iv(9)),iv(iv(10)),wv(iv(11)),wv(iv(13)),wv(iv(12)))
- return
- end
//GO.SYSIN DD lofort.f
echo ehg141.f 1>&2
- subroutine ehg141(trl,n,deg,k,d,nsing,dk,delta1,delta2)
- integer d,deg,dk,k,n,nsing
- external ehg176
- real ehg176
- real corx,delta1,delta2,trl,z
- real c(48), c1, c2, c3, c4
-c coef, d, deg, del
- data c /
- $ .2971620,.3802660,.5886043
- $,.4263766,.3346498,.6271053
- $,.5241198,.3484836,.6687687
- $,.6338795,.4076457,.7207693
- $,.1611761,.3091323,.4401023
- $,.2939609,.3580278,.5555741
- $,.3972390,.4171278,.6293196
- $,.4675173,.4699070,.6674802
- $,.2848308,.2254512,.2914126
- $,.5393624,.2517230,.3898970
- $,.7603231,.2969113,.4740130
- $,.9664956,.3629838,.5348889
- $,.2075670,.2822574,.2369957
- $,.3911566,.2981154,.3623232
- $,.5508869,.3501989,.4371032
- $,.7002667,.4291632,.4930370
- $ /
- if(deg.eq.0) dk=1
- if(deg.eq.1) dk=d+1
- if(deg.eq.2) dk=float((d+2)*(d+1))/2.
- corx=sqrt(k/float(n))
- z=(sqrt(k/trl)-corx)/(1-corx)
- if(nsing .eq. 0 .and. 1 .lt. z)
- $ call ehg184('Chernobyl! trL<k',trl,1,1)
- if(z .lt. 0) call ehg184('Chernobyl! trL>n',trl,1,1)
- z=min(1.0,max(0.0,z))
- c4=exp(ehg176(z))
- i=1+3*(min(d,4)-1+4*(deg-1))
- if(d.le.4)then
- c1=c(i)
- c2=c(i+1)
- c3=c(i+2)
- else
- c1=c(i)+(d-4)*(c(i)-c(i-3))
- c2=c(i+1)+(d-4)*(c(i+1)-c(i-2))
- c3=c(i+2)+(d-4)*(c(i+2)-c(i-1))
- endif
- delta1=n-trl*exp(c1*z**c2*(1-z)**c3*c4)
- i=i+24
- if(d.le.4)then
- c1=c(i)
- c2=c(i+1)
- c3=c(i+2)
- else
- c1=c(i)+(d-4)*(c(i)-c(i-3))
- c2=c(i+1)+(d-4)*(c(i+1)-c(i-2))
- c3=c(i+2)+(d-4)*(c(i+2)-c(i-1))
- endif
- delta2=n-trl*exp(c1*z**c2*(1-z)**c3*c4)
- return
- end
//GO.SYSIN DD ehg141.f
echo ehg176.f 1>&2
- real function ehg176(z)
- real z(*)
- integer d,vc,nv,nc
- integer a(17), c(2,17)
- integer hi(17), lo(17)
- real v(10,1)
- real vval(0:1,10)
- real xi(17)
- real ehg128
- data d,vc,nv,nc /1,2,10,17/
- data a(1) /1/
- data hi(1),lo(1),xi(1) /3,2,0.3705/
- data c(1,1) /1/
- data c(2,1) /2/
- data a(2) /1/
- data hi(2),lo(2),xi(2) /5,4,0.2017/
- data c(1,2) /1/
- data c(2,2) /3/
- data a(3) /1/
- data hi(3),lo(3),xi(3) /7,6,0.5591/
- data c(1,3) /3/
- data c(2,3) /2/
- data a(4) /1/
- data hi(4),lo(4),xi(4) /9,8,0.1204/
- data c(1,4) /1/
- data c(2,4) /4/
- data a(5) /1/
- data hi(5),lo(5),xi(5) /11,10,0.2815/
- data c(1,5) /4/
- data c(2,5) /3/
- data a(6) /1/
- data hi(6),lo(6),xi(6) /13,12,0.4536/
- data c(1,6) /3/
- data c(2,6) /5/
- data a(7) /1/
- data hi(7),lo(7),xi(7) /15,14,0.7132/
- data c(1,7) /5/
- data c(2,7) /2/
- data a(8) /0/
- data c(1,8) /1/
- data c(2,8) /6/
- data a(9) /0/
- data c(1,9) /6/
- data c(2,9) /4/
- data a(10) /0/
- data c(1,10) /4/
- data c(2,10) /7/
- data a(11) /0/
- data c(1,11) /7/
- data c(2,11) /3/
- data a(12) /0/
- data c(1,12) /3/
- data c(2,12) /8/
- data a(13) /0/
- data c(1,13) /8/
- data c(2,13) /5/
- data a(14) /0/
- data c(1,14) /5/
- data c(2,14) /9/
- data a(15) /1/
- data hi(15),lo(15),xi(15) /17,16,0.8751/
- data c(1,15) /9/
- data c(2,15) /2/
- data a(16) /0/
- data c(1,16) /9/
- data c(2,16) /10/
- data a(17) /0/
- data c(1,17) /10/
- data c(2,17) /2/
- data vval(0,1) /-9.0572E-2/
- data v(1,1) /-5.E-3/
- data vval(1,1) /4.4844/
- data vval(0,2) /-1.0856E-2/
- data v(2,1) /1.005/
- data vval(1,2) /-0.7736/
- data vval(0,3) /-5.3718E-2/
- data v(3,1) /0.3705/
- data vval(1,3) /-0.3495/
- data vval(0,4) /2.6152E-2/
- data v(4,1) /0.2017/
- data vval(1,4) /-0.7286/
- data vval(0,5) /-5.8387E-2/
- data v(5,1) /0.5591/
- data vval(1,5) /0.1611/
- data vval(0,6) /9.5807E-2/
- data v(6,1) /0.1204/
- data vval(1,6) /-0.7978/
- data vval(0,7) /-3.1926E-2/
- data v(7,1) /0.2815/
- data vval(1,7) /-0.4457/
- data vval(0,8) /-6.4170E-2/
- data v(8,1) /0.4536/
- data vval(1,8) /3.2813E-2/
- data vval(0,9) /-2.0636E-2/
- data v(9,1) /0.7132/
- data vval(1,9) /0.3350/
- data vval(0,10) /4.0172E-2/
- data v(10,1) /0.8751/
- data vval(1,10) /-4.1032E-2/
- ehg176=ehg128(z,d,nc,vc,a,xi,lo,hi,c,v,nv,vval)
- end
//GO.SYSIN DD ehg176.f
echo ehg196.f 1>&2
- subroutine ehg196(tau,d,f,trl)
- integer d,dka,dkb,execnt,tau
- real alpha,f,trl,trla,trlb
- external ehg197
- save execnt
- data execnt /0/
- execnt=execnt+1
- call ehg197(1,tau,d,f,dka,trla)
- call ehg197(2,tau,d,f,dkb,trlb)
- alpha=float(tau-dka)/float(dkb-dka)
- trl=(1-alpha)*trla+alpha*trlb
- return
- end
//GO.SYSIN DD ehg196.f
echo ehg197.f 1>&2
- subroutine ehg197(deg,tau,d,f,dk,trl)
- integer d,deg,dk,tau
- real trl, f
- dk = 0
- if(deg.eq.1) dk=d+1
- if(deg.eq.2) dk=float((d+2)*(d+1))/2.
- g1 = (-0.08125*d+0.13)*d+1.05
- trl = dk*(1+max(0.,(g1-f)/f))
- return
- end
//GO.SYSIN DD ehg197.f
echo lowesa.f 1>&2
- subroutine lowesa(trl,n,d,tau,nsing,delta1,delta2)
- integer d,dka,dkb,execnt,n,nsing,tau
- real alpha,d1a,d1b,d2a,d2b,delta1,delta2,trl
- external ehg141
- save execnt
- data execnt /0/
- execnt=execnt+1
- call ehg141(trl,n,1,tau,d,nsing,dka,d1a,d2a)
- call ehg141(trl,n,2,tau,d,nsing,dkb,d1b,d2b)
- alpha=float(tau-dka)/float(dkb-dka)
- delta1=(1-alpha)*d1a+alpha*d1b
- delta2=(1-alpha)*d2a+alpha*d2b
- return
- end
//GO.SYSIN DD lowesa.f
echo lowesc.f 1>&2
- subroutine lowesc(n,l,ll,trl,delta1,delta2)
- integer execnt,i,j,n
- real delta1,delta2,trl
- real l(n,n),ll(n,n)
- real sdot
- external sdot
- save execnt
- data execnt /0/
- execnt=execnt+1
-c compute $LL~=~(I-L)(I-L)'$
- do 3 i=1,n
- l(i,i)=l(i,i)-1
- 3 continue
- do 4 i=1,n
- do 5 j=1,i
- ll(i,j)=sdot(n,l(i,1),n,l(j,1),n)
- 5 continue
- 4 continue
- do 6 i=1,n
- do 7 j=i+1,n
- ll(i,j)=ll(j,i)
- 7 continue
- 6 continue
- do 8 i=1,n
- l(i,i)=l(i,i)+1
- 8 continue
-c accumulate first two traces
- trl=0
- delta1=0
- do 9 i=1,n
- trl=trl+l(i,i)
- delta1=delta1+ll(i,i)
- 9 continue
-c $delta sub 2 = "tr" LL sup 2$
- delta2=0
- do 10 i=1,n
- delta2=delta2+sdot(n,ll(i,1),n,ll(1,i),1)
- 10 continue
- return
- end
//GO.SYSIN DD lowesc.f
echo lowesp.f 1>&2
- subroutine lowesp(n,y,yhat,pwgts,rwgts,pi,ytilde)
- integer identi,execnt,i2,i3,i5,m,n
- integer pi(n)
- real c,i1,i4,mad
- real pwgts(n),rwgts(n),y(n),yhat(n),ytilde(n)
- external ehg106
- integer ifloor
- external ifloor
- save execnt
- data execnt /0/
-c Identity -> identi
- execnt=execnt+1
-c median absolute deviation
- do 3 i5=1,n
- ytilde(i5)=abs(y(i5)-yhat(i5))*sqrt(pwgts(i5))
- 3 continue
- do 4 identi=1,n
- pi(identi)=identi
- 4 continue
- m=ifloor(float(n)/2.)+1
- call ehg106(1,n,m,1,ytilde,pi,n)
- if((n-m)+1.lt.m)then
- call ehg106(1,m-1,m-1,1,ytilde,pi,n)
- mad=(ytilde(pi(m-1))+ytilde(pi(m)))/2
- else
- mad=ytilde(pi(m))
- end if
-c magic constant
- c=(6*mad)**2/5
- do 5 i5=1,n
- ytilde(i5)=1-((y(i5)-yhat(i5))**2*pwgts(i5))/c
- 5 continue
- do 6 i5=1,n
- ytilde(i5)=ytilde(i5)*sqrt(rwgts(i5))
- 6 continue
- if(n.le.0)then
- i4=0.
- else
- i3=n
- i1=ytilde(i3)
- do 7 i2=i3-1,1,-1
- i1=ytilde(i2)+i1
- 7 continue
- i4=i1
- end if
- c=n/i4
-c pseudovalues
- do 8 i5=1,n
- ytilde(i5)=yhat(i5)+(c*rwgts(i5))*(y(i5)-yhat(i5))
- 8 continue
- return
- end
//GO.SYSIN DD lowesp.f
echo lowesw.f 1>&2
- subroutine lowesw(res,n,rw,pi)
- integer identi,execnt,i,i1,n,nh
- integer pi(n)
- real cmad,rsmall
- real res(n),rw(n)
- external ehg106
- integer ifloor
- external ifloor
- real r1mach
- external r1mach
- save execnt
- data execnt /0/
-c Identity -> identi
- execnt=execnt+1
-c tranliterated from Devlin's ratfor
-c find median of absolute residuals
- do 3 i1=1,n
- rw(i1)=abs(res(i1))
- 3 continue
- do 4 identi=1,n
- pi(identi)=identi
- 4 continue
- nh=ifloor(float(n)/2.)+1
-c partial sort to find 6*mad
- call ehg106(1,n,nh,1,rw,pi,n)
- if((n-nh)+1.lt.nh)then
- call ehg106(1,nh-1,nh-1,1,rw,pi,n)
- cmad=3*(rw(pi(nh))+rw(pi(nh-1)))
- else
- cmad=6*rw(pi(nh))
- end if
- rsmall=r1mach(1)
- if(cmad.lt.rsmall)then
- do 5 i1=1,n
- rw(i1)=1
- 5 continue
- else
- do 6 i=1,n
- if(cmad*0.999.lt.rw(i))then
- rw(i)=0
- else
- if(cmad*0.001.lt.rw(i))then
- rw(i)=(1-(rw(i)/cmad)**2)**2
- else
- rw(i)=1
- end if
- end if
- 6 continue
- end if
- return
- end
//GO.SYSIN DD lowesw.f
echo lowesl.f 1>&2
- subroutine lowesl(iv,liv,lv,wv,m,z,l)
- integer execnt,m,n
- integer iv(*)
- real l(m,*),wv(*),z(m,1)
- external ehg182,ehg191
- save execnt
- data execnt /0/
- execnt=execnt+1
- if(.not.(iv(28).ne.172))then
- call ehg182(172)
- end if
- if(.not.(iv(28).eq.173))then
- call ehg182(173)
- end if
- if(.not.(iv(26).ne.iv(34)))then
- call ehg182(175)
- end if
- call ehg191(m,z,l,iv(2),iv(3),iv(19),iv(6),iv(17),iv(4),iv(iv(7)),
- $wv(iv(12)),iv(iv(10)),iv(iv(9)),iv(iv(8)),wv(iv(11)),iv(14),wv(iv(
- $24)),wv(iv(34)),iv(iv(25)))
- return
- end
//GO.SYSIN DD lowesl.f
echo ehg191.f 1>&2
- subroutine ehg191(m,z,l,d,n,nf,nv,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vv
- $al2,lf,lq)
- integer lq1,d,execnt,i,i1,i2,j,m,n,nc,ncmax,nf,nv,nvmax,p,vc
- integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax)
- real l(m,n),lf(0:d,nvmax,nf),v(nvmax,d),vval2(0:d,nvmax),xi(ncmax)
- $,z(m,d),zi(8)
- real ehg128
- external ehg128
- save execnt
- data execnt /0/
- execnt=execnt+1
- do 3 j=1,n
- do 4 i2=1,nv
- do 5 i1=0,d
- vval2(i1,i2)=0
- 5 continue
- 4 continue
- do 6 i=1,nv
-c linear search for i in Lq
- lq1=lq(i,1)
- lq(i,1)=j
- p=nf
-c top of while loop
- 7 if(.not.(lq(i,p).ne.j))goto 8
- p=p-1
- goto 7
-c bottom of while loop
- 8 lq(i,1)=lq1
- if(lq(i,p).eq.j)then
- do 9 i1=0,d
- vval2(i1,i)=lf(i1,i,p)
- 9 continue
- end if
- 6 continue
- do 10 i=1,m
- do 11 i1=1,d
- zi(i1)=z(i,i1)
- 11 continue
- l(i,j)=ehg128(zi,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vval2)
- 10 continue
- 3 continue
- return
- end
//GO.SYSIN DD ehg191.f
.
