1C 2C Modified from the SMART package by J.H. Friedman, 10/10/84 3C Main change is to add spline smoothing modified from BRUTO, 4C calling code written for smooth.spline in S. 5C 6C B.D. Ripley (ripley@stats.ox.ac.uk) 1994-7. 7C 8C 9 subroutine smart(m,mu,p,q,n, w,x,y,ww,smod,nsmod, 10 & sp,nsp,dp,ndp,edf) 11 12 integer m,mu,p,q,n, nsmod, nsp,ndp 13 double precision x(p,n),y(q,n),w(n),ww(q),smod(nsmod), 14 & sp(nsp),edf(m),dp(ndp) 15 smod(1)=m 16 smod(2)=p 17 smod(3)=q 18 smod(4)=n 19 call smart1(m,mu,p,q,n, w,x,y,ww, smod(6),smod(q+6), 20 & smod(q+7),smod(q+7+p*m),smod(q+7+m*(p+q)), 21 & smod(q+7+m*(p+q+n)),smod(q+7+m*(p+q+2*n)), 22 & sp,sp(q*n+1),sp(n*(q+15)+1),sp(n*(q+15)+q+1), 23 & dp,smod(5),edf) 24 return 25 end 26 27 subroutine smart1(m,mu,p,q,n, w,x,y,ww, yb,ys, 28 & a,b,f, 29 & t,asr, 30 & r,sc,bt,g, 31 & dp,flm,edf) 32 33 integer m,mu,p,q,n 34 double precision w(n),x(p,n),y(*),ww(q), yb(q), ys 35 double precision a(p,m),b(q,m),f(n,m),t(n,m), asr(15),asr1 36 double precision r(q,n),sc(n,15),bt(q),g(p,3) 37 double precision dp(*), flm,edf(m) 38C ^^^ really (ndb) of smart(.) 39 integer i,j,l, lm 40 double precision sw,s 41c Common Vars 42 double precision span,alpha,big 43 integer ifl,lf 44 common /pprpar/ ifl,lf,span,alpha,big 45 46 double precision conv, cutmin,fdel,cjeps 47 integer maxit,mitone, mitcj 48 common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj 49 50 sw=0d0 51 do j=1,n 52 sw=sw+w(j) 53 end do 54 do j=1,n 55 do i=1,q 56 r(i,j)=y(q*(j-1)+i) 57 end do 58 end do 59 do i=1,q 60 s=0d0 61 do j=1,n 62 s=s+w(j)*r(i,j) 63 end do 64 yb(i)=s/sw 65 end do 66c yb is vector of means 67 do j=1,n 68 do i=1,q 69 r(i,j)=r(i,j)-yb(i) 70 end do 71 end do 72 ys=0.d0 73 do i=1,q 74 s=0.d0 75 do j=1,n 76 s=s+w(j)*r(i,j)**2 77 end do 78 ys=ys+ww(i)*s/sw 79 end do 80 if(ys .gt. 0d0) goto 311 81c ys is the overall standard deviation -- quit if zero 82 return 83 84 311 continue 85 ys=sqrt(ys) 86 s=1.d0/ys 87 do j=1,n 88 do i=1,q 89 r(i,j)=r(i,j)*s 90 end do 91 end do 92 93c r is now standardized residuals 94c subfit adds up to m terms one at time; lm is the number fitted. 95 call subfit(m,p,q,n,w,sw,x,r,ww,lm,a,b,f,t,asr(1),sc,bt,g,dp,edf) 96 if(lf.le.0) go to 9999 97 call fulfit(lm,lf,p,q,n,w,sw,x,r,ww,a,b,f,t,asr,sc,bt,g,dp,edf) 98C REPEAT 99 371 continue 100 do l=1,lm 101 sc(l,1)=l+0.1d0 102 s=0d0 103 do i=1,q 104 s=s+ww(i)*abs(b(i,l)) 105 end do 106 sc(l,2)=-s 107 end do 108 call sort(sc(1,2),sc,1,lm) 109 do j=1,n 110 do i=1,q 111 r(i,j)=y(q*(j-1)+i) 112 end do 113 end do 114 115 do i=1,q 116 do j=1,n 117 r(i,j)=r(i,j)-yb(i) 118 s=0.d0 119 do l=1,lm 120 s=s+b(i,l)*f(j,l) 121 end do 122 r(i,j)=r(i,j)/ys-s 123 end do 124 end do 125 126 if(lm.le.mu) goto 9999 127c back to integer: 128 l=int(sc(lm,1)) 129 asr1=0d0 130 do j=1,n 131 do i=1,q 132 r(i,j)=r(i,j)+b(i,l)*f(j,l) 133 asr1=asr1+w(j)*ww(i)*r(i,j)**2 134 end do 135 end do 136 137 asr1=asr1/sw 138 asr(1)=asr1 139 if(l .ge. lm) goto 591 140 do i=1,p 141 a(i,l)=a(i,lm) 142 end do 143 do i=1,q 144 b(i,l)=b(i,lm) 145 end do 146 do j=1,n 147 f(j,l)=f(j,lm) 148 t(j,l)=t(j,lm) 149 end do 150 151 591 continue 152 lm=lm-1 153 call fulfit(lm,lf,p,q,n,w,sw,x,r,ww,a,b,f,t,asr,sc,bt,g,dp,edf) 154 goto 371 155C END REPEAT 156 9999 continue 157 flm=lm 158 return 159 end 160 161 subroutine subfit(m,p,q,n,w,sw,x,r,ww,lm,a,b,f,t,asr,sc, 162 & bt,g,dp,edf) 163c Args 164 integer m,p,q,n, lm 165 double precision w(n),sw, x(p,n),r(q,n),ww(q),a(p,m),b(q,m), 166 & f(n,m), t(n,m), asr(15), sc(n,15), bt(q), g(p,3), edf(m) 167 double precision dp(*) 168c Var 169 integer i,j,l, iflsv 170 double precision asrold 171c Common Vars 172 double precision span,alpha,big 173 integer ifl,lf 174 common /pprpar/ ifl,lf,span,alpha,big 175 176 double precision conv, cutmin,fdel,cjeps 177 integer maxit,mitone, mitcj 178 common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj 179 180 asr(1)=big 181 lm=0 182 do 100 l=1,m 183 call rchkusr() 184 lm=lm+1 185 asrold=asr(1) 186 call newb(lm,q,ww,b) 187c does 'edf' mean 'edf(1)' or 'edf(l)'? 188 call onetrm(0,p,q,n,w,sw,x,r,ww,a(1,lm),b(1,lm), 189 & f(1,lm),t(1,lm),asr(1),sc,g,dp,edf(1)) 190 do 20 j=1,n 191 do 10 i=1,q 192 r(i,j)=r(i,j)-b(i,lm)*f(j,lm) 193 10 continue 194 20 continue 195 if(lm.eq.1) goto 100 196 if(lf.gt.0) then 197 if(lm.eq.m) return 198 iflsv=ifl 199 ifl=0 200 call fulfit(lm,1,p,q,n,w,sw,x,r,ww,a,b,f,t,asr,sc,bt, 201 & g,dp, edf) 202 ifl=iflsv 203 endif 204 if(asr(1).le.0d0.or.(asrold-asr(1))/asrold.lt.conv) return 205100 continue 206 return 207 end 208 209 subroutine fulfit(lm,lbf,p,q,n,w,sw,x,r,ww,a,b,f,t, 210 & asr,sc,bt,g,dp,edf) 211c Args 212 integer lm,lbf,p,q,n 213 double precision w(n),sw,x(p,n),r(q,n),ww(q),a(p,lm),b(q,lm), 214 & f(n,lm),t(n,lm),asr(1+lm), sc(n,15),bt(q),g(p,3), edf(lm) 215 double precision dp(*) 216c Var 217 double precision asri, fsv, asrold 218 integer i,j,iter,lp,isv 219c Common Vars 220 double precision span,alpha,big 221 integer ifl,lf 222 common /pprpar/ ifl,lf,span,alpha,big 223 224 double precision conv, cutmin,fdel,cjeps 225 integer maxit,mitone, mitcj 226 common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj 227 228 if(lbf.le.0) return 229 asri=asr(1) 230 fsv=cutmin 231 isv=mitone 232 if(lbf .lt. 3) then 233 cutmin=1d0 234 mitone=lbf-1 235 endif 236 iter=0 237C Outer loop: 2381000 continue 239 asrold=asri 240 iter=iter+1 241 do 100 lp=1,lm 242 do 10 i=1,q 243 bt(i)=b(i,lp) 244 10 continue 245 do 20 i=1,p 246 g(i,3)=a(i,lp) 247 20 continue 248 do 35 j=1,n 249 do 30 i=1,q 250 r(i,j)=r(i,j)+bt(i)*f(j,lp) 251 30 continue 252 35 continue 253 254 call onetrm(1,p,q,n,w,sw,x,r,ww,g(1,3),bt,sc(1,14),sc(1,15), 255 & asri,sc,g,dp,edf(lp)) 256 if(asri .lt. asrold) then 257 do 40 i=1,q 258 b(i,lp)=bt(i) 259 40 continue 260 do 50 i=1,p 261 a(i,lp)=g(i,3) 262 50 continue 263 do 60 j=1,n 264 f(j,lp)=sc(j,14) 265 t(j,lp)=sc(j,15) 266 60 continue 267 else 268 asri=asrold 269 endif 270 do 85 j=1,n 271 do 80 i=1,q 272 r(i,j)=r(i,j)-b(i,lp)*f(j,lp) 273 80 continue 274 85 continue 275100 continue 276 if((iter .le. maxit) .and. ((asri .gt. 0d0) 277 & .and. ((asrold-asri)/asrold .ge. conv))) goto 1000 278 cutmin=fsv 279 mitone=isv 280 if(ifl .gt. 0) then 281 asr(1+lm) = asri 282 asr(1) = asri 283 endif 284 return 285 end 286 287 subroutine onetrm(jfl,p,q,n,w,sw,x,y,ww,a,b,f,t,asr, 288 & sc,g,dp,edf) 289c Args 290 integer jfl,p,q,n 291 double precision w(n),sw, x(p,n),y(q,n),ww(q),a(p),b(q),f(n),t(n), 292 & asr, sc(n,13),g(p,2), edf 293 double precision dp(*) 294c Var 295 double precision asrold,s 296 integer i,j,iter 297c Common Vars 298 double precision span,alpha,big 299 integer ifl,lf 300 common /pprpar/ ifl,lf,span,alpha,big 301 302 double precision conv, cutmin,fdel,cjeps 303 integer maxit,mitone, mitcj 304 common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj 305 306 iter=0 307 asr=big 308C REPEAT 3091000 continue 310 iter=iter+1 311 asrold=asr 312 do 11 j=1,n 313 s=0d0 314 do 21 i=1,q 315 s=s+ww(i)*b(i)*y(i,j) 316 21 continue 317 sc(j,13)=s 318 11 continue 319 call oneone(max0(jfl,iter-1),p,n,w,sw,sc(1,13),x,a,f,t, 320 & asr,sc,g,dp,edf) 321 do 31 i=1,q 322 s=0d0 323 do 41 j=1,n 324 s=s+w(j)*y(i,j)*f(j) 325 41 continue 326 b(i)=s/sw 327 31 continue 328 asr=0d0 329 do 51 i=1,q 330 s=0d0 331 do 61 j=1,n 332 s=s+w(j)*(y(i,j)-b(i)*f(j))**2 333 61 continue 334 asr=asr+ww(i)*s/sw 335 51 continue 336 if((q .ne. 1) .and. (iter .le. maxit) .and. (asr .gt. 0d0) 337 & .and. (asrold-asr)/asrold .ge. conv) goto 1000 338 return 339 end 340 341 subroutine oneone(ist,p,n, w,sw,y,x,a,f,t,asr,sc,g,dp,edf) 342c Args 343 integer ist,p,n 344 double precision w(n),sw,y(n),x(p,n),a(p),f(n),t(n),asr, 345 & sc(n,12), g(p,2), edf, dp(*) 346c Var 347 integer i,j,k,iter 348 double precision sml, s,v,cut,asrold 349c Common Vars 350 double precision span,alpha,big 351 integer ifl,lf 352 common /pprpar/ ifl,lf,span,alpha,big 353 354 double precision conv, cutmin,fdel,cjeps 355 integer maxit,mitone, mitcj 356 common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj 357 358 sml=1d0/big 359 if(ist .le. 0) then 360 if(p .le. 1) a(1)=1d0 361 do 10 j=1,n 362 sc(j,2)=1d0 363 10 continue 364 call pprdir(p,n,w,sw,y,x,sc(1,2),a,dp) 365 endif 366 s=0d0 367 do 20 i=1,p 368 g(i,1)=0d0 369 s=s+a(i)**2 370 20 continue 371 s=1d0/sqrt(s) 372 do 30 i=1,p 373 a(i)=a(i)*s 374 30 continue 375 iter=0 376 asr=big 377 cut=1d0 378C REPEAT ----------------------------- 379 100 continue 380 iter=iter+1 381 asrold=asr 382C REPEAT [inner loop] ----- 383 60 continue 384 s=0d0 385 do 70 i=1,p 386 g(i,2)=a(i)+g(i,1) 387 s=s+g(i,2)**2 388 70 continue 389 s=1.d0/sqrt(s) 390 do 80 i=1,p 391 g(i,2)=g(i,2)*s 392 80 continue 393 do 90 j=1,n 394 sc(j,1)=j+0.1d0 395 s=0.d0 396 do 91 i=1,p 397 s=s+g(i,2)*x(i,j) 398 91 continue 399 sc(j,11)=s 400 90 continue 401 call sort(sc(1,11),sc,1,n) 402 do 110 j=1,n 403 k=int(sc(j,1)) 404 sc(j,2)=y(k) 405 sc(j,3)=max(w(k),sml) 406 110 continue 407 call supsmu(n,sc(1,11),sc(1,2),sc(1,3),1,span,alpha, 408 & sc(1,12),sc(1,4), edf) 409 s=0d0 410 do 120 j=1,n 411 s=s+sc(j,3)*(sc(j,2)-sc(j,12))**2 412 120 continue 413 s=s/sw 414 if(s .lt. asr) goto 140 415 cut=cut*0.5d0 416 if(cut.lt.cutmin) goto 199 417 do 150 i=1,p 418 g(i,1)=g(i,1)*cut 419 150 continue 420 go to 60 421C -------- 422 140 continue 423 asr=s 424 cut=1d0 425 do 160 i=1,p 426 a(i)=g(i,2) 427 160 continue 428 do 170 j=1,n 429 k=int(sc(j,1)) 430 t(k)=sc(j,11) 431 f(k)=sc(j,12) 432 170 continue 433 if(asr.le.0d0.or.(asrold-asr)/asrold.lt.conv) goto 199 434 if(iter.gt.mitone.or.p.le.1) goto 199 435 call pprder(n,sc(1,11),sc(1,12),sc(1,3),fdel,sc(1,4),sc(1,5)) 436 do 180 j=1,n 437 k=int(sc(j,1)) 438 sc(j,5)=y(j)-f(j) 439 sc(k,6)=sc(j,4) 440 180 continue 441 call pprdir(p,n,w,sw,sc(1,5),x,sc(1,6),g,dp) 442 443 goto 100 444c-------------- 445 199 continue 446c-------------- 447 s=0d0 448 v=s 449 do 210 j=1,n 450 s=s+w(j)*f(j) 451 210 continue 452 s=s/sw 453 do 220 j=1,n 454 f(j)=f(j)-s 455 v=v+w(j)*f(j)**2 456 220 continue 457 if(v .gt. 0d0) then 458 v=1d0/sqrt(v/sw) 459 do 230 j=1,n 460 f(j)=f(j)*v 461 230 continue 462 endif 463 return 464 end 465 466 467 subroutine pprdir(p,n,w,sw,r,x,d,e,g) 468 469 integer p,n 470 double precision w(n),sw,r(n),x(p,n),d(n),e(p), g(*) 471 472 double precision s 473 integer i,j,k,l,m1,m2 474 475 double precision conv, cutmin,fdel,cjeps 476 integer maxit,mitone, mitcj 477 common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj 478 479 do 10 i=1,p 480 s=0d0 481 do 15 j=1,n 482 s=s+w(j)*d(j)*x(i,j) 483 15 continue 484 e(i)=s/sw 485 10 continue 486 k=0 487 m1=p*(p+1)/2 488 m2=m1+p 489 do 20 j=1,p 490 s=0d0 491 do 22 l=1,n 492 s=s+w(l)*r(l)*(d(l)*x(j,l)-e(j)) 493 22 continue 494 g(m1+j)=s/sw 495 do 25 i=1,j 496 s=0d0 497 do 27 l=1,n 498 s=s+w(l)*(d(l)*x(i,l)-e(i))*(d(l)*x(j,l)-e(j)) 499 27 continue 500 k=k+1 501 g(k)=s/sw 502 25 continue 503 20 continue 504 call ppconj(p,g,g(m1+1),g(m2+1),cjeps,mitcj,g(m2+p+1)) 505 do 30 i=1,p 506 e(i)=g(m2+i) 507 30 continue 508 return 509 end 510 511 subroutine ppconj(p,g,c,x,eps,maxit,sc) 512 integer p,maxit 513 double precision g(*),c(p),x(p),eps,sc(p,4) 514 515 integer i,j,im1,iter,nit 516 double precision beta,h,s,alpha,t 517 518 do 1 i=1,p 519 x(i)=0d0 520 sc(i,2)=0d0 521 1 continue 522 nit=0 523C REPEAT 52411321 continue 525 nit=nit+1 526 h=0d0 527 beta=0d0 528 do 11331 i=1,p 529 sc(i,4)=x(i) 530 s=g(i*(i-1)/2+i)*x(i) 531 im1=i-1 532 j=1 533 goto 11343 53411341 j=j+1 53511343 if(j.gt.im1) goto 11342 536 s=s+g(i*(i-1)/2+j)*x(j) 537 goto 11341 53811342 continue 539 j=i+1 540 goto 11353 54111351 j=j+1 54211353 if(j.gt.p) goto 11352 543 s=s+g(j*(j-1)/2+i)*x(j) 544 goto 11351 54511352 continue 546 sc(i,1)=s-c(i) 547 h=h+sc(i,1)**2 54811331 continue 549 if(h.le.0d0) goto 11322 550 do 11361 iter=1,p 551 do 11371 i=1,p 552 sc(i,2)=beta*sc(i,2)-sc(i,1) 55311371 continue 554 t=0d0 555 do 11381 i=1,p 556 s=g(i*(i-1)/2+i)*sc(i,2) 557 im1=i-1 558 j=1 559 goto 11393 56011391 j=j+1 56111393 if(j.gt.im1) goto 11392 562 s=s+g(i*(i-1)/2+j)*sc(j,2) 563 goto 11391 56411392 continue 565 j=i+1 566 goto 11403 56711401 j=j+1 56811403 if(j.gt.p) goto 11402 569 s=s+g(j*(j-1)/2+i)*sc(j,2) 570 goto 11401 57111402 continue 572 sc(i,3)=s 573 t=t+s*sc(i,2) 57411381 continue 575 alpha=h/t 576 s=0d0 577 do 11411 i=1,p 578 x(i)=x(i)+alpha*sc(i,2) 579 sc(i,1)=sc(i,1)+alpha*sc(i,3) 580 s=s+sc(i,1)**2 58111411 continue 582 if(s.le.0d0) goto 11362 583 beta=s/h 584 h=s 58511361 continue 58611362 continue 587 s=0d0 588 do 11421 i=1,p 589 s=dmax1(s,dabs(x(i)-sc(i,4))) 59011421 continue 591 if((s .ge. eps) .and. (nit .lt. maxit)) goto 11321 59211322 continue 593 return 594 end 595 596 subroutine pprder (n,x,s,w,fdel,d,sc) 597 integer n 598 double precision x(n),s(n),w(n), fdel, d(n),sc(n,3) 599 600 integer i,j,bl,el,bc,ec,br,er 601 double precision scale, del 602c 603c unnecessary initialization of bl el ec to keep g77 -Wall happy 604c 605 bl = 0 606 el = 0 607 ec = 0 608c 609 if(x(n) .gt. x(1)) goto 11441 610 do 11451 j=1,n 611 d(j)=0d0 61211451 continue 613 return 61411441 continue 615 i=n/4 616 j=3*i 617 scale=x(j)-x(i) 61811461 if(scale.gt.0d0) goto 11462 619 if(j.lt.n) j=j+1 620 if(i.gt.1) i=i-1 621 scale=x(j)-x(i) 622 goto 11461 62311462 continue 624 del=fdel*scale*2d0 625 do 11471 j=1,n 626 sc(j,1)=x(j) 627 sc(j,2)=s(j) 628 sc(j,3)=w(j) 62911471 continue 630 call pool (n,sc,sc(1,2),sc(1,3),del) 631 bc=0 632 br=bc 633 er=br 63411481 continue 635 br=er+1 636 er=br 63711491 if(er .ge. n) goto 11492 638 if(sc(br,1) .ne. sc(er+1,1)) goto 11511 639 er=er+1 640 goto 11521 64111511 continue 642 goto 11492 64311521 continue 644 goto 11491 64511492 continue 646 if(br .ne. 1) goto 11541 647 bl=br 648 el=er 649 goto 11481 65011541 continue 651 if(bc .ne. 0) goto 11561 652 bc=br 653 ec=er 654 do 11571 j=bl,el 655 d(j)=(sc(bc,2)-sc(bl,2))/(sc(bc,1)-sc(bl,1)) 65611571 continue 657 goto 11481 65811561 continue 659c sanity check needed for PR#13517 660 if(br.gt.n) call rexit('br is too large') 661 do 11581 j=bc,ec 662 d(j)=(sc(br,2)-sc(bl,2))/(sc(br,1)-sc(bl,1)) 66311581 continue 664 if(er .ne. n) goto 11601 665 do 11611 j=br,er 666 d(j)=(sc(br,2)-sc(bc,2))/(sc(br,1)-sc(bc,1)) 66711611 continue 668 goto 11482 66911601 continue 670 bl=bc 671 el=ec 672 bc=br 673 ec=er 674 goto 11481 67511482 continue 676 return 677 end 678 679 subroutine pool (n,x,y,w,del) 680 integer n 681 double precision x(n),y(n),w(n),del 682 683 integer i,bb,eb,br,er,bl,el 684 double precision px, py, pw 685 686 bb=0 687 eb=bb 68811621 if(eb.ge.n) goto 11622 689 bb=eb+1 690 eb=bb 69111631 if(eb .ge. n) goto 11632 692 if(x(bb) .ne. x(eb+1)) goto 11651 693 eb=eb+1 694 goto 11661 69511651 continue 696 goto 11632 69711661 continue 698 goto 11631 69911632 continue 700 if(eb .ge. n) goto 11681 701 if(x(eb+1)-x(eb) .ge. del) goto 11701 702 br=eb+1 703 er=br 70411711 if(er .ge. n) goto 11712 705 if(x(er+1) .ne. x(br)) goto 11731 706 er=er+1 707 goto 11741 70811731 continue 709 goto 11712 71011741 continue 711 goto 11711 71211712 continue 713C avoid bounds error: this was .and. but order is not guaranteed 714 if(er.lt.n) then 715 if(x(er+1)-x(er).lt.x(eb+1)-x(eb)) goto 11621 716 endif 717 eb=er 718 pw=w(bb)+w(eb) 719 px=(x(bb)*w(bb)+x(eb)*w(eb))/pw 720 py=(y(bb)*w(bb)+y(eb)*w(eb))/pw 721 do 11751 i=bb,eb 722 x(i)=px 723 y(i)=py 724 w(i)=pw 72511751 continue 72611701 continue 72711681 continue 72811761 continue 729 if(bb.le.1) goto 11762 730 if(x(bb)-x(bb-1).ge.del) goto 11762 731 bl=bb-1 732 el=bl 73311771 if(bl .le. 1) goto 11772 734 if(x(bl-1) .ne. x(el)) goto 11791 735 bl=bl-1 736 goto 11801 73711791 continue 738 goto 11772 73911801 continue 740 goto 11771 74111772 continue 742 bb=bl 743 pw=w(bb)+w(eb) 744 px=(x(bb)*w(bb)+x(eb)*w(eb))/pw 745 py=(y(bb)*w(bb)+y(eb)*w(eb))/pw 746 do 11811 i=bb,eb 747 x(i)=px 748 y(i)=py 749 w(i)=pw 75011811 continue 751 goto 11761 75211762 continue 753 goto 11621 75411622 continue 755 return 756 end 757 758 subroutine newb(lm,q,ww,b) 759 integer lm, q 760 double precision ww(q),b(q,lm) 761 762 integer i,lm1,l,l1 763 double precision s,t,sml 764c Common 765 double precision span,alpha,big 766 integer ifl,lf 767 common /pprpar/ ifl,lf,span,alpha,big 768 769 770 sml=1d0/big 771 if(q .ne. 1) goto 11831 772 b(1,lm)=1d0 773 return 77411831 continue 775 if(lm .ne. 1) goto 11851 776 do 11861 i=1,q 777 b(i,lm)=i 77811861 continue 779 return 78011851 continue 781 lm1=lm-1 782 do 11871 i=1,q 783 b(i,lm)=0d0 78411871 continue 785 t=0d0 786 do 11881 i=1,q 787 s=0d0 788 do 11891 l=1,lm1 789 s=s+abs(b(i,l)) 79011891 continue 791 b(i,lm)=s 792 t=t+s 79311881 continue 794 do 11901 i=1,q 795 b(i,lm)=ww(i)*(t-b(i,lm)) 79611901 continue 797 l1=1 798 if(lm.gt.q) l1=lm-q+1 799 do 11911 l=l1,lm1 800 s=0d0 801 t=s 802 do 11921 i=1,q 803 s=s+ww(i)*b(i,lm)*b(i,l) 804 t=t+ww(i)*b(i,l)**2 80511921 continue 806 s=s/sqrt(t) 807 do 11931 i=1,q 808 b(i,lm)=b(i,lm)-s*b(i,l) 80911931 continue 81011911 continue 811 do 11941 i=2,q 812 if(abs(b(i-1,lm)-b(i,lm)).gt.sml) return 81311941 continue 814 do 11951 i=1,q 815 b(i,lm)=i 81611951 continue 817 return 818 end 819 820 block data bkppr 821 822c Common Vars 823 double precision span,alpha,big 824 integer ifl,lf 825 common /pprpar/ ifl,lf,span,alpha,big 826 827 double precision conv, cutmin,fdel,cjeps 828 integer maxit,mitone, mitcj 829 common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj 830 831 double precision df, gcvpen 832 integer ismethod 833 logical trace 834 common /spsmooth/ df, gcvpen, ismethod, trace 835 836 data df, gcvpen, ismethod, trace /4d0, 1d0, 0, .false./ 837 838 data ifl,maxit, conv, mitone, cutmin, fdel, 839 & span,alpha, big, cjeps, mitcj, lf 840 & /6, 20, .005, 20, 0.1, 0.02, 841 & 0.0, 0.0,1.0e20,0.001, 1, 2/ 842 end 843 844 subroutine setppr(span1, alpha1, optlevel, ism, df1, gcvpen1) 845c Put 'parameters' into Common blocks 846 integer optlevel,ism 847 double precision span1,alpha1, df1, gcvpen1 848 849 double precision span,alpha,big 850 integer ifl,lf 851 common /pprpar/ ifl,lf,span,alpha,big 852 853 double precision df, gcvpen 854 integer ismethod 855 logical trace 856 common /spsmooth/ df, gcvpen, ismethod, trace 857 858 span = span1 859 lf = optlevel 860 alpha = alpha1 861 if(ism .ge. 0) then 862 ismethod = ism 863 trace = .false. 864 else 865 ismethod = -(ism+1) 866 trace = .true. 867 end if 868 df = df1 869 gcvpen = gcvpen1 870 return 871 end 872 873 subroutine fsort(mu,n,f,t,sp) 874c 875 integer mu, n 876 double precision f(n,mu),t(n,mu),sp(n,2) 877c 878 integer l,j,k 879 880 do 100 l=1,mu 881 do 10 j=1,n 882 sp(j,1)=j+0.1d0 883 sp(j,2)=f(j,l) 884 10 continue 885 call sort(t(1,l),sp,1,n) 886 do 20 j=1,n 887 k=int(sp(j,1)) 888 f(j,l)=sp(k,2) 889 20 continue 890 100 continue 891 return 892 end 893 894 subroutine pppred(np,x,smod,y,sc) 895 896 integer np 897 double precision x(np,*),y(np,*),smod(*), sc(*) 898 899 integer p,q, place,low,high, i,j,l,m,n, 900 + inp,ja,jb,jf,jt,jfl,jfh,jtl,jth, mu 901 double precision ys, s, t 902 903 m= int(smod(1)+0.1d0) 904 p= int(smod(2)+0.1d0) 905 q= int(smod(3)+0.1d0) 906 n= int(smod(4)+0.1d0) 907 mu=int(smod(5)+0.1d0) 908 ys=smod(q+6) 909 ja=q+6 910 jb=ja+p*m 911 jf=jb+m*q 912 jt=jf+n*m 913 call fsort(mu,n,smod(jf+1),smod(jt+1),sc) 914 do 100 inp = 1, np 915 ja=q+6 916 jb=ja+p*m 917 jf=jb+m*q 918 jt=jf+n*m 919 do 81 i=1,q 920 y(inp,i)=0d0 921 81 continue 922 do 91 l=1,mu 923 s=0d0 924 do 12201 j=1,p 925 s=s+smod(ja+j)*x(inp,j) 92612201 continue 927 if(s .gt. smod(jt+1)) goto 12221 928 place=1 929 go to 12230 93012221 continue 931 if(s .lt. smod(jt+n)) goto 12251 932 place=n 933 go to 12230 934 93512251 continue 936 low=0 937 high=n+1 938C WHILE 93912261 if(low+1.ge.high) goto 12262 940 place=(low+high)/2 941 t=smod(jt+place) 942 if(s.eq.t) goto 12230 943 if(s .lt. t) then 944 high=place 945 else 946 low=place 947 endif 948 goto 12261 949C END 95012262 continue 951 jfl=jf+low 952 jfh=jf+high 953 jtl=jt+low 954 jth=jt+high 955 t=smod(jfl)+(smod(jfh)-smod(jfl))*(s-smod(jtl)) / 956 & (smod(jth)-smod(jtl)) 957 go to 12300 95812230 continue 959 t=smod(jf+place) 96012300 continue 961 do 12311 i=1,q 962 y(inp,i)=y(inp,i)+smod(jb+i)*t 96312311 continue 964 ja=ja+p 965 jb=jb+q 966 jf=jf+n 967 jt=jt+n 968 91 continue 969 do 12321 i=1,q 970 y(inp,i)=ys*y(inp,i)+smod(i+5) 97112321 continue 972 100 continue 973 return 974 end 975 976c Called from R's supsmu() 977 subroutine setsmu (tr) 978 integer tr 979 980 double precision df, gcvpen 981 integer ismethod 982 logical trace 983 common /spsmooth/ df, gcvpen, ismethod, trace 984 985 ismethod = 0 986 trace = tr .ne. 0 987 return 988 end 989 990 subroutine supsmu (n,x,y,w,iper,span,alpha,smo,sc,edf) 991c 992c------------------------------------------------------------------ 993c 994c super smoother (Friedman, 1984). 995c 996c version 10/10/84 997c 998c coded and copywrite (c) 1984 by: 999c 1000c Jerome H. Friedman 1001c department of statistics 1002c and 1003c stanford linear accelerator center 1004c stanford university 1005c 1006c all rights reserved. 1007c 1008c 1009c input: 1010c n : number of observations (x,y - pairs). 1011c x(n) : ordered abscissa values. 1012c y(n) : corresponding ordinate (response) values. 1013c w(n) : weight for each (x,y) observation. 1014c iper : periodic variable flag. 1015c iper=1 => x is ordered interval variable. 1016c iper=2 => x is a periodic variable with values 1017c in the range (0.0,1.0) and period 1.0. 1018c span : smoother span (fraction of observations in window). 1019c span=0.0 <=> "cv" : automatic (variable) span selection. 1020c alpha : controls high frequency (small span) penality 1021c used with automatic span selection (bass tone control). 1022c (alpha.le.0.0 or alpha.gt.10.0 => no effect.) 1023c output: 1024c smo(n) : smoothed ordinate (response) values. 1025c scratch: 1026c sc(n,7) : internal working storage. 1027c 1028c note: 1029c for small samples (n < 40) or if there are substantial serial 1030c correlations between observations close in x - value, then 1031c a prespecified fixed span smoother (span > 0) should be 1032c used. reasonable span values are 0.2 to 0.4. 1033c 1034c------------------------------------------------------------------ 1035 1036c Args 1037 integer n, iper 1038 double precision x(n),y(n),w(n), smo(n),sc(n,7) 1039 double precision span, alpha, edf 1040c Var 1041 double precision sy,sw, a,h(n),f, scale,vsmlsq,resmin 1042 integer i,j, jper 1043 1044 double precision spans(3), big,sml,eps 1045 common /spans/ spans /consts/ big,sml,eps 1046 1047 double precision df, gcvpen 1048 integer ismethod 1049 logical trace 1050 common /spsmooth/ df, gcvpen, ismethod, trace 1051c Called from R's supsmu(), ismethod = 0, always (but not when called from ppr) 1052 1053 if (x(n).gt.x(1)) go to 30 1054c x(n) <= x(1) : boundary case: smo[.] := weighted mean( y ) 1055 sy=0d0 1056 sw=sy 1057 do 10 j=1,n 1058 sy=sy+w(j)*y(j) 1059 sw=sw+w(j) 1060 10 continue 1061 a=0d0 1062 if (sw.gt.0d0) a=sy/sw 1063 do 20 j=1,n 1064 smo(j)=a 1065 20 continue 1066 return 1067 1068C Normal Case 1069 30 continue 1070 if (ismethod .ne. 0) then ! possible only when called from ppr() 1071 call spline(n, x, y, w, smo, edf, sc) 1072 else 1073 i=n/4 1074 j=3*i 1075 scale=x(j)-x(i) ! = IQR(x) 1076 40 if (scale.gt.0d0) go to 50 1077 if (j.lt.n) j=j+1 1078 if (i.gt.1) i=i-1 1079 scale=x(j)-x(i) 1080 go to 40 1081 50 vsmlsq=(eps*scale)**2 1082 jper=iper 1083 if (iper.eq.2.and.(x(1).lt.0d0.or.x(n).gt.1d0)) jper=1 1084 if (jper.lt.1.or.jper.gt.2) jper=1 1085 if (span .gt. 0d0) then 1086 call smooth (n,x,y,w,span,jper,vsmlsq,smo,sc) 1087 return 1088 end if 1089C else "cv" (crossvalidation) from three spans[] 1090 do 70 i=1,3 1091 call smooth (n,x,y,w,spans(i),jper,vsmlsq, 1092 & sc(1,2*i-1),sc(1,7)) 1093 call smooth (n,x,sc(1,7),w,spans(2),-jper,vsmlsq, 1094 & sc(1,2*i),h) 1095 70 continue 1096 do 90 j=1,n 1097 resmin=big 1098 do 80 i=1,3 1099 if (sc(j,2*i).ge.resmin) go to 80 1100 resmin=sc(j,2*i) 1101 sc(j,7)=spans(i) 1102 80 continue 1103 if (alpha.gt.0d0 .and. alpha.le.10d0 .and. 1104 & resmin.lt.sc(j,6) .and. resmin.gt.0d0) 1105 & sc(j,7)= sc(j,7)+(spans(3)-sc(j,7)) * 1106 & max(sml,resmin/sc(j,6))**(10d0-alpha) 1107 90 continue 1108 1109 call smooth (n,x,sc(1,7),w,spans(2),-jper,vsmlsq,sc(1,2),h) 1110 do 110 j=1,n 1111 if (sc(j,2).le.spans(1)) sc(j,2)=spans(1) 1112 if (sc(j,2).ge.spans(3)) sc(j,2)=spans(3) 1113 f=sc(j,2)-spans(2) 1114 if (f.ge.0d0) go to 100 1115 f=-f/(spans(2)-spans(1)) 1116 sc(j,4)=(1d0-f)*sc(j,3)+f*sc(j,1) 1117 go to 110 1118 100 f=f/(spans(3)-spans(2)) 1119 sc(j,4)=(1d0-f)*sc(j,3)+f*sc(j,5) 1120 110 continue 1121 call smooth (n,x,sc(1,4),w,spans(1),-jper,vsmlsq,smo,h) 1122 edf = 0 1123 endif 1124 return 1125 end 1126 1127 subroutine smooth (n,x,y,w,span,iper,vsmlsq,smo,acvr) 1128c Args 1129 integer n, iper 1130 double precision x(n),y(n),w(n), span,vsmlsq, smo(n),acvr(n) 1131c Var 1132 integer i,j, in,out, jper,ibw,it, j0 1133 double precision xm,ym,var,cvar, fbw,fbo,xti,xto,tmp, a,h,sy,wt 1134 1135c will use 'trace': 1136 double precision df, gcvpen 1137 integer ismethod 1138 logical trace 1139 common /spsmooth/ df, gcvpen, ismethod, trace 1140 1141 xm=0d0 1142 ym=xm 1143 var=ym 1144 cvar=var 1145 fbw=cvar 1146 jper=iabs(iper) 1147 ibw=int(0.5d0*span*n+0.5d0) 1148 if (ibw.lt.2) ibw=2 1149 it=2*ibw+1 1150 if (it .gt. n) it = n 1151 do i=1,it 1152 j=i 1153 if (jper.eq.2) j=i-ibw-1 1154 if (j.ge.1) then 1155 xti=x(j) 1156 else ! if (j.lt.1) then 1157 j=n+j 1158 xti=x(j)-1d0 1159 end if 1160 wt=w(j) 1161 fbo=fbw 1162 fbw=fbw+wt 1163 if (fbw.gt.0d0) xm=(fbo*xm+wt*xti)/fbw 1164 if (fbw.gt.0d0) ym=(fbo*ym+wt*y(j))/fbw 1165 tmp=0d0 1166 if (fbo.gt.0d0) tmp=fbw*wt*(xti-xm)/fbo 1167 var =var +tmp*(xti-xm) 1168 cvar=cvar+tmp*(y(j)-ym) 1169 end do 1170 1171 do 80 j=1,n 1172 out=j-ibw-1 1173 in=j+ibw 1174 if ((jper.ne.2) .and. (out.lt.1.or.in.gt.n)) go to 60 1175 if (out .lt. 1) then 1176 out=n+out 1177 xto=x(out)-1d0 1178 xti=x(in) 1179 else if (in .gt. n) then 1180 in=in-n 1181 xti=x(in)+1d0 1182 xto=x(out) 1183 else 1184 xto=x(out) 1185 xti=x(in) 1186 end if 1187 wt=w(out) 1188 fbo=fbw 1189 fbw=fbw-wt 1190 tmp=0d0 1191 if (fbw.gt.0d0) tmp=fbo*wt*(xto-xm)/fbw 1192 var = var-tmp*(xto-xm) 1193 cvar=cvar-tmp*(y(out)-ym) 1194 if (fbw.gt.0d0) xm=(fbo*xm-wt*xto)/fbw 1195 if (fbw.gt.0d0) ym=(fbo*ym-wt*y(out))/fbw 1196 wt=w(in) 1197 fbo=fbw 1198 fbw=fbw+wt 1199 if (fbw.gt.0d0) xm=(fbo*xm+wt*xti)/fbw 1200 if (fbw.gt.0d0) ym=(fbo*ym+wt*y(in))/fbw 1201 tmp=0d0 1202 if (fbo.gt.0d0) tmp=fbw*wt*(xti-xm)/fbo 1203 var = var+tmp*(xti-xm) 1204 cvar=cvar+tmp*(y(in)-ym) 1205 60 a=0d0 1206 if (var.gt.vsmlsq) a=cvar/var 1207 smo(j)=a*(x(j)-xm)+ym 1208 if (iper.gt.0) then 1209 h=0d0 1210 if (fbw.gt.0d0) h=1d0/fbw 1211 if (var.gt.vsmlsq) h=h+(x(j)-xm)**2/var 1212 acvr(j)=0d0 1213 a=1d0-w(j)*h 1214 if (a.gt.0d0) then 1215 acvr(j)=abs(y(j)-smo(j))/a 1216 else if (j.gt.1) then 1217 acvr(j)=acvr(j-1) 1218 end if 1219 end if 1220 80 continue 1221 1222 if(trace) call smoothprt(span, iper, var, cvar) ! -> ./ksmooth.c 1223 1224c-- Recompute fitted values smo(j) as weighted mean for non-unique x(.) values: 1225 j=1 1226 90 j0=j 1227 sy=smo(j)*w(j) 1228 fbw=w(j) 1229 if (j.ge.n) go to 110 1230 100 if (x(j+1).le.x(j)) then 1231 j=j+1 1232 sy=sy+w(j)*smo(j) 1233 fbw=fbw+w(j) 1234 if (j.lt.n) go to 100 1235 end if 1236 110 if (j.gt.j0) then 1237 a=0d0 1238 if (fbw.gt.0d0) a=sy/fbw 1239 do i=j0,j 1240 smo(i)=a 1241 end do 1242 end if 1243 j=j+1 1244 if (j.le.n) go to 90 1245 return 1246 end 1247 1248 1249 block data bksupsmu 1250 double precision spans(3), big,sml,eps 1251 common /spans/ spans /consts/ big,sml,eps 1252 1253 data spans, big,sml,eps /0.05,0.2,0.5, 1.0e20,1.0e-7,1.0e-3/ 1254 end 1255c--------------------------------------------------------------- 1256c 1257c this sets the compile time (default) values for various 1258c internal parameters : 1259c 1260c spans : span values for the three running linear smoothers. 1261c spans(1) : tweeter span. 1262c spans(2) : midrange span. 1263c spans(3) : woofer span. 1264c (these span values should be changed only with care.) 1265c big : a large representable floating point number. 1266c sml : a small number. should be set so that (sml)**(10.0) does 1267c not cause floating point underflow. 1268c eps : used to numerically stabilize slope calculations for 1269c running linear fits. 1270c 1271c these parameter values can be changed by declaring the 1272c relevant labeled common in the main program and resetting 1273c them with executable statements. 1274c 1275 1276c Only for ppr(*, ismethod != 0): Compute "smoothing" spline 1277c (rather, a penalized regression spline with at most 15 (inner) knots): 1278c----------------------------------------------------------------- 1279c 1280 subroutine spline (n, x, y, w, smo, edf, sc) 1281c 1282c------------------------------------------------------------------ 1283c 1284c input: 1285c n : number of observations. 1286c x(n) : ordered abscissa values. 1287c y(n) : corresponding ordinate (response) values. 1288c w(n) : weight for each (x,y) observation. 1289c work space: 1290c sc(n,7) : used for dx(n), dy(n), dw(n), dsmo(n), lev(n) 1291c output: 1292c smo(n) : smoothed ordinate (response) values. 1293c edf : equivalent degrees of freedom 1294c 1295c------------------------------------------------------------------ 1296c Args 1297 integer n 1298 double precision x(n), y(n), w(n), smo(n), edf, sc(n,7) 1299 1300 call splineAA(n, x, y, w, smo, edf, 1301 + sc(n,1), ! dx 1302 + sc(n,2), ! dy 1303 + sc(n,3), ! dw 1304 + sc(n,4), ! dsmo 1305 + sc(n,5)) ! lev 1306 1307 return 1308 end 1309 1310 1311 subroutine splineAA(n, x, y, w, smo, edf, dx, dy, dw, dsmo, lev) 1312c 1313c Workhorse of spline() above 1314c------------------------------------------------------------------ 1315c 1316c Additional input variables (no extra output, work): 1317c dx : 1318c dy : 1319c dw : 1320c dsmo: 1321c lev : "leverages", i.e., diagonal entries S_{i,i} of the smoother matrix 1322 1323c 1324c------------------------------------------------------------------ 1325c Args 1326 integer n 1327 double precision x(n), y(n), w(n), smo(n), edf, 1328 + dx(n), dy(n), dw(n), dsmo(n), lev(n) 1329c Var 1330 double precision knot(29), coef(25), work((17+25)*25) 1331 double precision param(5), df1, lambda, crit, p, s 1332 integer iparms(4), i, nk, ip, ier 1333 1334 double precision df, gcvpen 1335 integer ismethod 1336 logical trace 1337 common /spsmooth/ df, gcvpen, ismethod, trace 1338 1339c__no-more__ if (n .gt. 2500) call bdrsplerr() 1340 do i = 1,n 1341 dx(i) = (x(i)-x(1))/(x(n)-x(1)) 1342 dy(i) = y(i) 1343 dw(i) = w(i) 1344 end do 1345 nk = min(n,15) 1346 knot(1) = dx(1) 1347 knot(2) = dx(1) 1348 knot(3) = dx(1) 1349 knot(4) = dx(1) 1350 knot(nk+1) = dx(n) 1351 knot(nk+2) = dx(n) 1352 knot(nk+3) = dx(n) 1353 knot(nk+4) = dx(n) 1354 do i = 5, nk 1355 p = (n-1)*real(i-4)/real(nk-3) 1356 ip = int(p) 1357 p = p-ip 1358 knot(i) = (1-p)*dx(ip+1) + p*dx(ip+2) 1359 end do 1360c call dblepr('knots', 5, knot, nk+4) 1361C iparms(1:2) := (icrit, ispar) for ./sbart.c 1362 if (ismethod .eq. 1) then 1363 iparms(1) = 3 1364 df1 = df 1365 else 1366 iparms(1) = 1 1367 df1 = 0d0 1368 endif 1369c 1370 iparms(2) = 0 ! ispar := 0 <==> estimate `spar' 1371 iparms(3) = 500 ! maxit = 500 1372 iparms(4) = 0 ! spar (!= lambda) 1373c 1374 param(1) = 0d0 ! = lspar : min{spar} 1375 param(2) = 1.5d0 ! = uspar : max{spar} 1376c tol for 'spar' estimation: 1377 param(3) = 1d-2 1378c 'eps' (~= 2^-12 = sqrt(2^-24) ?= sqrt(machine eps)) in ./sbart.c : 1379 param(4) = .000244 1380 1381 ier = 1 1382 call rbart(gcvpen,df1,dx,dy,dw,0.0d0,n,knot,nk,coef,dsmo,lev, 1383 & crit,iparms,lambda,param, work,4,1,ier) 1384 if(ier .gt. 0) call intpr('spline(.) TROUBLE:', 18, ier, 1) 1385 do i = 1,n 1386 smo(i) = dsmo(i) 1387 end do 1388c call dblepr('smoothed',8, dsmo, n) 1389 s = 0 1390 do i = 1, n 1391 s = s + lev(i) 1392 end do 1393 edf = s 1394 if(trace) call splineprt(df,gcvpen,ismethod, lambda, edf) 1395 return 1396 end 1397 1398 1399*********************************************************************** 1400 1401C=== This was 'sort()' in gamfit's mysort.f [or sortdi() in sortdi.f ] : 1402C 1403C=== FIXME: Translate to C and add to ../../../main/sort.c <<<<< 1404C 1405C a[] is double precision because the caller reuses a double (sometimes v[] itself!) 1406 subroutine sort (v,a, ii,jj) 1407c 1408c Puts into a the permutation vector which sorts v into 1409c increasing order. Only elements from ii to jj are considered. 1410c Arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements 1411c 1412c This is a modification of CACM algorithm #347 by R. C. Singleton, 1413c which is a modified Hoare quicksort. 1414c 1415 integer ii,jj 1416 double precision v(*), a(jj) 1417c 1418 integer iu(20),il(20) 1419 integer t,tt, m,i,j,ij,k,l 1420 double precision vt, vtt 1421 1422 m=1 1423 i=ii 1424 j=jj 1425 10 if (i.ge.j) go to 80 1426 20 k=i 1427 ij=(j+i)/2 1428 t=int(a(ij)) 1429 vt=v(ij) 1430 if (v(i).le.vt) go to 30 1431 a(ij)=a(i) 1432 a(i)=t 1433 t=int(a(ij)) 1434 v(ij)=v(i) 1435 v(i)=vt 1436 vt=v(ij) 1437 30 l=j 1438 if (v(j).ge.vt) go to 50 1439 a(ij)=a(j) 1440 a(j)=t 1441 t=int(a(ij)) 1442 v(ij)=v(j) 1443 v(j)=vt 1444 vt=v(ij) 1445 if (v(i).le.vt) go to 50 1446 a(ij)=a(i) 1447 a(i)=t 1448 t=int(a(ij)) 1449 v(ij)=v(i) 1450 v(i)=vt 1451 vt=v(ij) 1452 go to 50 1453 40 a(l)=a(k) 1454 a(k)=tt 1455 v(l)=v(k) 1456 v(k)=vtt 1457 50 l=l-1 1458 if (v(l).gt.vt) go to 50 1459 tt=int(a(l)) 1460 vtt=v(l) 1461 60 k=k+1 1462 if (v(k).lt.vt) go to 60 1463 if (k.le.l) go to 40 1464 if (l-i.le.j-k) go to 70 1465 il(m)=i 1466 iu(m)=l 1467 i=k 1468 m=m+1 1469 go to 90 1470 70 il(m)=k 1471 iu(m)=j 1472 j=l 1473 m=m+1 1474 go to 90 1475 80 m=m-1 1476 if (m.eq.0) return 1477 i=il(m) 1478 j=iu(m) 1479 90 if (j-i.gt.10) go to 20 1480 if (i.eq.ii) go to 10 1481 i=i-1 1482 100 i=i+1 1483 if (i.eq.j) go to 80 1484 t=int(a(i+1)) 1485 vt=v(i+1) 1486 if (v(i).le.vt) go to 100 1487 k=i 1488 110 a(k+1)=a(k) 1489 v(k+1)=v(k) 1490 k=k-1 1491 if (vt.lt.v(k)) go to 110 1492 a(k+1)=t 1493 v(k+1)=vt 1494 go to 100 1495 end 1496