1C 2C Copyright (C) 1998--2020 The R Core Team 3 4C The authors of this software are Cleveland, Grosse, and Shyu. 5C Copyright (c) 1989, 1992 by AT&T. 6C Permission to use, copy, modify, and distribute this software for any 7C purpose without fee is hereby granted, provided that this entire notice 8C is included in all copies of any software which is or includes a copy 9C or modification of this software and in all copies of the supporting 10C documentation for such software. 11C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 12C WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY 13C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 14C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 15 16C altered by B.D. Ripley to 17C 18C remove unused variables 19C make phi in ehg139 double precision to match calling sequence 20C pass integer not logical from C 21C 22C by M. Maechler, renaming ehg182() to loesswarn() 23 24C Note that loesswarn(errormsg_code) is in ./loessc.c 25 26 subroutine ehg126(d,n,vc,x,v,nvmax) 27 integer d,execnt,i,j,k,n,nvmax,vc 28 DOUBLE PRECISION machin,alpha,beta,mu,t 29 DOUBLE PRECISION v(nvmax,d), x(n,d) 30 31 DOUBLE PRECISION D1MACH 32 external D1MACH 33 save machin,execnt 34 data execnt /0/ 35c MachInf -> machin 36 execnt=execnt+1 37 if(execnt.eq.1)then 38c initialize d1mach(2) === DBL_MAX: 39 machin=D1MACH(2) 40 end if 41c fill in vertices for bounding box of $x$ 42c lower left, upper right 43 do 3 k=1,d 44 alpha=machin 45 beta=-machin 46 do 4 i=1,n 47 t=x(i,k) 48 alpha=min(alpha,t) 49 beta=max(beta,t) 50 4 continue 51c expand the box a little 52 mu=0.005D0*max(beta-alpha,1.d-10*max(DABS(alpha),DABS(beta))+ 53 + 1.d-30) 54 alpha=alpha-mu 55 beta=beta+mu 56 v(1,k)=alpha 57 v(vc,k)=beta 58 3 continue 59c remaining vertices 60 do 5 i=2,vc-1 61 j=i-1 62 do 6 k=1,d 63 v(i,k)=v(1+mod(j,2)*(vc-1),k) 64c Integer division would do here 65 j=INT(DBLE(j)/2.D0) 66 6 continue 67 5 continue 68 return 69 end 70 71 subroutine ehg125(p,nv,v,vhit,nvmax,d,k,t,r,s,f,l,u) 72 logical i1,i2,match 73 integer d,h,i,i3,j,k,m,mm,nv,nvmax,p,r,s 74 integer f(r,0:1,s),l(r,0:1,s),u(r,0:1,s),vhit(nvmax) 75 DOUBLE PRECISION t 76 DOUBLE PRECISION v(nvmax,d) 77 external loesswarn 78 h=nv 79 do 3 i=1,r 80 do 4 j=1,s 81 h=h+1 82 do 5 i3=1,d 83 v(h,i3)=v(f(i,0,j),i3) 84 5 continue 85 v(h,k)=t 86c check for redundant vertex 87 match=.false. 88 m=1 89c top of while loop 90 6 if(.not.match)then 91 i1=(m.le.nv) 92 else 93 i1=.false. 94 end if 95 if(.not.(i1))goto 7 96 match=(v(m,1).eq.v(h,1)) 97 mm=2 98c top of while loop 99 8 if(match)then 100 i2=(mm.le.d) 101 else 102 i2=.false. 103 end if 104 if(.not.(i2))goto 9 105 match=(v(m,mm).eq.v(h,mm)) 106 mm=mm+1 107 goto 8 108c bottom of while loop 109 9 m=m+1 110 goto 6 111c bottom of while loop 112 7 m=m-1 113 if(match)then 114 h=h-1 115 else 116 m=h 117 if(vhit(1).ge.0)then 118 vhit(m)=p 119 end if 120 end if 121 l(i,0,j)=f(i,0,j) 122 l(i,1,j)=m 123 u(i,0,j)=m 124 u(i,1,j)=f(i,1,j) 125 4 continue 126 3 continue 127 nv=h 128 if(.not.(nv.le.nvmax))then 129 call loesswarn(180) 130 end if 131 return 132 end 133 134 integer function ehg138(i,z,a,xi,lo,hi,ncmax) 135 logical i1 136 integer i,j,ncmax 137 integer a(ncmax),hi(ncmax),lo(ncmax) 138 DOUBLE PRECISION xi(ncmax),z(8) 139c descend tree until leaf or ambiguous 140 j=i 141c top of while loop 142 3 if(a(j).ne.0)then 143 i1=(z(a(j)).ne.xi(j)) 144 else 145 i1=.false. 146 end if 147 if(.not.(i1))goto 4 148 if(z(a(j)).le.xi(j))then 149 j=lo(j) 150 else 151 j=hi(j) 152 end if 153 goto 3 154c bottom of while loop 155 4 ehg138=j 156 return 157 end 158 159 subroutine ehg106(il,ir,k,nk,p,pi,n) 160 161c Partial sorting of p(1, il:ir) returning the sort indices pi() only 162c such that p(1, pi(k)) is correct 163 164c implicit none 165c Arguments 166c Input: 167 integer il,ir,k,nk,n 168 DOUBLE PRECISION p(nk,n) 169c using only p(1, pi(*)) 170c Output: 171 integer pi(n) 172 173c Variables 174 DOUBLE PRECISION t 175 integer i,ii,j,l,r 176 177c find the $k$-th smallest of $n$ elements 178c Floyd+Rivest, CACM Mar '75, Algorithm 489 179 l=il 180 r=ir 181c while (l < r ) 182 3 if(.not.(l.lt.r))goto 4 183c to avoid recursion, sophisticated partition deleted 184c partition $x sub {l..r}$ about $t$ 185 t=p(1,pi(k)) 186 i=l 187 j=r 188 ii=pi(l) 189 pi(l)=pi(k) 190 pi(k)=ii 191 if(t.lt.p(1,pi(r)))then 192 ii=pi(l) 193 pi(l)=pi(r) 194 pi(r)=ii 195 end if 196c top of while loop 197 5 if(.not.(i.lt.j))goto 6 198 ii=pi(i) 199 pi(i)=pi(j) 200 pi(j)=ii 201 i=i+1 202 j=j-1 203c top of while loop 204 7 if(.not.(p(1,pi(i)).lt.t))goto 8 205 i=i+1 206 goto 7 207c bottom of while loop 208 8 continue 209c top of while loop 210 9 if(.not.(t.lt.p(1,pi(j))))goto 10 211 j=j-1 212 goto 9 213c bottom of while loop 214 10 goto 5 215c bottom of while loop 216 6 if(p(1,pi(l)).eq.t)then 217 ii=pi(l) 218 pi(l)=pi(j) 219 pi(j)=ii 220 else 221 j=j+1 222 ii=pi(r) 223 pi(r)=pi(j) 224 pi(j)=ii 225 end if 226 if(j.le.k)then 227 l=j+1 228 end if 229 if(k.le.j)then 230 r=j-1 231 end if 232 goto 3 233c bottom of while loop 234 4 return 235 end 236 237 238 subroutine ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w, 239 + rcond,sing,sigma,u,e,dgamma,qraux,work,tol,dd,tdeg,cdeg,s) 240 integer column,d,dd,execnt,i,i3,i9,info,inorm2,j,jj,jpvt,k,kernel, 241 + n,nf,od,sing,tdeg 242 integer cdeg(8),psi(n) 243 double precision machep,f,i1,i10,i2,i4,i5,i6,i7,i8,rcond,rho,scal, 244 + tol 245 double precision g(15),sigma(15),u(15,15),e(15,15),b(nf,k), 246 + colnor(15),dist(n),eta(nf),dgamma(15),q(d),qraux(15),rw(n), 247 + s(0:od),w(nf),work(15),x(n,d),y(n) 248 249 integer idamax 250 double precision d1mach, ddot 251 252 external ehg106,loesswarn,ehg184,dqrdc,dqrsl,dsvdc 253 external idamax, d1mach, ddot 254 255 save machep,execnt 256 data execnt /0/ 257c colnorm -> colnor 258c E -> g 259c MachEps -> machep 260c V -> e 261c X -> b 262 execnt=execnt+1 263 if(execnt.eq.1)then 264c initialize d1mach(4) === 1 / DBL_EPSILON === 2^52 : 265 machep=d1mach(4) 266 end if 267c sort by distance 268 do 3 i3=1,n 269 dist(i3)=0 270 3 continue 271 do 4 j=1,dd 272 i4=q(j) 273 do 5 i3=1,n 274 dist(i3)=dist(i3)+(x(i3,j)-i4)**2 275 5 continue 276 4 continue 277 call ehg106(1,n,nf,1,dist,psi,n) 278 rho=dist(psi(nf))*max(1.d0,f) 279 if(rho .le. 0)then 280 call loesswarn(120) 281 end if 282c compute neighborhood weights 283 if(kernel.eq.2)then 284 do 6 i=1,nf 285 if(dist(psi(i)).lt.rho)then 286 i1=dsqrt(rw(psi(i))) 287 else 288 i1=0 289 end if 290 w(i)=i1 291 6 continue 292 else 293 do 7 i3=1,nf 294 w(i3)=dsqrt(dist(psi(i3))/rho) 295 7 continue 296 do 8 i3=1,nf 297 w(i3)=dsqrt(rw(psi(i3))*(1-w(i3)**3)**3) 298 8 continue 299 end if 300 if(dabs(w(idamax(nf,w,1))).eq.0)then 301 call ehg184('at ',q(1),dd,1) 302 call ehg184('radius ',rho,1,1) 303 if(.not..false.)then 304 call loesswarn(121) 305 end if 306 end if 307c fill design matrix 308 column=1 309 do 9 i3=1,nf 310 b(i3,column)=w(i3) 311 9 continue 312 if(tdeg.ge.1)then 313 do 10 j=1,d 314 if(cdeg(j).ge.1)then 315 column=column+1 316 i5=q(j) 317 do 11 i3=1,nf 318 b(i3,column)=w(i3)*(x(psi(i3),j)-i5) 319 11 continue 320 end if 321 10 continue 322 end if 323 if(tdeg.ge.2)then 324 do 12 j=1,d 325 if(cdeg(j).ge.1)then 326 if(cdeg(j).ge.2)then 327 column=column+1 328 i6=q(j) 329 do 13 i3=1,nf 330 b(i3,column)=w(i3)*(x(psi(i3),j)-i6)**2 331 13 continue 332 end if 333 do 14 jj=j+1,d 334 if(cdeg(jj).ge.1)then 335 column=column+1 336 i7=q(j) 337 i8=q(jj) 338 do 15 i3=1,nf 339 b(i3,column)=w(i3)*(x(psi(i3),j)-i7)*(x(psi(i3), 340 +jj)-i8) 341 15 continue 342 end if 343 14 continue 344 end if 345 12 continue 346 k=column 347 end if 348 do 16 i3=1,nf 349 eta(i3)=w(i3)*y(psi(i3)) 350 16 continue 351c equilibrate columns 352 do 17 j=1,k 353 scal=0 354 do 18 inorm2=1,nf 355 scal=scal+b(inorm2,j)**2 356 18 continue 357 scal=dsqrt(scal) 358 if(0.lt.scal)then 359 do 19 i3=1,nf 360 b(i3,j)=b(i3,j)/scal 361 19 continue 362 colnor(j)=scal 363 else 364 colnor(j)=1 365 end if 366 17 continue 367c singular value decomposition 368 call dqrdc(b,nf,nf,k,qraux,jpvt,work,0) 369 call dqrsl(b,nf,nf,k,qraux,eta,work,eta,eta,work,work,1000,info) 370 do 20 i9=1,k 371 do 21 i3=1,k 372 u(i3,i9)=0 373 21 continue 374 20 continue 375 do 22 i=1,k 376 do 23 j=i,k 377c FIXME: this has i = 3 vs bound 2 in a ggplot2 test 378 u(i,j)=b(i,j) 379 23 continue 380 22 continue 381 call dsvdc(u,15,k,k,sigma,g,u,15,e,15,work,21,info) 382 if(.not.(info.eq.0))then 383 call loesswarn(182) 384 end if 385 tol=sigma(1)*(100*machep) 386 rcond=min(rcond,sigma(k)/sigma(1)) 387 if(sigma(k).le.tol)then 388 sing=sing+1 389 if(sing.eq.1)then 390 call ehg184('pseudoinverse used at',q(1),d,1) 391 call ehg184('neighborhood radius',dsqrt(rho),1,1) 392 call ehg184('reciprocal condition number ',rcond,1,1) 393 else 394 if(sing.eq.2)then 395 call ehg184('There are other near singularities as well.' 396 +,rho,1,1) 397 end if 398 end if 399 end if 400c compensate for equilibration 401 do 24 j=1,k 402 i10=colnor(j) 403 do 25 i3=1,k 404 e(j,i3)=e(j,i3)/i10 405 25 continue 406 24 continue 407c solve least squares problem 408 do 26 j=1,k 409 if(tol.lt.sigma(j))then 410 i2=ddot(k,u(1,j),1,eta,1)/sigma(j) 411 else 412 i2=0.d0 413 end if 414 dgamma(j)=i2 415 26 continue 416 do 27 j=0,od 417c bug fix 2006-07-04 for k=1, od>1. (thanks btyner@gmail.com) 418 if(j.lt.k)then 419 s(j)=ddot(k,e(j+1,1),15,dgamma,1) 420 else 421 s(j)=0.d0 422 end if 423 27 continue 424 return 425 end 426 427 subroutine ehg131(x,y,rw,trl,diagl,kernel,k,n,d,nc,ncmax,vc,nv, 428 + nvmax,nf,f,a,c,hi,lo,pi,psi,v,vhit,vval,xi,dist,eta,b,ntol, 429 + fd,w,vval2,rcond,sing,dd,tdeg,cdeg,lq,lf,setlf) 430 logical setlf 431 integer identi,d,dd,i1,i2,j,k,kernel,n,nc,ncmax,nf,ntol,nv, 432 + nvmax,sing,tdeg,vc 433 integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),cdeg(8),hi(ncmax), 434 + lo(ncmax),pi(n),psi(n),vhit(nvmax) 435 double precision f,fd,rcond,trl 436 double precision lf(0:d,nvmax,nf),b(*),delta(8),diagl(n),dist(n), 437 + eta(nf),rw(n),v(nvmax,d),vval(0:d,nvmax),vval2(0:d,nvmax), 438 + w(nf),x(n,d),xi(ncmax),y(n) 439 external ehg126,loesswarn,ehg139,ehg124 440 double precision dnrm2 441 external dnrm2 442c Identity -> identi 443c X -> b 444 if(.not.(d.le.8))then 445 call loesswarn(101) 446 end if 447c build $k$-d tree 448 call ehg126(d,n,vc,x,v,nvmax) 449 nv=vc 450 nc=1 451 do 3 j=1,vc 452 c(j,nc)=j 453 vhit(j)=0 454 3 continue 455 do 4 i1=1,d 456 delta(i1)=v(vc,i1)-v(1,i1) 457 4 continue 458 fd=fd*dnrm2(d,delta,1) 459 do 5 identi=1,n 460 pi(identi)=identi 461 5 continue 462 call ehg124(1,n,d,n,nv,nc,ncmax,vc,x,pi,a,xi,lo,hi,c,v,vhit,nvmax, 463 +ntol,fd,dd) 464c smooth 465 if(trl.ne.0)then 466 do 6 i2=1,nv 467 do 7 i1=0,d 468 vval2(i1,i2)=0 469 7 continue 470 6 continue 471 end if 472 call ehg139(v,nvmax,nv,n,d,nf,f,x,pi,psi,y,rw,trl,kernel,k,dist, 473 + dist,eta,b,d,w,diagl,vval2,nc,vc,a,xi,lo,hi,c,rcond, 474 + sing,dd,tdeg,cdeg,lq,lf,setlf,vval) 475 return 476 end 477 478c called from lowese() only : 479 subroutine ehg133(d,vc,nvmax,ncmax, a,c,hi,lo, v,vval,xi,m,z,s) 480 integer d,vc,nvmax,ncmax, m 481 integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax) 482 double precision v(nvmax,d),vval(0:d,nvmax),xi(ncmax),z(m,d),s(m) 483c Var 484 double precision delta(8) 485 integer i,i1 486 487 double precision ehg128 488 external ehg128 489 490 do 3 i=1,m 491 do 4 i1=1,d 492 delta(i1)=z(i,i1) 493 4 continue 494 s(i)=ehg128(delta,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vval) 495 3 continue 496 return 497 end 498 499 subroutine ehg140(iw,i,j) 500 integer i,j 501 integer iw(i) 502 iw(i)=j 503 return 504 end 505 506 subroutine ehg141(trl,n,deg,k,d,nsing,dk,delta1,delta2) 507 double precision trl,delta1,delta2 508 integer n,deg,k,d,nsing,dk 509 510 double precision c(48), c1, c2, c3, c4, corx,z,zz(1) 511 integer i 512 513 external ehg176 514 double precision ehg176 515 516c coef, d, deg, del 517 data c / .2971620d0,.3802660d0,.5886043d0,.4263766d0,.3346498d0, 518 +.6271053d0,.5241198d0,.3484836d0,.6687687d0,.6338795d0,.4076457d0, 519 +.7207693d0,.1611761d0,.3091323d0,.4401023d0,.2939609d0,.3580278d0, 520 +.5555741d0,.3972390d0,.4171278d0,.6293196d0,.4675173d0,.4699070d0, 521 +.6674802d0,.2848308d0,.2254512d0,.2914126d0,.5393624d0,.2517230d0, 522 +.3898970d0,.7603231d0,.2969113d0,.4740130d0,.9664956d0,.3629838d0, 523 +.5348889d0,.2075670d0,.2822574d0,.2369957d0,.3911566d0,.2981154d0, 524 +.3623232d0,.5508869d0,.3501989d0,.4371032d0,.7002667d0,.4291632d0, 525 +.4930370d0 / 526 527 if(deg.eq.0) dk=1 528 if(deg.eq.1) dk=d+1 529 if(deg.eq.2) dk=int(dble((d+2)*(d+1))/2.d0) 530 corx=dsqrt(k/dble(n)) 531 z=(dsqrt(k/trl)-corx)/(1-corx) 532 if(nsing .eq. 0 .and. 1 .lt. z) then 533 call ehg184('Chernobyl! trL<k',trl,1,1) 534 else if(z .lt. 0) then 535 call ehg184('Chernobyl! trL>n',trl,1,1) 536 endif 537 z=min(1.0d0,max(0.0d0,z)) 538c R fix 539 zz(1)=z 540 c4=dexp(ehg176(zz)) 541 i=1+3*(min(d,4)-1+4*(deg-1)) 542 if(d.le.4)then 543 c1=c(i) 544 c2=c(i+1) 545 c3=c(i+2) 546 else 547 c1=c(i)+(d-4)*(c(i)-c(i-3)) 548 c2=c(i+1)+(d-4)*(c(i+1)-c(i-2)) 549 c3=c(i+2)+(d-4)*(c(i+2)-c(i-1)) 550 endif 551 delta1=n-trl*dexp(c1*z**c2*(1-z)**c3*c4) 552 i=i+24 553 if(d.le.4)then 554 c1=c(i) 555 c2=c(i+1) 556 c3=c(i+2) 557 else 558 c1=c(i)+(d-4)*(c(i)-c(i-3)) 559 c2=c(i+1)+(d-4)*(c(i+1)-c(i-2)) 560 c3=c(i+2)+(d-4)*(c(i+2)-c(i-1)) 561 endif 562 delta2=n-trl*dexp(c1*z**c2*(1-z)**c3*c4) 563 return 564 end 565 566 subroutine lowesc(n,l,ll,trl,delta1,delta2) 567 integer i,j,n 568 double precision delta1,delta2,trl 569 double precision l(n,n),ll(n,n) 570 double precision ddot 571 external ddot 572c compute $LL~=~(I-L)(I-L)'$ 573 do 3 i=1,n 574 l(i,i)=l(i,i)-1 575 3 continue 576 do 4 i=1,n 577 do 5 j=1,i 578 ll(i,j)=ddot(n,l(i,1),n,l(j,1),n) 579 5 continue 580 4 continue 581 do 6 i=1,n 582 do 7 j=i+1,n 583 ll(i,j)=ll(j,i) 584 7 continue 585 6 continue 586 do 8 i=1,n 587 l(i,i)=l(i,i)+1 588 8 continue 589c accumulate first two traces 590 trl=0 591 delta1=0 592 do 9 i=1,n 593 trl=trl+l(i,i) 594 delta1=delta1+ll(i,i) 595 9 continue 596c $delta sub 2 = "tr" LL sup 2$ 597 delta2=0 598 do 10 i=1,n 599 delta2=delta2+ddot(n,ll(i,1),n,ll(1,i),1) 600 10 continue 601 return 602 end 603 604 subroutine ehg169(d,vc,nc,ncmax,nv,nvmax,v,a,xi,c,hi,lo) 605 integer d,vc,nc,ncmax,nv,nvmax 606 integer a(ncmax), c(vc,ncmax), hi(ncmax), lo(ncmax) 607 DOUBLE PRECISION v(nvmax,d),xi(ncmax) 608 609 integer novhit(1),i,j,k,mc,mv,p 610 external ehg125,loesswarn 611 integer ifloor 612 external ifloor 613 614c as in bbox 615c remaining vertices 616 do 3 i=2,vc-1 617 j=i-1 618 do 4 k=1,d 619 v(i,k)=v(1+mod(j,2)*(vc-1),k) 620 j=ifloor(DBLE(j)/2.D0) 621 4 continue 622 3 continue 623c as in ehg131 624 mc=1 625 mv=vc 626 novhit(1)=-1 627 do 5 j=1,vc 628 c(j,mc)=j 629 5 continue 630c as in rbuild 631 p=1 632c top of while loop 633 6 if(.not.(p.le.nc))goto 7 634 if(a(p).ne.0)then 635 k=a(p) 636c left son 637 mc=mc+1 638 lo(p)=mc 639c right son 640 mc=mc+1 641 hi(p)=mc 642 call ehg125(p,mv,v,novhit,nvmax,d,k,xi(p),2**(k-1),2**(d-k), 643 + c(1,p),c(1,lo(p)),c(1,hi(p))) 644 end if 645 p=p+1 646 goto 6 647c bottom of while loop 648 7 if(.not.(mc.eq.nc))then 649 call loesswarn(193) 650 end if 651 if(.not.(mv.eq.nv))then 652 call loesswarn(193) 653 end if 654 return 655 end 656 657 DOUBLE PRECISION function ehg176(z) 658c 659 DOUBLE PRECISION z(*) 660c 661 integer d,vc,nv,nc 662 integer a(17), c(2,17) 663 integer hi(17), lo(17) 664 DOUBLE PRECISION v(10,1) 665 DOUBLE PRECISION vval(0:1,10) 666 DOUBLE PRECISION xi(17) 667 double precision ehg128 668 external ehg128 669 670 data d,vc,nv,nc /1,2,10,17/ 671 data a(1) /1/ 672 data hi(1),lo(1),xi(1) /3,2,0.3705D0/ 673 data c(1,1) /1/ 674 data c(2,1) /2/ 675 data a(2) /1/ 676 data hi(2),lo(2),xi(2) /5,4,0.2017D0/ 677 data c(1,2) /1/ 678 data c(2,2) /3/ 679 data a(3) /1/ 680 data hi(3),lo(3),xi(3) /7,6,0.5591D0/ 681 data c(1,3) /3/ 682 data c(2,3) /2/ 683 data a(4) /1/ 684 data hi(4),lo(4),xi(4) /9,8,0.1204D0/ 685 data c(1,4) /1/ 686 data c(2,4) /4/ 687 data a(5) /1/ 688 data hi(5),lo(5),xi(5) /11,10,0.2815D0/ 689 data c(1,5) /4/ 690 data c(2,5) /3/ 691 data a(6) /1/ 692 data hi(6),lo(6),xi(6) /13,12,0.4536D0/ 693 data c(1,6) /3/ 694 data c(2,6) /5/ 695 data a(7) /1/ 696 data hi(7),lo(7),xi(7) /15,14,0.7132D0/ 697 data c(1,7) /5/ 698 data c(2,7) /2/ 699 data a(8) /0/ 700 data c(1,8) /1/ 701 data c(2,8) /6/ 702 data a(9) /0/ 703 data c(1,9) /6/ 704 data c(2,9) /4/ 705 data a(10) /0/ 706 data c(1,10) /4/ 707 data c(2,10) /7/ 708 data a(11) /0/ 709 data c(1,11) /7/ 710 data c(2,11) /3/ 711 data a(12) /0/ 712 data c(1,12) /3/ 713 data c(2,12) /8/ 714 data a(13) /0/ 715 data c(1,13) /8/ 716 data c(2,13) /5/ 717 data a(14) /0/ 718 data c(1,14) /5/ 719 data c(2,14) /9/ 720 data a(15) /1/ 721 data hi(15),lo(15),xi(15) /17,16,0.8751D0/ 722 data c(1,15) /9/ 723 data c(2,15) /2/ 724 data a(16) /0/ 725 data c(1,16) /9/ 726 data c(2,16) /10/ 727 data a(17) /0/ 728 data c(1,17) /10/ 729 data c(2,17) /2/ 730 data vval(0,1) /-9.0572D-2/ 731 data v(1,1) /-5.D-3/ 732 data vval(1,1) /4.4844D0/ 733 data vval(0,2) /-1.0856D-2/ 734 data v(2,1) /1.005D0/ 735 data vval(1,2) /-0.7736D0/ 736 data vval(0,3) /-5.3718D-2/ 737 data v(3,1) /0.3705D0/ 738 data vval(1,3) /-0.3495D0/ 739 data vval(0,4) /2.6152D-2/ 740 data v(4,1) /0.2017D0/ 741 data vval(1,4) /-0.7286D0/ 742 data vval(0,5) /-5.8387D-2/ 743 data v(5,1) /0.5591D0/ 744 data vval(1,5) /0.1611D0/ 745 data vval(0,6) /9.5807D-2/ 746 data v(6,1) /0.1204D0/ 747 data vval(1,6) /-0.7978D0/ 748 data vval(0,7) /-3.1926D-2/ 749 data v(7,1) /0.2815D0/ 750 data vval(1,7) /-0.4457D0/ 751 data vval(0,8) /-6.4170D-2/ 752 data v(8,1) /0.4536D0/ 753 data vval(1,8) /3.2813D-2/ 754 data vval(0,9) /-2.0636D-2/ 755 data v(9,1) /0.7132D0/ 756 data vval(1,9) /0.3350D0/ 757 data vval(0,10) /4.0172D-2/ 758 data v(10,1) /0.8751D0/ 759 data vval(1,10) /-4.1032D-2/ 760 ehg176=ehg128(z,d,nc,vc,a,xi,lo,hi,c,v,nv,vval) 761 end 762 763 subroutine lowesa(trl,n,d,tau,nsing,delta1,delta2) 764 integer n,d,tau,nsing 765 double precision trl, delta1,delta2 766 767 integer dka,dkb 768 double precision alpha,d1a,d1b,d2a,d2b 769 external ehg141 770 771 call ehg141(trl,n,1,tau,d,nsing,dka,d1a,d2a) 772 call ehg141(trl,n,2,tau,d,nsing,dkb,d1b,d2b) 773 alpha=dble(tau-dka)/dble(dkb-dka) 774 delta1=(1-alpha)*d1a+alpha*d1b 775 delta2=(1-alpha)*d2a+alpha*d2b 776 return 777 end 778 779 subroutine ehg191(m,z,l,d,n,nf,nv,ncmax,vc,a,xi,lo,hi,c,v,nvmax, 780 + vval2,lf,lq) 781c Args 782 integer m,d,n,nf,nv,ncmax,nvmax,vc 783 double precision z(m,d), l(m,n), xi(ncmax), v(nvmax,d), 784 + vval2(0:d,nvmax), lf(0:d,nvmax,nf) 785 integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),lo(ncmax),hi(ncmax) 786c Var 787 integer lq1,i,i1,i2,j,p 788 double precision zi(8) 789 double precision ehg128 790 external ehg128 791 792 do 3 j=1,n 793 do 4 i2=1,nv 794 do 5 i1=0,d 795 vval2(i1,i2)=0 796 5 continue 797 4 continue 798 do 6 i=1,nv 799c linear search for i in Lq 800 lq1=lq(i,1) 801 lq(i,1)=j 802 p=nf 803c top of while loop 804 7 if(.not.(lq(i,p).ne.j))goto 8 805 p=p-1 806 goto 7 807c bottom of while loop 808 8 lq(i,1)=lq1 809 if(lq(i,p).eq.j)then 810 do 9 i1=0,d 811 vval2(i1,i)=lf(i1,i,p) 812 9 continue 813 end if 814 6 continue 815 do 10 i=1,m 816 do 11 i1=1,d 817 zi(i1)=z(i,i1) 818 11 continue 819 l(i,j)=ehg128(zi,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vval2) 820 10 continue 821 3 continue 822 return 823 end 824 825 subroutine ehg196(tau,d,f,trl) 826 integer d,dka,dkb,tau 827 double precision alpha,f,trl,trla,trlb 828 external ehg197 829 call ehg197(1,d,f,dka,trla) 830 call ehg197(2,d,f,dkb,trlb) 831 alpha=dble(tau-dka)/dble(dkb-dka) 832 trl=(1-alpha)*trla+alpha*trlb 833 return 834 end 835 836 subroutine ehg197(deg,d,f,dk,trl) 837 integer deg,d,dk 838 double precision f, trl 839 840 double precision g1 841 dk = 0 842 if(deg.eq.1) dk=d+1 843 if(deg.eq.2) dk=int(dble((d+2)*(d+1))/2.d0) 844 g1 = (-0.08125d0*d+0.13d0)*d+1.05d0 845 trl = dk*(1+max(0.d0,(g1-f)/f)) 846 return 847 end 848 849 subroutine ehg192(y,d,n,nf,nv,nvmax,vval,lf,lq) 850 integer d,i,i1,i2,j,n,nf,nv,nvmax 851 integer lq(nvmax,nf) 852 DOUBLE PRECISION i3 853 DOUBLE PRECISION lf(0:d,nvmax,nf),vval(0:d,nvmax),y(n) 854 855 do 3 i2=1,nv 856 do 4 i1=0,d 857 vval(i1,i2)=0 858 4 continue 859 3 continue 860 do 5 i=1,nv 861 do 6 j=1,nf 862 i3=y(lq(i,j)) 863 do 7 i1=0,d 864 vval(i1,i)=vval(i1,i)+i3*lf(i1,i,j) 865 7 continue 866 6 continue 867 5 continue 868 return 869 end 870 871 DOUBLE PRECISION function ehg128(z,d,ncmax,vc,a,xi,lo,hi,c,v, 872 + nvmax,vval) 873 874c implicit none 875c Args 876 integer d,ncmax,nvmax,vc 877 integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax) 878 DOUBLE PRECISION z(d),xi(ncmax),v(nvmax,d), vval(0:d,nvmax) 879c Vars 880 logical i2,i3,i4,i5,i6,i7,i8,i9,i10 881 integer i,i1,i11,i12,ig,ii,j,lg,ll,m,nt,ur 882 integer t(20) 883 DOUBLE PRECISION ge,gn,gs,gw,gpe,gpn,gps,gpw,h,phi0,phi1, 884 + psi0,psi1,s,sew,sns,v0,v1,xibar 885 DOUBLE PRECISION g(0:8,256),g0(0:8),g1(0:8) 886 887 external loesswarn,ehg184 888c locate enclosing cell 889 nt=1 890 t(nt)=1 891 j=1 892c top of while loop 893 3 if(.not.(a(j).ne.0))goto 4 894 nt=nt+1 895 if(z(a(j)).le.xi(j))then 896 i1=lo(j) 897 else 898 i1=hi(j) 899 end if 900 t(nt)=i1 901 if(.not.(nt.lt.20))then 902 call loesswarn(181) 903 end if 904 j=t(nt) 905 goto 3 906c bottom of while loop 907 4 continue 908c tensor 909 do 5 i12=1,vc 910 do 6 i11=0,d 911 g(i11,i12)=vval(i11,c(i12,j)) 912 6 continue 913 5 continue 914 lg=vc 915 ll=c(1,j) 916 ur=c(vc,j) 917 do 7 i=d,1,-1 918 h=(z(i)-v(ll,i))/(v(ur,i)-v(ll,i)) 919 if(h.lt.-.001D0)then 920 call ehg184('eval ',z(1),d,1) 921 call ehg184('lowerlimit ',v(ll,1),d,nvmax) 922 else 923 if(1.001D0.lt.h)then 924 call ehg184('eval ',z(1),d,1) 925 call ehg184('upperlimit ',v(ur,1),d,nvmax) 926 end if 927 end if 928 if(-.001D0.le.h)then 929 i2=(h.le.1.001D0) 930 else 931 i2=.false. 932 end if 933 if(.not.i2)then 934 call loesswarn(122) 935 end if 936 lg=int(DBLE(lg)/2.D0) 937 do 8 ig=1,lg 938c Hermite basis 939 phi0=(1-h)**2*(1+2*h) 940 phi1=h**2*(3-2*h) 941 psi0=h*(1-h)**2 942 psi1=h**2*(h-1) 943 g(0,ig)=phi0*g(0,ig) + phi1*g(0,ig+lg) + 944 + (psi0*g(i,ig)+psi1*g(i,ig+lg)) * (v(ur,i)-v(ll,i)) 945 do 9 ii=1,i-1 946 g(ii,ig)=phi0*g(ii,ig)+phi1*g(ii,ig+lg) 947 9 continue 948 8 continue 949 7 continue 950 s=g(0,1) 951c blending 952 if(d.eq.2)then 953c ----- North ----- 954 v0=v(ll,1) 955 v1=v(ur,1) 956 do 10 i11=0,d 957 g0(i11)=vval(i11,c(3,j)) 958 10 continue 959 do 11 i11=0,d 960 g1(i11)=vval(i11,c(4,j)) 961 11 continue 962 xibar=v(ur,2) 963 m=nt-1 964c top of while loop 965 12 if(m.eq.0)then 966 i4=.true. 967 else 968 if(a(t(m)).eq.2)then 969 i3=(xi(t(m)).eq.xibar) 970 else 971 i3=.false. 972 end if 973 i4=i3 974 end if 975 if(.not.(.not.i4))goto 13 976 m=m-1 977c voidp junk 978 goto 12 979c bottom of while loop 980 13 if(m.ge.1)then 981 m=hi(t(m)) 982c top of while loop 983 14 if(.not.(a(m).ne.0))goto 15 984 if(z(a(m)).le.xi(m))then 985 m=lo(m) 986 else 987 m=hi(m) 988 end if 989 goto 14 990c bottom of while loop 991 15 if(v0.lt.v(c(1,m),1))then 992 v0=v(c(1,m),1) 993 do 16 i11=0,d 994 g0(i11)=vval(i11,c(1,m)) 995 16 continue 996 end if 997 if(v(c(2,m),1).lt.v1)then 998 v1=v(c(2,m),1) 999 do 17 i11=0,d 1000 g1(i11)=vval(i11,c(2,m)) 1001 17 continue 1002 end if 1003 end if 1004 h=(z(1)-v0)/(v1-v0) 1005c Hermite basis 1006 phi0=(1-h)**2*(1+2*h) 1007 phi1=h**2*(3-2*h) 1008 psi0=h*(1-h)**2 1009 psi1=h**2*(h-1) 1010 gn=phi0*g0(0)+phi1*g1(0)+(psi0*g0(1)+psi1*g1(1))*(v1-v0) 1011 gpn=phi0*g0(2)+phi1*g1(2) 1012c ----- South ----- 1013 v0=v(ll,1) 1014 v1=v(ur,1) 1015 do 18 i11=0,d 1016 g0(i11)=vval(i11,c(1,j)) 1017 18 continue 1018 do 19 i11=0,d 1019 g1(i11)=vval(i11,c(2,j)) 1020 19 continue 1021 xibar=v(ll,2) 1022 m=nt-1 1023c top of while loop 1024 20 if(m.eq.0)then 1025 i6=.true. 1026 else 1027 if(a(t(m)).eq.2)then 1028 i5=(xi(t(m)).eq.xibar) 1029 else 1030 i5=.false. 1031 end if 1032 i6=i5 1033 end if 1034 if(.not.(.not.i6))goto 21 1035 m=m-1 1036c voidp junk 1037 goto 20 1038c bottom of while loop 1039 21 if(m.ge.1)then 1040 m=lo(t(m)) 1041c top of while loop 1042 22 if(.not.(a(m).ne.0))goto 23 1043 if(z(a(m)).le.xi(m))then 1044 m=lo(m) 1045 else 1046 m=hi(m) 1047 end if 1048 goto 22 1049c bottom of while loop 1050 23 if(v0.lt.v(c(3,m),1))then 1051 v0=v(c(3,m),1) 1052 do 24 i11=0,d 1053 g0(i11)=vval(i11,c(3,m)) 1054 24 continue 1055 end if 1056 if(v(c(4,m),1).lt.v1)then 1057 v1=v(c(4,m),1) 1058 do 25 i11=0,d 1059 g1(i11)=vval(i11,c(4,m)) 1060 25 continue 1061 end if 1062 end if 1063 h=(z(1)-v0)/(v1-v0) 1064c Hermite basis 1065 phi0=(1-h)**2*(1+2*h) 1066 phi1=h**2*(3-2*h) 1067 psi0=h*(1-h)**2 1068 psi1=h**2*(h-1) 1069 gs=phi0*g0(0)+phi1*g1(0)+(psi0*g0(1)+psi1*g1(1))*(v1-v0) 1070 gps=phi0*g0(2)+phi1*g1(2) 1071c ----- East ----- 1072 v0=v(ll,2) 1073 v1=v(ur,2) 1074 do 26 i11=0,d 1075 g0(i11)=vval(i11,c(2,j)) 1076 26 continue 1077 do 27 i11=0,d 1078 g1(i11)=vval(i11,c(4,j)) 1079 27 continue 1080 xibar=v(ur,1) 1081 m=nt-1 1082c top of while loop 1083 28 if(m.eq.0)then 1084 i8=.true. 1085 else 1086 if(a(t(m)).eq.1)then 1087 i7=(xi(t(m)).eq.xibar) 1088 else 1089 i7=.false. 1090 end if 1091 i8=i7 1092 end if 1093 if(.not.(.not.i8))goto 29 1094 m=m-1 1095c voidp junk 1096 goto 28 1097c bottom of while loop 1098 29 if(m.ge.1)then 1099 m=hi(t(m)) 1100c top of while loop 1101 30 if(.not.(a(m).ne.0))goto 31 1102 if(z(a(m)).le.xi(m))then 1103 m=lo(m) 1104 else 1105 m=hi(m) 1106 end if 1107 goto 30 1108c bottom of while loop 1109 31 if(v0.lt.v(c(1,m),2))then 1110 v0=v(c(1,m),2) 1111 do 32 i11=0,d 1112 g0(i11)=vval(i11,c(1,m)) 1113 32 continue 1114 end if 1115 if(v(c(3,m),2).lt.v1)then 1116 v1=v(c(3,m),2) 1117 do 33 i11=0,d 1118 g1(i11)=vval(i11,c(3,m)) 1119 33 continue 1120 end if 1121 end if 1122 h=(z(2)-v0)/(v1-v0) 1123c Hermite basis 1124 phi0=(1-h)**2*(1+2*h) 1125 phi1=h**2*(3-2*h) 1126 psi0=h*(1-h)**2 1127 psi1=h**2*(h-1) 1128 ge=phi0*g0(0)+phi1*g1(0)+(psi0*g0(2)+psi1*g1(2))*(v1-v0) 1129 gpe=phi0*g0(1)+phi1*g1(1) 1130c ----- West ----- 1131 v0=v(ll,2) 1132 v1=v(ur,2) 1133 do 34 i11=0,d 1134 g0(i11)=vval(i11,c(1,j)) 1135 34 continue 1136 do 35 i11=0,d 1137 g1(i11)=vval(i11,c(3,j)) 1138 35 continue 1139 xibar=v(ll,1) 1140 m=nt-1 1141c top of while loop 1142 36 if(m.eq.0)then 1143 i10=.true. 1144 else 1145 if(a(t(m)).eq.1)then 1146 i9=(xi(t(m)).eq.xibar) 1147 else 1148 i9=.false. 1149 end if 1150 i10=i9 1151 end if 1152 if(.not.(.not.i10))goto 37 1153 m=m-1 1154c voidp junk 1155 goto 36 1156c bottom of while loop 1157 37 if(m.ge.1)then 1158 m=lo(t(m)) 1159c top of while loop 1160 38 if(.not.(a(m).ne.0))goto 39 1161 if(z(a(m)).le.xi(m))then 1162 m=lo(m) 1163 else 1164 m=hi(m) 1165 end if 1166 goto 38 1167c bottom of while loop 1168 39 if(v0.lt.v(c(2,m),2))then 1169 v0=v(c(2,m),2) 1170 do 40 i11=0,d 1171 g0(i11)=vval(i11,c(2,m)) 1172 40 continue 1173 end if 1174 if(v(c(4,m),2).lt.v1)then 1175 v1=v(c(4,m),2) 1176 do 41 i11=0,d 1177 g1(i11)=vval(i11,c(4,m)) 1178 41 continue 1179 end if 1180 end if 1181 h=(z(2)-v0)/(v1-v0) 1182c Hermite basis 1183 phi0=(1-h)**2*(1+2*h) 1184 phi1=h**2*(3-2*h) 1185 psi0=h*(1-h)**2 1186 psi1=h**2*(h-1) 1187 gw=phi0*g0(0)+phi1*g1(0)+(psi0*g0(2)+psi1*g1(2))*(v1-v0) 1188 gpw=phi0*g0(1)+phi1*g1(1) 1189c NS 1190 h=(z(2)-v(ll,2))/(v(ur,2)-v(ll,2)) 1191c Hermite basis 1192 phi0=(1-h)**2*(1+2*h) 1193 phi1=h**2*(3-2*h) 1194 psi0=h*(1-h)**2 1195 psi1=h**2*(h-1) 1196 sns=phi0*gs+phi1*gn+(psi0*gps+psi1*gpn)*(v(ur,2)-v(ll,2)) 1197c EW 1198 h=(z(1)-v(ll,1))/(v(ur,1)-v(ll,1)) 1199c Hermite basis 1200 phi0=(1-h)**2*(1+2*h) 1201 phi1=h**2*(3-2*h) 1202 psi0=h*(1-h)**2 1203 psi1=h**2*(h-1) 1204 sew=phi0*gw+phi1*ge+(psi0*gpw+psi1*gpe)*(v(ur,1)-v(ll,1)) 1205 s=(sns+sew)-s 1206 end if 1207 ehg128=s 1208 return 1209 end 1210 1211 integer function ifloor(x) 1212 DOUBLE PRECISION x 1213 ifloor=int(x) 1214 if(ifloor.gt.x) ifloor=ifloor-1 1215 end 1216 1217c DSIGN is unused, causes conflicts on some platforms 1218c DOUBLE PRECISION function DSIGN(a1,a2) 1219c DOUBLE PRECISION a1, a2 1220c DSIGN=DABS(a1) 1221c if(a2.ge.0)DSIGN=-DSIGN 1222c end 1223 1224 1225c ehg136() is the workhorse of lowesf(.) 1226c n = number of observations 1227c m = number of x values at which to evaluate 1228c f = span 1229c nf = min(n, floor(f * n)) 1230 subroutine ehg136(u,lm,m,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b, 1231 + od,o,ihat,w,rcond,sing,dd,tdeg,cdeg,s) 1232 integer identi,d,dd,i,i1,ihat,info,j,k,kernel,l,lm,m,n,nf, 1233 + od,sing,tdeg 1234 integer cdeg(8),psi(n) 1235 double precision f,i2,rcond,scale,tol 1236 double precision o(m,n),sigma(15),e(15,15),g(15,15),b(nf,k), 1237 $ dist(n),eta(nf),dgamma(15),q(8),qraux(15),rw(n),s(0:od,m), 1238 $ u(lm,d),w(nf),work(15),x(n,d),y(n) 1239 1240 external ehg127,loesswarn,dqrsl 1241 double precision ddot 1242 external ddot 1243 1244c V -> g 1245c U -> e 1246c Identity -> identi 1247c L -> o 1248c X -> b 1249 if(k .gt. nf-1) call loesswarn(104) 1250 if(k .gt. 15) call loesswarn(105) 1251 do 3 identi=1,n 1252 psi(identi)=identi 1253 3 continue 1254 do 4 l=1,m 1255 do 5 i1=1,d 1256 q(i1)=u(l,i1) 1257 5 continue 1258 call ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w, 1259 + rcond,sing,sigma,e,g,dgamma,qraux,work,tol,dd,tdeg,cdeg, 1260 + s(0,l)) 1261 if(ihat.eq.1)then 1262c $L sub {l,l} = 1263c V sub {1,:} SIGMA sup {+} U sup T 1264c (Q sup T W e sub i )$ 1265 if(.not.(m.eq.n))then 1266 call loesswarn(123) 1267 end if 1268c find $i$ such that $l = psi sub i$ 1269 i=1 1270c top of while loop 1271 6 if(.not.(l.ne.psi(i)))goto 7 1272 i=i+1 1273 if(.not.(i.lt.nf))then 1274 call loesswarn(123) 1275c next line is not in current dloess 1276 goto 7 1277 end if 1278 goto 6 1279c bottom of while loop 1280 7 do 8 i1=1,nf 1281 eta(i1)=0 1282 8 continue 1283 eta(i)=w(i) 1284c $eta = Q sup T W e sub i$ 1285 call dqrsl(b,nf,nf,k,qraux,eta,eta,eta,eta,eta,eta,1000, 1286 + info) 1287c $gamma = U sup T eta sub {1:k}$ 1288 do 9 i1=1,k 1289 dgamma(i1)=0 1290 9 continue 1291 do 10 j=1,k 1292 i2=eta(j) 1293 do 11 i1=1,k 1294 dgamma(i1)=dgamma(i1)+i2*e(j,i1) 1295 11 continue 1296 10 continue 1297c $gamma = SIGMA sup {+} gamma$ 1298 do 12 j=1,k 1299 if(tol.lt.sigma(j))then 1300 dgamma(j)=dgamma(j)/sigma(j) 1301 else 1302 dgamma(j)=0.d0 1303 end if 1304 12 continue 1305c voidp junk 1306c voidp junk 1307 o(l,1)=ddot(k,g(1,1),15,dgamma,1) 1308 else 1309 if(ihat.eq.2)then 1310c $L sub {l,:} = 1311c V sub {1,:} SIGMA sup {+} 1312c ( U sup T Q sup T ) W $ 1313 do 13 i1=1,n 1314 o(l,i1)=0 1315 13 continue 1316 do 14 j=1,k 1317 do 15 i1=1,nf 1318 eta(i1)=0 1319 15 continue 1320 do 16 i1=1,k 1321 eta(i1)=e(i1,j) 1322 16 continue 1323 call dqrsl(b,nf,nf,k,qraux,eta,eta,work,work,work,work 1324 + ,10000,info) 1325 if(tol.lt.sigma(j))then 1326 scale=1.d0/sigma(j) 1327 else 1328 scale=0.d0 1329 end if 1330 do 17 i1=1,nf 1331 eta(i1)=eta(i1)*(scale*w(i1)) 1332 17 continue 1333 do 18 i=1,nf 1334 o(l,psi(i))=o(l,psi(i))+g(1,j)*eta(i) 1335 18 continue 1336 14 continue 1337 end if 1338 end if 1339 4 continue 1340 return 1341 end 1342 1343c called from lowesb() ... compute fit ..?..?... 1344c somewhat similar to ehg136 1345 subroutine ehg139(v,nvmax,nv,n,d,nf,f,x,pi,psi,y,rw,trl,kernel,k, 1346 + dist,phi,eta,b,od,w,diagl,vval2,ncmax,vc,a,xi,lo,hi,c, 1347 + rcond,sing,dd,tdeg,cdeg,lq,lf,setlf,s) 1348 logical setlf 1349 integer identi,d,dd,i,i2,i3,i5,i6,ii,ileaf,info,j,k,kernel, 1350 + l,n,ncmax,nf,nleaf,nv,nvmax,od,sing,tdeg,vc 1351 integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),cdeg(8),hi(ncmax), 1352 + leaf(256),lo(ncmax),pi(n),psi(n) 1353 DOUBLE PRECISION f,i1,i4,i7,rcond,scale,term,tol,trl 1354 DOUBLE PRECISION lf(0:d,nvmax,nf),sigma(15),u(15,15),e(15,15), 1355 + b(nf,k),diagl(n),dist(n),eta(nf),DGAMMA(15),q(8),qraux(15), 1356 + rw(n),s(0:od,nv),v(nvmax,d),vval2(0:d,nv),w(nf),work(15), 1357 + x(n,d),xi(ncmax),y(n),z(8) 1358 DOUBLE PRECISION phi(n) 1359 1360 external ehg127,loesswarn,DQRSL,ehg137 1361 DOUBLE PRECISION ehg128 1362 external ehg128 1363 DOUBLE PRECISION DDOT 1364 external DDOT 1365 1366c V -> e 1367c Identity -> identi 1368c X -> b 1369c l2fit with trace(L) 1370 if(k .gt. nf-1) call loesswarn(104) 1371 if(k .gt. 15) call loesswarn(105) 1372 if(trl.ne.0)then 1373 do 3 i5=1,n 1374 diagl(i5)=0 1375 3 continue 1376 do 4 i6=1,nv 1377 do 5 i5=0,d 1378 vval2(i5,i6)=0 1379 5 continue 1380 4 continue 1381 end if 1382 do 6 identi=1,n 1383 psi(identi)=identi 1384 6 continue 1385 do 7 l=1,nv 1386 do 8 i5=1,d 1387 q(i5)=v(l,i5) 1388 8 continue 1389 call ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w, 1390 + rcond,sing,sigma,u,e,DGAMMA,qraux,work,tol,dd,tdeg,cdeg, 1391 + s(0,l)) 1392 if(trl.ne.0)then 1393c invert $psi$ 1394 do 9 i5=1,n 1395 phi(i5)=0 1396 9 continue 1397 do 10 i=1,nf 1398 phi(psi(i))=i 1399 10 continue 1400 do 11 i5=1,d 1401 z(i5)=v(l,i5) 1402 11 continue 1403 call ehg137(z,leaf,nleaf,d,ncmax,a,xi, 1404 + lo,hi) 1405 do 12 ileaf=1,nleaf 1406 do 13 ii=lo(leaf(ileaf)),hi(leaf(ileaf)) 1407 i=int(phi(pi(ii))) 1408 if(i.ne.0)then 1409 if(.not.(psi(i).eq.pi(ii)))then 1410 call loesswarn(194) 1411 end if 1412 do 14 i5=1,nf 1413 eta(i5)=0 1414 14 continue 1415 eta(i)=w(i) 1416c $eta = Q sup T W e sub i$ 1417 call DQRSL(b,nf,nf,k,qraux,eta,work,eta,eta,work, 1418 + work,1000,info) 1419 do 15 j=1,k 1420 if(tol.lt.sigma(j))then 1421 i4=DDOT(k,u(1,j),1,eta,1)/sigma(j) 1422 else 1423 i4=0.D0 1424 end if 1425 DGAMMA(j)=i4 1426 15 continue 1427 do 16 j=1,d+1 1428c bug fix 2006-07-15 for k=1, od>1. (thanks btyner@gmail.com) 1429 if(j.le.k)then 1430 vval2(j-1,l)=DDOT(k,e(j,1),15,DGAMMA,1) 1431 else 1432 vval2(j-1,l)=0.0d0 1433 end if 1434 16 continue 1435 do 17 i5=1,d 1436 z(i5)=x(pi(ii),i5) 1437 17 continue 1438 term=ehg128(z,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax, 1439 + vval2) 1440 diagl(pi(ii))=diagl(pi(ii))+term 1441 do 18 i5=0,d 1442 vval2(i5,l)=0 1443 18 continue 1444 end if 1445 13 continue 1446 12 continue 1447 end if 1448 if(setlf)then 1449c $Lf sub {:,l,:} = V SIGMA sup {+} U sup T Q sup T W$ 1450 if(.not.(k.ge.d+1))then 1451 call loesswarn(196) 1452 end if 1453 do 19 i5=1,nf 1454 lq(l,i5)=psi(i5) 1455 19 continue 1456 do 20 i6=1,nf 1457 do 21 i5=0,d 1458 lf(i5,l,i6)=0 1459 21 continue 1460 20 continue 1461 do 22 j=1,k 1462 do 23 i5=1,nf 1463 eta(i5)=0 1464 23 continue 1465 do 24 i5=1,k 1466 eta(i5)=u(i5,j) 1467 24 continue 1468 call DQRSL(b,nf,nf,k,qraux,eta,eta,work,work,work,work, 1469 + 10000,info) 1470 if(tol.lt.sigma(j))then 1471 scale=1.D0/sigma(j) 1472 else 1473 scale=0.D0 1474 end if 1475 do 25 i5=1,nf 1476 eta(i5)=eta(i5)*(scale*w(i5)) 1477 25 continue 1478 do 26 i=1,nf 1479 i7=eta(i) 1480 do 27 i5=0,d 1481 if(i5.lt.k)then 1482 lf(i5,l,i)=lf(i5,l,i)+e(1+i5,j)*i7 1483 else 1484 lf(i5,l,i)=0 1485 end if 1486 27 continue 1487 26 continue 1488 22 continue 1489 end if 1490 7 continue 1491 if(trl.ne.0)then 1492 if(n.le.0)then 1493 trl=0.D0 1494 else 1495 i3=n 1496 i1=diagl(i3) 1497 do 28 i2=i3-1,1,-1 1498 i1=diagl(i2)+i1 1499 28 continue 1500 trl=i1 1501 end if 1502 end if 1503 return 1504 end 1505 1506 subroutine lowesb(xx,yy,ww,diagl,infl,iv,wv) 1507 integer infl 1508 integer iv(*) 1509 DOUBLE PRECISION xx(*),yy(*),ww(*),diagl(*),wv(*) 1510c Var 1511 DOUBLE PRECISION trl 1512 logical setlf 1513 1514 integer ifloor 1515 external ifloor 1516 external ehg131,loesswarn,ehg183 1517 1518 if(.not.(iv(28).ne.173))then 1519 call loesswarn(174) 1520 end if 1521 if(iv(28).ne.172)then 1522 if(.not.(iv(28).eq.171))then 1523 call loesswarn(171) 1524 end if 1525 end if 1526 iv(28)=173 1527 if(infl.ne.0)then 1528 trl=1.D0 1529 else 1530 trl=0.D0 1531 end if 1532 setlf=(iv(27).ne.iv(25)) 1533 call ehg131(xx,yy,ww,trl,diagl,iv(20),iv(29),iv(3),iv(2),iv(5), 1534 + iv(17),iv(4),iv(6),iv(14),iv(19),wv(1),iv(iv(7)),iv(iv(8)), 1535 + iv(iv(9)),iv(iv(10)),iv(iv(22)),iv(iv(27)),wv(iv(11)), 1536 + iv(iv(23)),wv(iv(13)),wv(iv(12)),wv(iv(15)),wv(iv(16)), 1537 + wv(iv(18)),ifloor(iv(3)*wv(2)),wv(3),wv(iv(26)),wv(iv(24)), 1538 + wv(4),iv(30),iv(33),iv(32),iv(41),iv(iv(25)),wv(iv(34)), 1539 + setlf) 1540 if(iv(14).lt.iv(6)+DBLE(iv(4))/2.D0)then 1541 call ehg183('k-d tree limited by memory; nvmax=', 1542 + iv(14),1,1) 1543 else 1544 if(iv(17).lt.iv(5)+2)then 1545 call ehg183('k-d tree limited by memory. ncmax=', 1546 + iv(17),1,1) 1547 end if 1548 end if 1549 return 1550 end 1551 1552c lowesd() : Initialize iv(*) and v(1:4) 1553c ------ called only by loess_workspace() in ./loessc.c 1554 subroutine lowesd(iv, liv,lv, v, d,n,f,ideg,nf,nvmax, setlf) 1555 integer liv,lv, d,n, ideg,nf,nvmax, setlf 1556c setlf {Rboolean}: if(true) need L [nf x nvmax] matrices 1557 integer iv(liv) 1558 double precision f, v(lv) 1559 1560c had nf = min(n,ifloor(n*f)) 1561 1562 integer bound,i,i1,i2,j,ncmax,vc 1563 external loesswarn 1564 integer ifloor 1565 external ifloor 1566c 1567c unnecessary initialization of i1 to keep g77 -Wall happy 1568c 1569 i1 = 0 1570cc version -> versio 1571c if(.not.(versio.eq.106))then 1572c call loesswarn(100) 1573c end if 1574 iv(28)=171 1575 iv(2)=d 1576 iv(3)=n 1577 vc=2**d 1578 iv(4)=vc 1579 if(.not.(0.lt.f))then 1580 call loesswarn(120) 1581 end if 1582 1583 iv(19)=nf 1584 iv(20)=1 1585 if(ideg.eq.0)then 1586 i1=1 1587 else 1588 if(ideg.eq.1)then 1589 i1=d+1 1590 else 1591 if(ideg.eq.2)then 1592 i1=int(dble((d+2)*(d+1))/2.d0) 1593 end if 1594 end if 1595 end if 1596 iv(29)=i1 1597 iv(21)=1 1598 iv(14)=nvmax 1599 ncmax=nvmax 1600 iv(17)=ncmax 1601 iv(30)=0 1602 iv(32)=ideg 1603 if(.not.(ideg.ge.0))then 1604 call loesswarn(195) 1605 end if 1606 if(.not.(ideg.le.2))then 1607 call loesswarn(195) 1608 end if 1609 iv(33)=d 1610 do 3 i2=41,49 1611 iv(i2)=ideg 1612 3 continue 1613 iv(7)=50 1614 iv(8)=iv(7)+ncmax 1615 iv(9)=iv(8)+vc*ncmax 1616 iv(10)=iv(9)+ncmax 1617 iv(22)=iv(10)+ncmax 1618c initialize permutation 1619 j=iv(22)-1 1620 do 4 i=1,n 1621 iv(j+i)=i 1622 4 continue 1623 iv(23)=iv(22)+n 1624 iv(25)=iv(23)+nvmax 1625 if(setlf.ne.0)then 1626 iv(27)=iv(25)+nvmax*nf 1627 else 1628 iv(27)=iv(25) 1629 end if 1630 bound=iv(27)+n 1631 if(.not.(bound-1.le.liv))then 1632 call loesswarn(102) 1633 end if 1634 iv(11)=50 1635 iv(13)=iv(11)+nvmax*d 1636 iv(12)=iv(13)+(d+1)*nvmax 1637 iv(15)=iv(12)+ncmax 1638 iv(16)=iv(15)+n 1639 iv(18)=iv(16)+nf 1640 iv(24)=iv(18)+iv(29)*nf 1641 iv(34)=iv(24)+(d+1)*nvmax 1642 if(setlf.ne.0)then 1643 iv(26)=iv(34)+(d+1)*nvmax*nf 1644 else 1645 iv(26)=iv(34) 1646 end if 1647 bound=iv(26)+nf 1648 if(.not.(bound-1.le.lv))then 1649 call loesswarn(103) 1650 end if 1651 v(1)=f 1652 v(2)=0.05d0 1653 v(3)=0.d0 1654 v(4)=1.d0 1655 return 1656 end 1657 1658 subroutine lowese(iv,wv,m,z,s) 1659 integer m 1660 integer iv(*) 1661 double precision s(m),wv(*),z(m,1) 1662 1663 external ehg133,loesswarn 1664 1665 if(.not.(iv(28).ne.172))then 1666 call loesswarn(172) 1667 end if 1668 if(.not.(iv(28).eq.173))then 1669 call loesswarn(173) 1670 end if 1671 call ehg133(iv(2),iv(4),iv(14),iv(17),iv(iv(7)),iv(iv( 1672 +8)),iv(iv(9)),iv(iv(10)),wv(iv(11)),wv(iv(13)),wv(iv(12)),m,z,s) 1673 return 1674 end 1675 1676c "direct" (non-"interpolate") fit aka predict() : 1677 subroutine lowesf(xx,yy,ww,iv,wv,m,z,l,ihat,s) 1678 integer m,ihat 1679c m = number of x values at which to evaluate 1680 integer iv(*) 1681 double precision xx(*),yy(*),ww(*),wv(*),z(m,1),l(m,*),s(m) 1682 1683 logical i1 1684 1685 external loesswarn,ehg136 1686 if(171.le.iv(28))then 1687 i1=(iv(28).le.174) 1688 else 1689 i1=.false. 1690 end if 1691 if(.not.i1)then 1692 call loesswarn(171) 1693 end if 1694 iv(28)=172 1695 if(.not.(iv(14).ge.iv(19)))then 1696 call loesswarn(186) 1697 end if 1698 1699c do the work; in ehg136() give the argument names as they are there: 1700c ehg136(u,lm,m, n, d, nf, f, x, psi, y ,rw, 1701 call ehg136(z,m,m,iv(3),iv(2),iv(19),wv(1),xx,iv(iv(22)),yy,ww, 1702c kernel, k, dist, eta, b, od,o,ihat, 1703 + iv(20),iv(29),wv(iv(15)),wv(iv(16)),wv(iv(18)),0,l,ihat, 1704c w, rcond,sing, dd, tdeg,cdeg, s) 1705 + wv(iv(26)),wv(4),iv(30),iv(33),iv(32),iv(41),s) 1706 return 1707 end 1708 1709c Called either from loess_raw() only for case [surf_stat = "interpolate/exact"], or 1710c from loess_ise() {used only when 'se = TRUE' and surface = "interpolate"} 1711 subroutine lowesl(iv,wv,m,z,l) 1712 integer m 1713 integer iv(*) 1714 double precision l(m,*), wv(*), z(m,1) 1715 1716 external loesswarn,ehg191 1717 1718 if(.not.(iv(28).ne.172))then 1719 call loesswarn(172) 1720 end if 1721 if(.not.(iv(28).eq.173))then 1722 call loesswarn(173) 1723 end if 1724 if(.not.(iv(26).ne.iv(34)))then 1725 call loesswarn(175) 1726 end if 1727 call ehg191(m,z,l,iv(2),iv(3),iv(19),iv(6),iv(17),iv(4),iv(iv(7)), 1728 + wv(iv(12)),iv(iv(10)),iv(iv(9)),iv(iv(8)),wv(iv(11)),iv(14), 1729 + wv(iv(24)),wv(iv(34)),iv(iv(25))) 1730 return 1731 end 1732 1733c Not used 1734 subroutine lowesr(yy,iv,wv) 1735 integer iv(*) 1736 DOUBLE PRECISION yy(*),wv(*) 1737 1738 external loesswarn,ehg192 1739 if(.not.(iv(28).ne.172))then 1740 call loesswarn(172) 1741 end if 1742 if(.not.(iv(28).eq.173))then 1743 call loesswarn(173) 1744 end if 1745 call ehg192(yy,iv(2),iv(3),iv(19),iv(6),iv(14),wv(iv(13)), 1746 + wv(iv(34)),iv(iv(25))) 1747 return 1748 end 1749 1750c Update Robustness Weights -- called via .Fortran() from R's simpleLoess() in ../R/loess.R 1751c 1752 subroutine lowesw(res,n,rw,pi) 1753c Tranliterated from Devlin's ratfor 1754 1755c implicit none 1756c Args 1757 integer n 1758 double precision res(n), rw(n) 1759 integer pi(n) 1760c Var 1761 integer identi,i,i1,nh 1762 double precision cmad,rsmall 1763 1764 integer ifloor 1765 double precision d1mach 1766 1767 external ehg106 1768 external ifloor 1769 external d1mach 1770 1771c Identity -> identi 1772 1773c find median of absolute residuals 1774 do 3 i1=1,n 1775 rw(i1)=dabs(res(i1)) 1776 3 continue 1777 do 4 identi=1,n 1778 pi(identi)=identi 1779 4 continue 1780 nh=ifloor(dble(n)/2.d0)+1 1781c partial sort to find 6*mad 1782 call ehg106(1,n,nh,1,rw,pi,n) 1783 if((n-nh)+1.lt.nh)then 1784 call ehg106(1,nh-1,nh-1,1,rw,pi,n) 1785 cmad=3*(rw(pi(nh))+rw(pi(nh-1))) 1786 else 1787 cmad=6*rw(pi(nh)) 1788 end if 1789 rsmall=d1mach(1) 1790 if(cmad.lt.rsmall)then 1791 do 5 i1=1,n 1792 rw(i1)=1 1793 5 continue 1794 else 1795 do 6 i=1,n 1796 if(cmad*0.999d0.lt.rw(i))then 1797 rw(i)=0 1798 else 1799 if(cmad*0.001d0.lt.rw(i))then 1800 rw(i)=(1-(rw(i)/cmad)**2)**2 1801 else 1802 rw(i)=1 1803 end if 1804 end if 1805 6 continue 1806 end if 1807 return 1808 end 1809 1810c Compute Pseudovalues -- called via .Fortran() from R's simpleLoess() in ../R/loess.R 1811c in the case of robustness iterations 1812 subroutine lowesp(n,y,yhat,pwgts,rwgts,pi,ytilde) 1813 integer n 1814 integer pi(n) 1815 double precision y(n),yhat(n),pwgts(n),rwgts(n),ytilde(n) 1816c Var 1817 double precision c,i1,i4,mad 1818 integer i2,i3,i,m 1819 1820 external ehg106 1821 integer ifloor 1822 external ifloor 1823c median absolute deviation (using partial sort): 1824 do 3 i=1,n 1825 ytilde(i)=dabs(y(i)-yhat(i))*dsqrt(pwgts(i)) 1826 pi(i) = i 1827 3 continue 1828 m=ifloor(dble(n)/2.d0)+1 1829 call ehg106(1,n,m,1,ytilde,pi,n) 1830 if((n-m)+1.lt.m)then 1831 call ehg106(1,m-1,m-1,1,ytilde,pi,n) 1832 mad=(ytilde(pi(m-1))+ytilde(pi(m)))/2 1833 else 1834 mad=ytilde(pi(m)) 1835 end if 1836c magic constant 1837 c=(6*mad)**2/5 1838 do 5 i=1,n 1839 ytilde(i)= 1 - ((y(i)-yhat(i))**2 * pwgts(i))/c 1840 5 continue 1841 do 6 i=1,n 1842 ytilde(i)=ytilde(i)*dsqrt(rwgts(i)) 1843 6 continue 1844 if(n.le.0)then 1845 i4=0.d0 1846 else 1847 i3=n 1848 i1=ytilde(i3) 1849 do 7 i2=i3-1,1,-1 1850 i1=ytilde(i2)+i1 1851 7 continue 1852 i4=i1 1853 end if 1854 c=n/i4 1855c pseudovalues 1856 do 8 i=1,n 1857 ytilde(i)=yhat(i) + (c*rwgts(i))*(y(i)-yhat(i)) 1858 8 continue 1859 return 1860 end 1861 1862 subroutine ehg124(ll,uu,d,n,nv,nc,ncmax,vc,x,pi,a,xi,lo,hi,c,v, 1863 + vhit,nvmax,fc,fd,dd) 1864 1865 integer ll,uu,d,n,nv,nc,ncmax,vc,nvmax,fc,dd 1866 integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax),pi(n),vhit(nvmax) 1867 DOUBLE PRECISION fd, v(nvmax,d),x(n,d),xi(ncmax) 1868 1869 logical i1,i2,leaf 1870 integer i4,inorm2,k,l,m,p,u, upper, lower, check, offset 1871 DOUBLE PRECISION diam,diag(8),sigma(8) 1872 1873 external ehg125,ehg106,ehg129 1874 integer IDAMAX 1875 external IDAMAX 1876 p=1 1877 l=ll 1878 u=uu 1879 lo(p)=l 1880 hi(p)=u 1881c top of while loop 1882 3 if(.not.(p.le.nc))goto 4 1883 do 5 i4=1,dd 1884 diag(i4)=v(c(vc,p),i4)-v(c(1,p),i4) 1885 5 continue 1886 diam=0 1887 do 6 inorm2=1,dd 1888 diam=diam+diag(inorm2)**2 1889 6 continue 1890 diam=DSQRT(diam) 1891 if((u-l)+1.le.fc)then 1892 i1=.true. 1893 else 1894 i1=(diam.le.fd) 1895 end if 1896 if(i1)then 1897 leaf=.true. 1898 else 1899 if(ncmax.lt.nc+2)then 1900 i2=.true. 1901 else 1902 i2=(nvmax.lt.nv+DBLE(vc)/2.D0) 1903 end if 1904 leaf=i2 1905 end if 1906 if(.not.leaf)then 1907 call ehg129(l,u,dd,x,pi,n,sigma) 1908 k=IDAMAX(dd,sigma,1) 1909 m=int(DBLE(l+u)/2.D0) 1910 call ehg106(l,u,m,1,x(1,k),pi,n) 1911 1912c all ties go with hi son 1913c top of while loop 1914c bug fix from btyner@gmail.com 2006-07-20 1915 offset = 0 1916 7 if(((m+offset).ge.u).or.((m+offset).lt.l))goto 8 1917 if(offset .lt. 0)then 1918 lower = l 1919 check = m + offset 1920 upper = check 1921 else 1922 lower = m + offset + 1 1923 check = lower 1924 upper = u 1925 end if 1926 call ehg106(lower,upper,check,1,x(1,k),pi,n) 1927 if(x(pi(m + offset),k).eq.x(pi(m+offset+1),k))then 1928 offset = -offset 1929 if(offset .ge. 0) then 1930 offset = offset + 1 1931 end if 1932 goto 7 1933 else 1934 m = m + offset 1935 goto 8 1936 end if 1937 1938c bottom of while loop 1939 8 if(v(c(1,p),k).eq.x(pi(m),k))then 1940 leaf=.true. 1941 else 1942 leaf=(v(c(vc,p),k).eq.x(pi(m),k)) 1943 end if 1944 end if 1945 if(leaf)then 1946 a(p)=0 1947 else 1948 a(p)=k 1949 xi(p)=x(pi(m),k) 1950c left son 1951 nc=nc+1 1952 lo(p)=nc 1953 lo(nc)=l 1954 hi(nc)=m 1955c right son 1956 nc=nc+1 1957 hi(p)=nc 1958 lo(nc)=m+1 1959 hi(nc)=u 1960 call ehg125(p,nv,v,vhit,nvmax,d,k,xi(p),2**(k-1),2**(d-k), 1961 + c(1,p),c(1,lo(p)),c(1,hi(p))) 1962 end if 1963 p=p+1 1964 l=lo(p) 1965 u=hi(p) 1966 goto 3 1967c bottom of while loop 1968 4 return 1969 end 1970 1971 subroutine ehg129(l,u,d,x,pi,n,sigma) 1972 integer d,execnt,i,k,l,n,u 1973 integer pi(n) 1974 DOUBLE PRECISION machin,alpha,beta,t 1975 DOUBLE PRECISION sigma(d),x(n,d) 1976 DOUBLE PRECISION D1MACH 1977 external D1MACH 1978 save machin,execnt 1979 data execnt /0/ 1980c MachInf -> machin 1981 execnt=execnt+1 1982 if(execnt.eq.1)then 1983c initialize d1mach(2) === DBL_MAX: 1984 machin=D1MACH(2) 1985 end if 1986 do 3 k=1,d 1987 alpha=machin 1988 beta=-machin 1989 do 4 i=l,u 1990 t=x(pi(i),k) 1991 alpha=min(alpha,x(pi(i),k)) 1992 beta=max(beta,t) 1993 4 continue 1994 sigma(k)=beta-alpha 1995 3 continue 1996 return 1997 end 1998 1999c {called only from ehg127} purpose...?... 2000 subroutine ehg137(z,leaf,nleaf,d,ncmax,a,xi,lo,hi) 2001 integer d,nleaf 2002 integer leaf(256),a(ncmax),hi(ncmax),lo(ncmax),pstack(20) 2003 DOUBLE PRECISION z(d),xi(ncmax) 2004 2005 integer p,stackt 2006 2007 external loesswarn 2008c stacktop -> stackt 2009c find leaf cells affected by $z$ 2010 stackt=0 2011 p=1 2012 nleaf=0 2013c top of while loop 2014 3 if(.not.(0.lt.p))goto 4 2015 if(a(p).eq.0)then 2016c leaf 2017 nleaf=nleaf+1 2018 leaf(nleaf)=p 2019c Pop 2020 if(stackt.ge.1)then 2021 p=pstack(stackt) 2022 else 2023 p=0 2024 end if 2025 stackt=max(0,stackt-1) 2026 else 2027 if(z(a(p)).eq.xi(p))then 2028c Push 2029 stackt=stackt+1 2030 if(.not.(stackt.le.20))then 2031 call loesswarn(187) 2032 end if 2033 pstack(stackt)=hi(p) 2034 p=lo(p) 2035 else 2036 if(z(a(p)).le.xi(p))then 2037 p=lo(p) 2038 else 2039 p=hi(p) 2040 end if 2041 end if 2042 end if 2043 goto 3 2044c bottom of while loop 2045 4 if(.not.(nleaf.le.256))then 2046 call loesswarn(185) 2047 end if 2048 return 2049 end 2050 2051C-- For Error messaging, call the "a" routines at the bottom of ./loessc.c : 2052 subroutine ehg183(s, i, n, inc) 2053 character s*(*) 2054 integer i, n, inc 2055 call ehg183a(s, len(s), i, n, inc) 2056 end 2057 2058 subroutine ehg184(s, x, n, inc) 2059 character s*(*) 2060 double precision x 2061 integer n, inc 2062 call ehg184a(s, len(s), x, n, inc) 2063 end 2064