c c DART/HYESS c c This Fortran program (collection of subroutines) implements the recursive c covering approach to local learning described in c c Friedman, J. H. (1996). Local learning based on recursive covering. c (ftp://stat.stanford.edu/pub/friedman/dart.ps.Z) c c The documentation for using this program is described in c c Friedman, J. H. (1996). DART/HYESS users guide. c (ftp://stat.stanford.edu/pub/friedman/dart.doc.ps.Z) c c c Coded and copywrite (1996) by Jerome H. Friedman. c c subroutine dart (no,ni,x,y,w,mxm,mxt,itr,rtr,trm,ms) 1 integer itr(*),ms(*) 2 real x(*),y(*),w(*),rtr(*),trm(*) 3 call dart1(no,ni,x,y,w,itr,rtr,trm,ms,ms(2*no*ni+1), ms(2*ni*(no+ 4 11)+3),mxm,mxt) 5 return 6 end 7 subroutine ans (p,x,itr,rtr,trm,y) 8 integer p,itr(2,*) 9 real x(p),rtr(3,*),trm(*) 10 data lin /2/ 11 k=1 12 if(itr(1,1).eq.0.and.itr(2,1).eq.0) k=0 13 1 if(k.le.0) go to 7 14 if(rtr(1,k) .le. 0.5) go to 2 15 s=x(int(rtr(1,k)+.1)) 16 go to 4 17 2 l=abs(rtr(1,k))+.1 18 s=trm(l+p+1) 19 do 3 j=1,p 20 s=s+trm(l+j)*x(j) 21 3 continue 22 4 if(s .ge. 0.5*(rtr(2,k)+rtr(3,k))) go to 5 23 k=itr(1,k) 24 go to 1 25 5 k=itr(2,k) 26 go to 1 27 7 k=iabs(k) 28 if((lin .ne. 0) .and. (lin .ne. 2)) go to 8 29 y=trm(k+1) 30 return 31 8 s=trm(k+p+3) 32 do 9 j=1,p 33 s=s+trm(k+j+2)*x(j) 34 9 continue 35 y=s 36 return 37 entry stlin4(irg) 38 lin=irg 39 return 40 end 41 subroutine ave (p,x,itr,rtr,trm,ya) 42 parameter(mxisk=100) 43 integer p,itr(2,*),ksk(mxisk) 44 real x(p),rtr(3,*),trm(*) 45 double precision rsk(mxisk),pn,pr 46 logical cmod 47 data lin /2/ 48 isk=0 49 k=1 50 if(itr(1,1).eq.0.and.itr(2,1).eq.0) k=0 51 ya=0.0 52 wt=ya 53 pn=1.d0 54 cmod=lin.eq.0.or.lin.eq.2 55 1 if(k .le. 0) go to 9 56 isk=isk+1 57 ksk(isk)=k 58 rsk(isk)=pn 59 if(rtr(1,k) .le. 0.5) go to 2 60 s=x(int(rtr(1,k)+.1)) 61 go to 4 62 2 l=abs(rtr(1,k))+.1 63 s=trm(l+p+1) 64 do 3 j=1,p 65 s=s+trm(l+j)*x(j) 66 3 continue 67 4 s1=rtr(2,k) 68 s2=rtr(3,k) 69 d=s2-s1 70 if(d .gt. 0.0) go to 5 71 pr=0.d0 72 if(s.gt.s1) pr=1.d0 73 go to 6 74 5 pr=amax1(0.0,amin1(1.0,(s-s1)/d)) 75 6 if(s .ge. s2) go to 7 76 pn=pn*(1.d0-pr) 77 k=itr(1,k) 78 go to 1 79 7 pn=pn*pr 80 k=itr(2,k) 81 go to 1 82 9 l=iabs(k) 83 if(.not.(cmod)) go to 10 84 s=trm(l+1) 85 go to 12 86 10 s=trm(l+p+3) 87 do 11 j=1,p 88 s=s+trm(l+j+2)*x(j) 89 11 continue 90 12 wt=wt+pn 91 ya=ya+pn*s 92 13 if(isk.eq.0) go to 19 93 if(itr(1,ksk(isk)) .ne. k) go to 17 94 kp=ksk(isk) 95 if(rtr(1,kp) .le. 0.5) go to 14 96 s=x(int(rtr(1,kp)+.1)) 97 go to 16 98 14 l=abs(rtr(1,kp))+.1 99 s=trm(l+p+1) 100 do 15 j=1,p 101 s=s+trm(l+j)*x(j) 102 15 continue 103 16 if(s.gt.rtr(2,kp)) go to 19 104 k=kp 105 isk=isk-1 106 go to 13 107 17 k=ksk(isk) 108 isk=isk-1 109 go to 13 110 19 if(isk.eq.0) go to 22 111 s1=rtr(2,kp) 112 s2=rtr(3,kp) 113 d=s2-s1 114 if(d .gt. 0.0) go to 20 115 pr=0.d0 116 if(s.gt.s1) pr=1.d0 117 go to 21 118 20 pr=amax1(0.0,amin1(1.0,(s-s1)/d)) 119 21 pn=rsk(isk)*pr 120 k=itr(2,kp) 121 go to 1 122 22 ya=ya/wt 123 return 124 entry stlin3(irg) 125 lin=irg 126 return 127 end 128 subroutine start (no,ni,x,y,w,itr,rtr,itrm,rtrm,dtrm,mxm,mxi,mxr,m 129 1xd,ms) 130 integer itr(*),itrm(*),ms(*) 131 real x(*),y(*),w(*),rtr(*),rtrm(*) 132 double precision dtrm(*) 133 call start1(no,ni,x,y,w,itr,rtr,itrm,rtrm,dtrm,ms,ms(2*no*ni+1), 134 1ms(2*ni*(no+1)+3),mxm,mxi,mxr,mxd) 135 return 136 end 137 subroutine hyess(xp,n,p,x,y,w,itr,rtr,itrm,rtrm,dtrm,yh,mt) 138 parameter (mxit=400) 139 integer p,itr(*),itrm(*),mt(*) 140 real xp(p+1),x(n,p+3),y(n),w(n),rtr(*),rtrm(*),yp(mxit),fn(mxit) 141 double precision dtrm(*),s,esqm,drg 142 data lin,ntn,esqm /2,10,1.d-6/ 143 call ktrm(p,xp,itr,rtr,rtrm,ii) 144 ir=itrm(ii) 145 id=itrm(ii+1) 146 ii=ii+2 147 yp(1)=rtrm(ir) 148 i3=ii+2*p+2 149 if((itrm(i3)-itrm(i3+1) .ge. ntn) .and. (dtrm(id+2) .ge. esqm)) go 150 1 to 3 151 if((lin .ne. 0) .and. (lin .ne. 2)) go to 1 152 it=1 153 return 154 1 il=id+itrm(i3+5)-1 155 s=dtrm(il+p+1) 156 do 2 j=1,p 157 s=s+dtrm(il+j)*xp(j) 158 2 continue 159 yp(1)=s 160 it=1 161 return 162 3 nt=itrm(ii+2*p) 163 i1=ii+2*p+8 164 i2=i1+2*nt*p 165 i4=i2+nt 166 fn(1)=nt 167 call hyess1 (xp,0,n,p,x,y,w,nt,itrm(i1),itrm(ii),itrm(i4),itrm(i2) 168 1, rtrm(ir+1),dtrm(id),yp,fn,it,mt,mt(2*n*(p+1)+1),mt(2*n*(p+1)+2* 169 1p+5)) 170 if(it.gt.mxit) call diag(14,it,mxit) 171 yh=yp(it) 172 return 173 entry stlin8(irg) 174 lin=irg 175 return 176 entry sthyn1(irg) 177 ntn=irg 178 return 179 entry stesq2(drg) 180 esqm=drg 181 return 182 end 183 subroutine xvalid (n,p,x,y,w,itr,rtr,itrm,rtrm,dtrm,nks,erm,err,ln 184 1,mt) 185 parameter (maxdim=100, mxit=400) 186 integer p,itr(*),itrm(*),mt(*) 187 real xp(maxdim+1),x(n,p+3),y(n),w(n),rtr(*),rtrm(*),err(2,*) 188 real yp(mxit),fn(mxit),er2(mxit) 189 double precision dtrm(*) 190 data lnx,ntn,big,jx,thr /40,10,9.9e30,1,-1.0/ 191 if(jx.gt.n/2) call diag(15,jx,n/2) 192 ln=1 193 err(1,1)=ntn 194 ftn=(float(n-1)/ntn)**(1.0/(lnx-1)) 195 1 if(err(1,ln).gt.n-1) go to 2 196 ln=ln+1 197 err(1,ln)=ftn*err(1,ln-1) 198 go to 1 199 2 ln=min0(ln,lnx) 200 fnx=n-1 201 err(1,ln)=fnx 202 do 3 i=1,ln 203 err(2,i)=0.0 204 er2(i)=err(2,i) 205 3 continue 206 ix=0 207 4 ix=ix+jx 208 if(ix.gt.n) go to 8 209 do 5 j=1,p 210 xp(j)=x(ix,j) 211 5 continue 212 call ktrm(p,xp,itr,rtr,rtrm,ii) 213 ir=itrm(ii) 214 id=itrm(ii+1) 215 ii=ii+2 216 yp(1)=rtrm(ir) 217 nt=itrm(ii+2*p) 218 i1=ii+2*p+8 219 i2=i1+2*nt*p 220 i4=i2+nt 221 fn(1)=nt-1 222 fnx=amin1(fnx,fn(1)) 223 call hyess1 (xp,ix,n,p,x,y,w,nt,itrm(i1),itrm(ii),itrm(i4),itrm(i2 224 1), rtrm(ir+1),dtrm(id),yp,fn,it,mt,mt(2*n*(p+1)+1),mt(2*n*(p+1)+2 225 1*p+5)) 226 if(it.gt.mxit) call diag(14,it,mxit) 227 k=1 228 do 7 i=it,1,-1 229 yc=yp(i) 230 if(thr.gt.0.0) yc=sign(0.5,yp(i)-thr)+0.5 231 eri=w(ix)*abs(y(ix)-yc) 232 6 if(err(1,k).gt.fn(i)) go to 7 233 err(2,k)=err(2,k)+eri 234 er2(k)=er2(k)+w(ix) 235 k=k+1 236 if(k.gt.ln) go to 7 237 go to 6 238 7 continue 239 go to 4 240 8 erm=big 241 nkm=ntn 242 do 10 i=1,ln 243 if((er2(i) .gt. 0.0) .and. (err(1,i) .le. fnx)) go to 9 244 err(2,i)=big 245 go to 10 246 9 err(2,i)=err(2,i)/er2(i) 247 if(err(2,i) .gt. erm) go to 10 248 erm=err(2,i) 249 nks=err(1,i)+0.5 250 10 continue 251 return 252 entry setnxv(irg) 253 lnx=min0(mxit,max0(1,irg)) 254 return 255 entry setkxv(irg) 256 jx=max0(1,irg) 257 return 258 entry setcls(arg) 259 thr=arg 260 return 261 entry sthyn(irg) 262 ntn=irg 263 return 264 end 265 subroutine dart1 (n,p,x,y,w,itr,rtr,trm,m,mp,ll,mxit,mxitm) 266 parameter (iflg=999999, mxisk=100, mxko=100000, maxid=50000) 267 integer p,itr(2,mxit),m(2,n,p),mp(2,p+1),ll(n),ksk(4,mxisk),osk(mx 268 1ko) 269 real x(n,p),y(n),w(n),rtr(*),trm(*) 270 double precision dp(maxid) 271 save it,itm 272 it=0 273 isk=it 274 itm=isk 275 ko=1 276 id=ko 277 ls=id 278 ksk(2,1)=ls 279 ksk(4,1)=ksk(2,1) 280 call strt(0) 281 call initl(n,p,x,y,w,m,mp,dp,nd,ll) 282 if(nd.gt.maxid) call diag(5,nd,maxid) 283 1 call split(n,p,x,y,w,dp(id),m,mp,it+1,itm,rtr,trm(itm+1),nt,ll,m1, 284 1m2) 285 if(m1 .le. 0) go to 4 286 it=it+1 287 if(it.ge.mxit) call diag(1,it,mxit) 288 itm=itm+nt 289 if(itm.ge.mxitm) call diag(4,itm,mxitm) 290 if(isk .le. 0) go to 2 291 itr(2,ksk(1,isk))=iflg 292 itr(ls,ksk(1,isk))=it 293 2 isk=isk+1 294 if(isk.gt.mxisk) call diag(2,isk,mxisk) 295 ksk(1,isk)=it 296 ksk(3,isk)=m1 297 do 3 i=1,m2 298 osk(ko)=ll(i) 299 ko=ko+1 300 if(ko.gt.mxko) call diag(3,ko,mxko) 301 3 continue 302 ksk(2,isk+1)=ko 303 call delete(n,p,ll(m1+1),m2-m1,m,mp) 304 call mdldel(n,p,x,y,w,dp(id),ll(m1+1),m2-m1,dp(id+nd)) 305 id=id+nd 306 ksk(4,isk+1)=id 307 if(id.gt.maxid) call diag(5,id,maxid) 308 ls=1 309 go to 1 310 4 if(isk .gt. 0) go to 5 311 itr(1,1)=0 312 itr(2,1)=itr(1,1) 313 go to 9 314 5 itr(2,ksk(1,isk))=iflg 315 itr(ls,ksk(1,isk))=-itm 316 itm=itm+nt 317 if(itm.ge.mxitm) call diag(4,itm,mxitm) 318 m1=ksk(3,isk) 319 ko=ksk(2,isk) 320 m2=ksk(2,isk+1)-ko 321 id=ksk(4,isk) 322 6 if(itr(2,ksk(1,isk)).eq.iflg) go to 7 323 call insert(n,p,osk(ko),m1,m,mp) 324 isk=isk-1 325 m1=ksk(3,isk) 326 ko=ksk(2,isk) 327 m2=ksk(2,isk+1)-ko 328 id=ksk(4,isk) 329 if(isk.le.0) go to 7 330 go to 6 331 7 if(isk.le.0) go to 9 332 call insert(n,p,osk(ko+m1),m2-m1,m,mp) 333 call delete(n,p,osk(ko),m1,m,mp) 334 call mdldel(n,p,x,y,w,dp(id),osk(ko),m1,dp(id+nd)) 335 ko=ksk(2,isk+1) 336 id=ksk(4,isk+1) 337 ls=2 338 go to 1 339 9 return 340 entry dtsize(irg1,irg2) 341 irg1=it+1 342 irg2=itm 343 return 344 end 345 subroutine start1 (n,p,x,y,w,itr,rtr,itrm,rtrm,dtrm,m,mp,ll,mxit,m 346 1xi,mxr,mxd) 347 parameter (iflg=999999, mxisk=100, mxko=100000, maxid=50000) 348 integer p,itr(2,mxit),m(2,n,p),mp(2,p+1),ll(n),ksk(4,mxisk),osk(mx 349 1ko),itrm(mxi) 350 real x(n,p+2),y(n),w(n),rtr(*),rtrm(mxr) 351 double precision dtrm(mxd),dp(maxid) 352 save it,ii,ir,jd 353 it=0 354 isk=it 355 ko=1 356 id=ko 357 ls=id 358 ksk(2,1)=ls 359 ksk(4,1)=ksk(2,1) 360 ii=ksk(4,1) 361 ir=ii 362 jd=ir 363 call strt(1) 364 call initl(n,p,x,y,w,m,mp,dp,nd,ll) 365 if(nd.gt.maxid) call diag(5,nd,maxid) 366 1 call split(n,p,x,y,w,dp(id),m,mp,it+1,ir-1,rtr,rtrm(ir),nt,ll,m1,m 367 12) 368 if(m1 .le. 0) go to 4 369 it=it+1 370 if(it.ge.mxit) call diag(1,it,mxit) 371 ir=ir+nt 372 if(ir.ge.mxr) call diag(11,ir,mxr) 373 if(isk .le. 0) go to 2 374 itr(2,ksk(1,isk))=iflg 375 itr(ls,ksk(1,isk))=it 376 2 isk=isk+1 377 if(isk.gt.mxisk) call diag(2,isk,mxisk) 378 ksk(1,isk)=it 379 ksk(3,isk)=m1 380 do 3 i=1,m2 381 osk(ko)=ll(i) 382 ko=ko+1 383 if(ko.gt.mxko) call diag(3,ko,mxko) 384 3 continue 385 ksk(2,isk+1)=ko 386 call delete(n,p,ll(m1+1),m2-m1,m,mp) 387 call mdldel(n,p,x,y,w,dp(id),ll(m1+1),m2-m1,dp(id+nd)) 388 id=id+nd 389 ksk(4,isk+1)=id 390 if(id.gt.maxid) call diag(5,id,maxid) 391 ls=1 392 go to 1 393 4 itr(2,ksk(1,isk))=iflg 394 itr(ls,ksk(1,isk))=-ii 395 call trmnl(n,p,x,y,w,m,mp,dp(id),ii,ir,jd,itrm,rtrm,dtrm,ll) 396 if(ii.ge.mxi) call diag(12,ii,mxi) 397 if(ir.ge.mxr) call diag(11,ir,mxr) 398 if(jd.ge.mxd) call diag(13,jd,mxd) 399 if(isk .gt. 0) go to 5 400 itr(1,1)=0 401 itr(2,1)=itr(1,1) 402 go to 9 403 5 m1=ksk(3,isk) 404 ko=ksk(2,isk) 405 m2=ksk(2,isk+1)-ko 406 id=ksk(4,isk) 407 6 if(itr(2,ksk(1,isk)).eq.iflg) go to 7 408 call insert(n,p,osk(ko),m1,m,mp) 409 isk=isk-1 410 m1=ksk(3,isk) 411 ko=ksk(2,isk) 412 m2=ksk(2,isk+1)-ko 413 id=ksk(4,isk) 414 if(isk.le.0) go to 7 415 go to 6 416 7 if(isk.le.0) go to 9 417 call insert(n,p,osk(ko+m1),m2-m1,m,mp) 418 call delete(n,p,osk(ko),m1,m,mp) 419 call mdldel(n,p,x,y,w,dp(id),osk(ko),m1,dp(id+nd)) 420 ko=ksk(2,isk+1) 421 id=ksk(4,isk+1) 422 ls=2 423 go to 1 424 9 return 425 entry stsize(irg1,irg2,irg3,irg4) 426 irg1=it+1 427 irg2=ii 428 irg3=ir 429 irg4=jd 430 return 431 end 432 subroutine hyess1 (xp,ix,n,p,x,y,w,nt,mm,mmp,lll,mq,sp,dpp,yp,fn,i 433 1t,m,mp,ll) 434 parameter(maxdim=100) 435 integer p,pp,mm(2,nt,p),mmp(2,p+4),lll(nt),mq(nt),m(2,nt,p+1),mp(2 436 1,p+2),ll(nt) 437 real xp(p+1),x(n,p+3),y(n),w(n),sp(nt),yp(*),fn(*) 438 double precision dpp(*),dp((maxdim+2)*(maxdim+1)+3),a(maxdim+1) 439 data kn /0/ 440 m1=mmp(2,p+2) 441 do 2 j=1,p 442 mp(1,j)=mmp(1,j) 443 mp(2,j)=mmp(2,j) 444 do 1 i=1,nt 445 m(1,i,j)=mm(1,i,j) 446 m(2,i,j)=mm(2,i,j) 447 1 continue 448 2 continue 449 mp(1,p+1)=mmp(1,p+1) 450 mp(2,p+1)=mmp(2,p+1) 451 it=1 452 lx=0 453 pp=p 454 if(kn .le. 0) go to 3 455 call knn(xp,n,p,x,nt,m,mp,mq,ll) 456 pp=p+1 457 3 if(ix .le. 0) go to 4 458 call mdldel(n,p,x,y,w,dpp,ix,1,dp) 459 yp(1)=dp(1)/dp(2) 460 call bfnd(nt,ix,x(ix,1),mq,x,lx) 461 if(lx.le.0) call diag(16,lx,0) 462 call delete(nt,pp,lx,1,m,mp) 463 4 if(mmp(1,p+3) .gt. 0) go to 6 464 nx=0 465 do 5 i=1,nt 466 if(i.eq.lx) go to 5 467 nx=nx+1 468 x(nx,pp+2)=sp(i) 469 ll(nx)=lll(i) 470 5 continue 471 6 call seldel(xp,n,p,x,nx,m1, mmp(1,p+3),mmp(2,p+3),mmp(1,p+4),dpp( 472 1mmp(2,p+4)+1),nt,m,mp,mq,ll,yp(1)) 473 if(m1.le.0) return 474 call delete(nt,pp,ll,m1,m,mp) 475 do 7 i=1,m1 476 ll(i)=mq(ll(i)) 477 7 continue 478 if(ix .ne. 0) go to 8 479 call mdldel(n,p,x,y,w,dpp,ll,m1,dp) 480 go to 9 481 8 call mdldel(n,p,x,y,w,dp,ll,m1,dp) 482 9 it=it+1 483 call cut(n,p,x,y,w,dp,nt,m,mp,mq,ll,nu,m1,js,j1,j2,a,yp(it)) 484 call seldel(xp,n,p,x,nu,m1,js,j1,j2,a,nt,m,mp,mq,ll,yp(it)) 485 fn(it)=nu 486 if(m1.le.0) go to 11 487 call delete(nt,pp,ll,m1,m,mp) 488 do 10 i=1,m1 489 ll(i)=mq(ll(i)) 490 10 continue 491 call mdldel(n,p,x,y,w,dp,ll,m1,dp) 492 go to 9 493 11 return 494 entry setknn(irg) 495 kn=irg 496 call stkn1(kn) 497 call stkn2(kn) 498 call stkn3(kn) 499 return 500 end 501 subroutine split (n,p,x,y,w,dp,m,mp,jt,itm,rtr,trm,nt,ll,m1,m2) 502 parameter (maxdim=100) 503 integer p,m(2,n,p),mp(2,p+1),ll(*),ma(2) 504 real x(n,p+2),y(n),w(n),rtr(3,*),trm(*) 505 double precision dp(*),s,esqm,drg,a(maxdim+1),ea(2),ym 506 logical term,lmod,cmod 507 data thl,lin,trf,ntn,esqm,istr /1.e-6,2,0.25,10,1.d-6,0/ 508 nu=mp(1,p+1) 509 m1=min0(nu/2,max0(1,int(trf*nu+0.5))) 510 ym=dp(1)/dp(2) 511 crm=0.0 512 js=1 513 term=nu-m1.lt.ntn.or.dp(3).lt.esqm 514 lmod=lin.eq.1.or.lin.eq.3 515 cmod=.not.lmod 516 if(lin .le. 0) go to 2 517 call soln(p+1,dp(4),dp((p+1)**2+4),a) 518 if(term) go to 12 519 call search(2,1,n,m,mp,nu,m2,ll) 520 call proj(n,p,x,m2,ll,a,x(1,p+2)) 521 if(.not.(lmod)) go to 2 522 s=0.d0 523 do 1 i=1,m2 524 k=ll(i) 525 x(k,p+1)=y(k)-x(k,p+2) 526 s=s+w(k)*x(k,p+1)**2 527 1 continue 528 if(s.lt.thl*dp(3)) term=.true. 529 2 if(term) go to 12 530 if(lin .lt. 2) go to 4 531 call psort(x(1,p+2),ll,1,m2) 532 if(x(ll(m2),p+2) .le. x(ll(1),p+2)) go to 4 533 js=0 534 if(.not.(cmod)) go to 3 535 call eavl(n,x(1,p+2),y,w,m1,m2,ll,ea,ma) 536 crm=abs(ea(1)-ym)+abs(ea(2)-ym) 537 j1=ma(1) 538 j2=ma(2) 539 go to 4 540 3 call eavl(n,x(1,p+2),x(1,p+1),w,m1,m2,ll,ea,ma) 541 crm=abs(ea(1))+abs(ea(2)) 542 j1=ma(1) 543 j2=ma(2) 544 4 do 7 j=1,p 545 if(x(mp(1,j),j).le.x(mp(2,j),j)) go to 7 546 if(.not.(cmod)) go to 5 547 call eav(j,n,x(1,j),y,w,m1,m,mp,ea,ma) 548 cr=abs(ea(1)-ym)+abs(ea(2)-ym) 549 go to 6 550 5 call eav(j,n,x(1,j),x(1,p+1),w,m1,m,mp,ea,ma) 551 cr=abs(ea(1))+abs(ea(2)) 552 6 if(cr .lt. crm) go to 7 553 crm=cr 554 js=j 555 j1=ma(1) 556 j2=ma(2) 557 7 continue 558 if(crm.le.0.0) go to 12 559 if(j1 .le. nu-j2) go to 8 560 j=j1 561 j1=nu-j2 562 j2=nu-j 563 8 if(js .le. 0) go to 9 564 rtr(1,jt)=js 565 call search(2,js,n,m,mp,j1,k,ll) 566 rtr(2,jt)=x(ll(k),js) 567 call search(1,js,n,m,mp,j2,m2,ll(k+1)) 568 m1=k 569 m2=m1+m2 570 rtr(3,jt)=x(ll(m2),js) 571 nt=0 572 return 573 9 nh=m2-j2 574 rtr(1,jt)=-itm 575 nt=p+1 576 do 10 j=1,nt 577 trm(j)=a(j) 578 10 continue 579 rtr(2,jt)=x(ll(j1),p+2) 580 rtr(3,jt)=x(ll(nh+1),p+2) 581 do 11 i=1,j2 582 ll(i+j1)=ll(i+nh) 583 11 continue 584 m1=j1 585 m2=j1+j2 586 return 587 12 m1=0 588 m2=m1 589 if(istr.eq.1) return 590 trm(1)=ym 591 trm(2)=dp(3)/dp(2) 592 nt=2 593 if(.not.(lmod)) go to 14 594 do 13 j=1,p+1 595 trm(j+2)=a(j) 596 13 continue 597 nt=nt+p+1 598 14 return 599 entry setdtf(arg) 600 trf=arg 601 return 602 entry setthl(arg) 603 thl=arg 604 call stthl(thl) 605 return 606 entry setdtn(irg) 607 ntn=max0(irg,3) 608 return 609 entry stesqm(drg) 610 esqm=drg 611 return 612 entry stlin2(irg) 613 lin=irg 614 return 615 entry strt(irg) 616 istr=irg 617 return 618 end 619 subroutine eav (j,n,x,y,w,na,m,mp,ea,ma) 620 integer m(2,n,*),mp(2,*),ma(2) 621 real x(n),y(n),w(n) 622 double precision ea(2),s,t,ss,tt 623 do 5 ms=1,2 624 ks=mod(ms,2)+1 625 k=mp(ks,j) 626 l=0 627 s=0.d0 628 t=s 629 xs=x(k) 630 1 l=l+1 631 s=s+w(k)*y(k) 632 t=t+w(k) 633 k=m(ks,k,j) 634 if(k.eq.0) go to 3 635 if(x(k) .eq. xs) go to 1 636 ls=l 637 ss=s 638 tt=t 639 xs=x(k) 640 if(l.ge.na) go to 3 641 go to 1 642 3 if((k .ne. 0) .and. (iabs(ls-na) .gt. iabs(l-na))) go to 4 643 ea(ms)=ss/tt 644 ma(ms)=ls 645 go to 5 646 4 ea(ms)=s/t 647 ma(ms)=l 648 5 continue 649 return 650 end 651 subroutine eavl (n,x,y,w,na,nt,ll,ea,ma) 652 integer ll(nt),ma(2) 653 real x(n),y(n),w(n) 654 double precision ea(2),s,t,ss,tt 655 do 8 ms=1,2 656 if(ms .ne. 1) go to 1 657 ks=1 658 js=ks 659 go to 2 660 1 ks=nt 661 js=-1 662 2 l=0 663 s=0.d0 664 t=s 665 k=ll(ks) 666 xs=x(k) 667 3 l=l+1 668 s=s+w(k)*y(k) 669 t=t+w(k) 670 ks=ks+js 671 if((ks .ge. 1) .and. (ks .le. nt)) go to 4 672 k=0 673 go to 6 674 4 k=ll(ks) 675 if(x(k) .eq. xs) go to 3 676 ls=l 677 ss=s 678 tt=t 679 xs=x(k) 680 if(l.ge.na) go to 6 681 go to 3 682 6 if((k .ne. 0) .and. (iabs(ls-na) .gt. iabs(l-na))) go to 7 683 ea(ms)=ss/tt 684 ma(ms)=ls 685 go to 8 686 7 ea(ms)=s/t 687 ma(ms)=l 688 8 continue 689 return 690 end 691 subroutine proj (n,p,x,nl,ll,a,z) 692 integer p,ll(nl) 693 real x(n,p),z(n) 694 double precision a(p+1),s 695 do 2 i=1,nl 696 s=a(p+1) 697 do 1 j=1,p 698 s=s+a(j)*x(ll(i),j) 699 1 continue 700 z(ll(i))=s 701 2 continue 702 return 703 end 704 subroutine mdlbgn (n,p,x,y,w,nd,dp) 705 parameter (maxdim=100) 706 integer p 707 real x(n,p),y(n),w(n) 708 double precision dp(*),s 709 double precision ds(maxdim+1,maxdim+1) 710 data thr,lin /1.e-6,2/ 711 nd=3 712 do 1 i=1,3 713 dp(i)=0.d0 714 1 continue 715 do 2 i=1,n 716 dp(1)=dp(1)+w(i)*y(i) 717 dp(2)=dp(2)+w(i) 718 2 continue 719 if(dp(2).le.0.d0) call diag(6,int(dp(2)),0) 720 s=dp(1)/dp(2) 721 do 3 i=1,n 722 dp(3)=dp(3)+w(i)*(y(i)-s)**2 723 3 continue 724 call stesqm(thr*dp(3)) 725 call stesq1(thr*dp(3)) 726 if(lin.le.0) return 727 call xtxty(n,p,x,y,w,ds,dp((p+1)**2+4)) 728 call inv(p+1,ds,dp(4),ie) 729 if(ie.ne.0) call diag(8,ie,0) 730 nd=nd+(p+2)*(p+1) 731 return 732 entry setthr(arg) 733 thr=arg 734 return 735 entry setlin(irg) 736 lin=irg 737 call stlin1(lin) 738 call stlin2(lin) 739 call stlin3(lin) 740 call stlin4(lin) 741 call stlin5(lin) 742 call stlin6(lin) 743 return 744 end 745 subroutine mdldel (n,p,x,y,w,dp,ld,nd,dq) 746 parameter (maxdim=100) 747 integer p,ld(nd) 748 real x(n,p),y(n),w(n) 749 double precision dp(*),dq(*),s,s0,dl(maxdim+1),ds(maxdim+1),st,yp 750 data lin /2/ 751 do 1 i=1,3 752 dq(i)=dp(i) 753 1 continue 754 if(lin .le. 0) go to 3 755 nu=3+(p+2)*(p+1) 756 do 2 i=4,nu 757 dq(i)=dp(i) 758 2 continue 759 3 do 5 i=1,nd 760 l=ld(i) 761 s0=dq(2) 762 s=dq(1)/dq(2) 763 dq(2)=dq(2)-w(l) 764 if(dq(2).le.0.d0) call diag(6,int(dq(2)),0) 765 dq(1)=dq(1)-w(l)*y(l) 766 dq(3)=dq(3)-s0*w(l)*(y(l)-s)**2/dq(2) 767 if(lin.le.0) go to 5 768 st=1.d0 769 if(w(l).ne.1.0) st=sqrt(w(l)) 770 do 4 j=1,p 771 dl(j)=st*x(l,j) 772 4 continue 773 dl(p+1)=st 774 yp=st*y(l) 775 call dd(p+1,dq(4),dq((p+1)**2+4),dl,yp,ds,ie) 776 if(ie.ne.0) call diag(9,ie,0) 777 5 continue 778 return 779 entry stlin1(irg) 780 lin=irg 781 return 782 end 783 subroutine llist (n,p,x,w,m,mi,mp) 784 integer p,m(2,n,p),mi(n),mp(2,p+1) 785 real x(n,p),w(n) 786 do 3 j=1,p 787 do 1 i=1,n 788 mi(i)=i 789 1 continue 790 call psort(x(1,j),mi,1,n) 791 mp(2,j)=mi(1) 792 mp(1,j)=mi(n) 793 m(1,mi(1),j)=0 794 m(2,mi(n),j)=m(1,mi(1),j) 795 m(2,mi(1),j)=mi(2) 796 m(1,mi(n),j)=mi(n-1) 797 do 2 i=2,n-1 798 m(1,mi(i),j)=mi(i-1) 799 m(2,mi(i),j)=mi(i+1) 800 2 continue 801 3 continue 802 mp(1,p+1)=n 803 k=0 804 do 4 i=1,n 805 if(w(i).gt.0.0) go to 4 806 k=k+1 807 mi(k)=i 808 4 continue 809 if(k.gt.0) call delete(n,p,mi,k,m,mp) 810 return 811 end 812 subroutine delete (n,p,ld,nd,m,mp) 813 integer p,ld(nd),m(2,n,p),mp(2,p+1) 814 do 5 i=1,nd 815 mp(1,p+1)=mp(1,p+1)-1 816 l=ld(i) 817 do 4 j=1,p 818 k1=m(1,l,j) 819 k2=m(2,l,j) 820 if(k1 .le. 0) go to 1 821 m(2,k1,j)=k2 822 go to 2 823 1 mp(2,j)=k2 824 2 if(k2 .le. 0) go to 3 825 m(1,k2,j)=k1 826 go to 4 827 3 mp(1,j)=k1 828 4 continue 829 5 continue 830 return 831 end 832 subroutine insert (n,p,li,ni,m,mp) 833 integer p,li(ni),m(2,n,p),mp(2,p+1) 834 do 6 i=ni,1,-1 835 mp(1,p+1)=mp(1,p+1)+1 836 l=li(i) 837 do 5 j=1,p 838 if(mp(1,p+1) .ne. 1) go to 1 839 m(1,l,j)=0 840 m(2,l,j)=m(1,l,j) 841 mp(1,j)=l 842 mp(2,j)=mp(1,j) 843 go to 5 844 1 k1=m(1,l,j) 845 k2=m(2,l,j) 846 if(k1 .le. 0) go to 2 847 m(2,k1,j)=l 848 go to 3 849 2 mp(2,j)=l 850 3 if(k2 .le. 0) go to 4 851 m(1,k2,j)=l 852 go to 5 853 4 mp(1,j)=l 854 5 continue 855 6 continue 856 return 857 end 858 subroutine search (id,jp,n,m,mp,ns,nl,ls) 859 integer m(2,n,*),mp(2,*),ls(ns) 860 k=mp(id,jp) 861 nl=0 862 1 if(nl.ge.ns) go to 2 863 nl=nl+1 864 ls(nl)=k 865 if(m(id,k,jp).eq.0) go to 2 866 k=m(id,k,jp) 867 go to 1 868 2 return 869 end 870 subroutine eav1 (j,n,x,y,w,na,m,mp,mq,ea,ma) 871 integer m(2,n,*),mp(2,*),mq(n),ma(2) 872 real x(*),y(*),w(*) 873 double precision ea(2),s,t,ss,tt 874 do 5 ms=1,2 875 ks=mod(ms,2)+1 876 k=mp(ks,j) 877 l=0 878 s=0.d0 879 t=s 880 kk=mq(k) 881 xs=x(kk) 882 1 l=l+1 883 s=s+w(kk)*y(kk) 884 t=t+w(kk) 885 k=m(ks,k,j) 886 if(k.eq.0) go to 3 887 kk=mq(k) 888 if(x(kk) .eq. xs) go to 1 889 ls=l 890 ss=s 891 tt=t 892 xs=x(kk) 893 if(l.ge.na) go to 3 894 go to 1 895 3 if((k .ne. 0) .and. (iabs(ls-na) .gt. iabs(l-na))) go to 4 896 ea(ms)=ss/tt 897 ma(ms)=ls 898 go to 5 899 4 ea(ms)=s/t 900 ma(ms)=l 901 5 continue 902 return 903 end 904 subroutine eavl1 (n,x,y,w,na,nt,ll,mq,ea,ma) 905 integer ll(nt),mq(n),ma(2) 906 real x(nt),y(*),w(*) 907 double precision ea(2),s,t,ss,tt 908 do 8 ms=1,2 909 if(ms .ne. 1) go to 1 910 ks=1 911 js=ks 912 go to 2 913 1 ks=nt 914 js=-1 915 2 l=0 916 s=0.d0 917 t=s 918 k=ll(ks) 919 kk=mq(k) 920 xs=x(k) 921 3 l=l+1 922 s=s+w(kk)*y(kk) 923 t=t+w(kk) 924 ks=ks+js 925 if((ks .ge. 1) .and. (ks .le. nt)) go to 4 926 k=0 927 go to 6 928 4 k=ll(ks) 929 kk=mq(k) 930 if(x(k) .eq. xs) go to 3 931 ls=l 932 ss=s 933 tt=t 934 xs=x(k) 935 if(l.ge.na) go to 6 936 go to 3 937 6 if((k .ne. 0) .and. (iabs(ls-na) .gt. iabs(l-na))) go to 7 938 ea(ms)=ss/tt 939 ma(ms)=ls 940 go to 8 941 7 ea(ms)=s/t 942 ma(ms)=l 943 8 continue 944 return 945 end 946 subroutine ktrm (p,x,itr,rtr,trm,k) 947 integer p,itr(2,*) 948 real x(p),rtr(3,*),trm(*) 949 k=1 950 if(itr(1,1).eq.0.and.itr(2,1).eq.0) return 951 1 if(k.le.0) go to 7 952 if(rtr(1,k) .le. 0.5) go to 2 953 s=x(int(rtr(1,k)+.1)) 954 go to 4 955 2 l=abs(rtr(1,k))+.1 956 s=trm(l+p+1) 957 do 3 j=1,p 958 s=s+trm(l+j)*x(j) 959 3 continue 960 4 if(s .ge. 0.5*(rtr(2,k)+rtr(3,k))) go to 5 961 k=itr(1,k) 962 go to 1 963 5 k=itr(2,k) 964 go to 1 965 7 k=iabs(k) 966 return 967 end 968 subroutine proj1 (n,p,x,nl,ll,mq,a,z) 969 integer p,ll(nl),mq(n) 970 real x(n,p),z(nl) 971 double precision a(p+1),s 972 do 2 i=1,nl 973 s=a(p+1) 974 k=mq(ll(i)) 975 do 1 j=1,p 976 s=s+a(j)*x(k,j) 977 1 continue 978 z(ll(i))=s 979 2 continue 980 return 981 end 982 subroutine cut (n,p,x,y,w,dp,nt,m,mp,mq,ll,nu,m1,js,j1,j2,a,yp) 983 integer p,pp,m(2,nt,p+1),mp(2,p+2),ll(nt),mq(nt),ma(2) 984 real x(n,p+3),y(n),w(n) 985 double precision dp(*),s,esqm,drg,a(p+1),ea(2),ym 986 logical term,lmod,cmod 987 data thl,lin,trf,ntn,esqm,kn /1.e-6,2,0.1,10,1.d-6,0/ 988 pp=p 989 if(kn.gt.0) pp=p+1 990 nu=mp(1,pp+1) 991 m1=min0(nu/2,max0(1,int(trf*nu+0.5))) 992 ym=dp(1)/dp(2) 993 yp=ym 994 crm=0.0 995 js=1 996 term=nu-m1.lt.ntn.or.dp(3).lt.esqm 997 lmod=lin.eq.1.or.lin.eq.3 998 cmod=.not.lmod 999 if(lin .le. 0) go to 2 1000 call soln(p+1,dp(4),dp((p+1)**2+4),a) 1001 if(term) go to 9 1002 call search(2,1,nt,m,mp,nu,m2,ll) 1003 call proj1(n,p,x,m2,ll,mq,a,x(1,pp+2)) 1004 if(.not.(lmod)) go to 2 1005 s=0.d0 1006 do 1 i=1,m2 1007 k=ll(i) 1008 kk=mq(k) 1009 x(kk,pp+1)=y(kk)-x(k,pp+2) 1010 s=s+w(kk)*x(kk,pp+1)**2 1011 1 continue 1012 if(s.lt.thl*dp(3)) term=.true. 1013 2 if(term) go to 9 1014 if(lin .lt. 2) go to 4 1015 call psort(x(1,pp+2),ll,1,m2) 1016 if(x(ll(m2),pp+2) .le. x(ll(1),pp+2)) go to 4 1017 js=0 1018 if(.not.(cmod)) go to 3 1019 call eavl1(nt,x(1,pp+2),y,w,m1,m2,ll,mq,ea,ma) 1020 crm=abs(ea(1)-ym)+abs(ea(2)-ym) 1021 j1=ma(1) 1022 j2=ma(2) 1023 go to 4 1024 3 call eavl1(nt,x(1,pp+2),x(1,pp+1),w,m1,m2,ll,mq,ea,ma) 1025 crm=abs(ea(1))+abs(ea(2)) 1026 j1=ma(1) 1027 j2=ma(2) 1028 4 do 7 j=1,pp 1029 if(x(mq(mp(1,j)),j).le.x(mq(mp(2,j)),j)) go to 7 1030 if(.not.(cmod)) go to 5 1031 call eav1(j,nt,x(1,j),y,w,m1,m,mp,mq,ea,ma) 1032 cr=abs(ea(1)-ym)+abs(ea(2)-ym) 1033 go to 6 1034 5 call eav1(j,nt,x(1,j),x(1,pp+1),w,m1,m,mp,mq,ea,ma) 1035 cr=abs(ea(1))+abs(ea(2)) 1036 6 if(cr .lt. crm) go to 7 1037 crm=cr 1038 js=j 1039 j1=ma(1) 1040 j2=ma(2) 1041 7 continue 1042 if(crm .gt. 0.0) go to 8 1043 term=.true. 1044 go to 9 1045 8 if(j1 .le. nu-j2) go to 9 1046 j=j1 1047 j1=nu-j2 1048 j2=nu-j 1049 9 if(term) m1=0 1050 return 1051 entry sethtf(arg) 1052 trf=arg 1053 return 1054 entry stthl(arg) 1055 thl=arg 1056 return 1057 entry sethyn(irg) 1058 ntn=max0(3,irg) 1059 call sthyn(ntn) 1060 call sthyn1(ntn) 1061 return 1062 entry stesq1(drg) 1063 esqm=drg 1064 call stesq2(esqm) 1065 return 1066 entry stlin5(irg) 1067 lin=irg 1068 call stlin6(lin) 1069 call stlin7(lin) 1070 call stlin8(lin) 1071 return 1072 entry stkn1(irg) 1073 kn=irg 1074 return 1075 end 1076 subroutine seldel (xp,n,p,x,nu,m1,js,j1,j2,a,nt,m,mp,mq,ll,yp) 1077 integer p,pp,m(2,nt,p+1),mp(2,p+2),mq(nt),ll(*) 1078 real xp(p+1),x(n,p+3) 1079 double precision a(p+1) 1080 data lin,kn /2,0/ 1081 pp=p 1082 if(kn.gt.0) pp=p+1 1083 jr=0 1084 if(m1.eq.0) go to 7 1085 if(js .le. 0) go to 3 1086 call search(2,js,nt,m,mp,j1,m1,ll) 1087 call search(1,js,nt,m,mp,j2,m2,ll(m1+1)) 1088 if(xp(js) .ge. 0.5*(x(mq(ll(m1)),js)+x(mq(ll(m1+m2)),js))) go to 7 1089 do 1 i=1,m2 1090 ll(i)=ll(i+m1) 1091 1 continue 1092 m1=m2 1093 go to 7 1094 3 jr=pp+2 1095 nh=nu-j2 1096 sx=a(p+1) 1097 do 4 j=1,p 1098 sx=sx+a(j)*xp(j) 1099 4 continue 1100 if(sx .lt. 0.5*(x(ll(j1),jr)+x(ll(nh+1),jr))) go to 5 1101 m1=j1 1102 go to 7 1103 5 do 6 i=1,j2 1104 ll(i)=ll(i+nh) 1105 6 continue 1106 m1=j2 1107 7 if((lin .ne. 1) .and. (lin .ne. 3)) go to 10 1108 if(jr .ne. pp+2) go to 8 1109 yp=sx 1110 go to 10 1111 8 yp=a(p+1) 1112 do 9 j=1,p 1113 yp=yp+a(j)*xp(j) 1114 9 continue 1115 10 return 1116 entry stlin6(irg) 1117 lin=irg 1118 return 1119 entry stkn2(irg) 1120 kn=irg 1121 return 1122 end 1123 subroutine cmprs (n,p,m,mp,nt,mo,mpo,mq,l1,l2) 1124 integer p,m(2,n,p),mp(2,p+1),mo(2,nt,p),mpo(2,p+1),mq(nt),l1(n),l2 1125 1(nt) 1126 call search(2,1,n,m,mp,nt,mt,mq) 1127 do 1 i=1,nt 1128 l1(mq(i))=i 1129 1 continue 1130 mpo(2,1)=1 1131 mpo(1,1)=nt 1132 mo(1,1,1)=0 1133 mo(2,nt,1)=mo(1,1,1) 1134 mo(2,1,1)=2 1135 mo(1,nt,1)=nt-1 1136 do 2 i=2,nt-1 1137 mo(1,i,1)=i-1 1138 mo(2,i,1)=i+1 1139 2 continue 1140 do 4 j=2,p 1141 call search(2,j,n,m,mp,nt,mt,l2) 1142 mpo(2,j)=l1(l2(1)) 1143 mpo(1,j)=l1(l2(nt)) 1144 mo(1,l1(l2(1)),j)=0 1145 mo(2,l1(l2(nt)),j)=mo(1,l1(l2(1)),j) 1146 mo(2,l1(l2(1)),j)=l1(l2(2)) 1147 mo(1,l1(l2(nt)),j)=l1(l2(nt-1)) 1148 do 3 i=2,nt-1 1149 mo(1,l1(l2(i)),j)=l1(l2(i-1)) 1150 mo(2,l1(l2(i)),j)=l1(l2(i+1)) 1151 3 continue 1152 4 continue 1153 mpo(1,p+1)=nt 1154 return 1155 end 1156 subroutine trmnl (n,p,x,y,w,m,mp,dp,ii,ir,id,itrm,rtrm,dtrm,ll) 1157 parameter(maxdim=100) 1158 integer p,itrm(*),m(2,n,p),mp(2,p+1),ll(*) 1159 real x(n,p+2),y(n),w(n),rtrm(*) 1160 double precision dp(*),dtrm(*),a(maxdim+1) 1161 data lin,kn /2,0/ 1162 itrm(ii)=ir 1163 itrm(ii+1)=id 1164 ii=ii+2 1165 nt=mp(1,p+1) 1166 i1=ii+2*p+8 1167 i2=i1+2*nt*p 1168 i3=ii+2*p+2 1169 i4=i2+nt 1170 call cmprs(n,p,m,mp,nt,itrm(i1),itrm(ii),itrm(i2),ll,ll(n+1)) 1171 jd=3 1172 if(lin.gt.0) jd=jd+(p+2)*(p+1) 1173 itrm(i3+5)=jd 1174 do 1 i=1,jd 1175 dtrm(id)=dp(i) 1176 id=id+1 1177 1 continue 1178 kn0=kn 1179 call stkn1(0) 1180 call cut(n,p,x,y,w,dp,nt,itrm(i1),itrm(ii),itrm(i2),itrm(i4), itr 1181 1m(i3),itrm(i3+1),itrm(i3+2),itrm(i3+3),itrm(i3+4),a,rtrm(ir)) 1182 call stkn1(kn0) 1183 ir=ir+1 1184 ii=ii+2*p*(nt+1)+nt+8 1185 if(itrm(i3+2) .gt. 0) go to 3 1186 ii=ii+nt 1187 do 2 i=1,nt 1188 rtrm(ir)=x(i,p+2) 1189 ir=ir+1 1190 2 continue 1191 3 if(lin .le. 0) go to 5 1192 do 4 i=1,p+1 1193 dtrm(id)=a(i) 1194 id=id+1 1195 4 continue 1196 5 return 1197 entry stlin7(irg) 1198 lin=irg 1199 return 1200 entry stkn3(irg) 1201 kn=irg 1202 return 1203 end 1204 subroutine bfnd (n,ix,x,mq,xq,lx) 1205 integer mq(n) 1206 real xq(n) 1207 ll=0 1208 lh=n+1 1209 lx=0 1210 1 if(ll+1.ge.lh) go to 12 1211 l=(ll+lh)/2 1212 if(mq(l) .ne. ix) go to 2 1213 lx=l 1214 return 1215 2 if(x .ne. xq(mq(l))) go to 9 1216 k=l+1 1217 3 if(k.gt.n) go to 5 1218 if(xq(mq(k)).gt.x) go to 5 1219 if(mq(k) .ne. ix) go to 4 1220 lx=k 1221 return 1222 4 k=k+1 1223 go to 3 1224 5 k=l-1 1225 6 if(k.lt.1) go to 8 1226 if(xq(mq(k)).lt.x) go to 8 1227 if(mq(k) .ne. ix) go to 7 1228 lx=k 1229 return 1230 7 k=k-1 1231 go to 6 1232 8 return 1233 9 if(x .ge. xq(mq(l))) go to 10 1234 lh=l 1235 go to 1 1236 10 ll=l 1237 go to 1 1238 12 return 1239 end 1240 subroutine initl (n,p,x,y,w,m,mp,dp,nd,ll) 1241 parameter (maxdim=100) 1242 integer p,m(2,n,p),mp(2,p+1),ll(n) 1243 real x(n,p),y(n),w(n) 1244 double precision dp(*),dq((maxdim+2)*(maxdim+1)+3) 1245 save dq 1246 data nd0 /0/ 1247 if(p.gt.maxdim) call diag(7,p,maxdim) 1248 if(nd0 .ne. 0) go to 1 1249 call llist(n,p,x,w,m,ll,mp) 1250 call mdlbgn(n,p,x,y,w,nd0,dq) 1251 1 nd=nd0 1252 do 2 i=1,nd0 1253 dp(i)=dq(i) 1254 2 continue 1255 return 1256 entry new(idm) 1257 nd0=0 1258 return 1259 end 1260 subroutine knn (xp,n,p,x,nt,m,mp,mq,ll) 1261 parameter (maxdim=100) 1262 integer p,pp1,pp2,m(2,nt,p+1),mp(2,p+2),mq(nt),ll(nt) 1263 real xp(p+1),x(n,p+3),sc(maxdim) 1264 data scl /0.25/ 1265 pp1=p+1 1266 pp2=p+2 1267 xp(pp1)=0.0 1268 m1=amax1(0.6,scl*nt)+0.5 1269 do 2 j=1,p 1270 if(scl .gt. 0.0) go to 1 1271 sc(j)=1.0 1272 go to 2 1273 1 call search(1,j,nt,m,mp,m1,j1,ll) 1274 sc(j)=x(mq(ll(j1)),j) 1275 call search(2,j,nt,m,mp,m1,j1,ll) 1276 sc(j)=sc(j)-x(mq(ll(j1)),j) 1277 if(sc(j).le.0.0) sc(j)=1.0 1278 sc(j)=1.0/sc(j) 1279 2 continue 1280 do 4 i=1,nt 1281 k=mq(i) 1282 s=0.0 1283 do 3 j=1,p 1284 s=s+((xp(j)-x(k,j))*sc(j))**2 1285 3 continue 1286 x(i,pp2)=s 1287 x(i,pp1)=1.0 1288 4 continue 1289 call llist(nt,1,x(1,pp2),x(1,pp1),m(1,1,pp1),ll,mp(1,pp1)) 1290 do 5 i=1,nt 1291 x(mq(i),pp1)=x(i,pp2) 1292 5 continue 1293 return 1294 entry setscl(arg) 1295 scl=arg 1296 return 1297 end 1298 subroutine xtxty (n,p,x,y,w,a,b) 1299 integer p,pp1 1300 real x(n,p),y(n),w(n) 1301 double precision a(p+1,*),b(*),s 1302 data rdg /0.01/ 1303 pp1=p+1 1304 do 5 j=1,p 1305 do 2 i=1,j 1306 s=0.d0 1307 do 1 l=1,n 1308 s=s+w(l)*x(l,i)*x(l,j) 1309 1 continue 1310 a(i,j)=s 1311 2 continue 1312 s=0.d0 1313 do 3 i=1,n 1314 s=s+w(i)*x(i,j) 1315 3 continue 1316 a(j,pp1)=s 1317 s=0.d0 1318 do 4 i=1,n 1319 s=s+w(i)*y(i)*x(i,j) 1320 4 continue 1321 b(j)=s 1322 5 continue 1323 s=0.d0 1324 do 6 i=1,n 1325 s=s+w(i)*y(i) 1326 6 continue 1327 b(pp1)=s 1328 s=0.d0 1329 do 7 i=1,n 1330 s=s+w(i) 1331 7 continue 1332 a(pp1,pp1)=s 1333 if(rdg.le.0.0) return 1334 s=float(pp1)*rdg**2/n 1335 do 8 j=1,p 1336 a(j,j)=a(j,j)+s*max(0.d0,a(j,j)-a(j,pp1)**2/a(pp1,pp1)) 1337 8 continue 1338 return 1339 entry setrdg(arg) 1340 rdg=arg 1341 return 1342 end 1343 subroutine inv (pp1,a,ai,ie) 1344 integer pp1 1345 double precision a(pp1,pp1),ai(pp1,pp1) 1346 ie=1 1347 call spofa(a,pp1,pp1,ie) 1348 if(ie.ne.0) return 1349 do 2 j=1,pp1 1350 do 1 i=1,pp1 1351 ai(i,j)=0.d0 1352 1 continue 1353 ai(j,j)=1.d0 1354 call sposl(a,pp1,pp1,ai(1,j)) 1355 2 continue 1356 return 1357 end 1358 subroutine dd (pp1,ai,b,r,y,d,ie) 1359 integer pp1 1360 double precision ai(pp1,pp1),b(pp1),r(pp1),y,d(pp1),s,t 1361 ie=0 1362 t=0.d0 1363 do 2 j=1,pp1 1364 s=0.d0 1365 do 1 i=1,pp1 1366 s=s+ai(i,j)*r(i) 1367 1 continue 1368 d(j)=s 1369 t=t+s*r(j) 1370 2 continue 1371 t=1.d0-t 1372 if(t .ne. 0.d0) go to 3 1373 ie=1 1374 return 1375 3 do 5 j=1,pp1 1376 b(j)=b(j)-y*r(j) 1377 do 4 i=1,j 1378 ai(i,j)=ai(i,j)+d(i)*d(j)/t 1379 ai(j,i)=ai(i,j) 1380 4 continue 1381 5 continue 1382 return 1383 end 1384 subroutine soln (pp1,ai,b,a) 1385 integer pp1 1386 double precision ai(pp1,pp1),b(pp1),a(pp1),s 1387 do 2 j=1,pp1 1388 s=0.d0 1389 do 1 i=1,pp1 1390 s=s+ai(i,j)*b(i) 1391 1 continue 1392 a(j)=s 1393 2 continue 1394 return 1395 end 1396 subroutine diag (ifl,i,mx) 1397 if((ifl .ge. 1) .and. (ifl .le. 16)) go to 1 1398 write(6,'('' SUBROUTINE DIAG: ifl = ,''i10,'' outside [1,16].'')' 1399 1) ifl 1400 stop 1401 1 go to (2,3,4,5,6,7,8,9,10,11,12,13,14,15, 16,17),ifl 1402 2 write(6,'('' User (input) dimension mxm ='',i10,'' is too small.'' 1403 1)') mx 1404 stop 1405 3 write(6,'('' SUBROUTINE DART1 or START1: Internal parameter mxisk 1406 1= '',i10, '' is too small.'')') mx 1407 go to 18 1408 4 write(6,'('' SUBROUTINE DART1 or START1: Internal parameter mxko = 1409 1 '',i10, '' is too small.'')') mx 1410 go to 18 1411 5 write(6,'('' User (input) dimension mxt ='',i10,'' is too small.'' 1412 1)') mx 1413 stop 1414 6 write(6,'('' SUBROUTINE DART1 or START1: Internal parameter maxid 1415 1= '',i10, '' is too small.'')') mx 1416 go to 18 1417 7 write(6,'('' Sum of observation weights not positive.'')') 1418 stop 1419 8 write(6,'('' SUBs XVALID,HYESS1,SPLIT,MDLBGN,MDLDEL,TRMNL,INITL,KN 1420 1N:'',/, ''Internal parameter maxdim = '',i5,'' MUST BE SET T 1421 1O '',I5)') mx,i 1422 go to 18 1423 9 write(6,'('' Singular input variable covariance matrix.'',/, 1424 1 '' Try increasing value of ridge (weight decay) para 1425 1meter rdg.'')') 1426 stop 1427 10 write(6,'('' Numerical error in updating covariance matrix.'',/, 1428 1 '' Try increasing value of ridge (weight decay) para 1429 1meter rdg.'')') 1430 stop 1431 11 write(6,'('' Impossible.'')') 1432 stop 1433 12 write(6,'('' User (input) dimension mxr ='',i10,'' is too small.'' 1434 1)') mx 1435 stop 1436 13 write(6,'('' User (input) dimension mxi ='',i10,'' is too small.'' 1437 1)') mx 1438 stop 1439 14 write(6,'('' User (input) dimension mxd ='',i10,'' is too small.'' 1440 1)') mx 1441 stop 1442 15 write(6,'('' SUBROUTINE HYESS or XVALID: Internal parameter mxit = 1443 1 '',i10, '' is too small.'')') mx 1444 go to 18 1445 16 write(6,'('' User (input) parameter kx ='',i10, 1446 1 '' must be less than'',i10)') i,mx 1447 stop 1448 17 write(6,'('' SUBROUTINE HYESS1: internal error - should not happen 1449 1.'')') 1450 stop 1451 18 write(6, '('' Increase its value in the parameter statement(s) an 1452 1d recompile.'')') 1453 stop 1454 end 1455 subroutine spofa(a,m,n,info) double precision a(m,*),s,t,u j1 = info do 30 j = j1, n info=j s = 0.0d0 jm1 = j - 1 if (jm1 .lt. 1) go to 20 do 10 k = 1, jm1 u=0.0 km1=k-1 if(km1.le.0) go to 40 do 50 i=1,km1 u=u+a(i,k)*a(i,j) 50 continue 40 continue t = a(k,j) - u t = t/a(k,k) a(k,j) = t s = s + t*t 10 continue 20 continue s = a(j,j) - s if (s .le. 0.0d0) return a(j,j) = dsqrt(s) 30 continue info=0 return end subroutine sposl(a,m,n,b) double precision a(m,*),b(*),t do 10 k = 1, n t = 0.0 km1=k-1 if(km1.le.0) go to 30 do 40 i=1,km1 t=t+a(i,k)*b(i) 40 continue 30 continue b(k) = (b(k) - t)/a(k,k) 10 continue do 20 kb=1,n k=n+1-kb b(k)=b(k)/a(k,k) t=-b(k) km1=k-1 if(km1.le.0) go to 50 if(t.eq.0.0) go to 50 do 60 i=1,km1 b(i)=b(i)+t*a(i,k) 60 continue 50 continue 20 continue return end subroutine psort (v,a,ii,jj) c c puts into a the permutation vector which sorts v into c increasing order. the array v is not modified. c only elements from ii to jj are considered. c arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements c c this is a modification of cacm algorithm #347 by r. c. singleton, c which is a modified hoare quicksort. c dimension a(jj),v(jj),iu(20),il(20) integer t,tt integer a real v m=1 i=ii j=jj 10 if (i.ge.j) go to 80 20 k=i ij=(j+i)/2 t=a(ij) vt=v(t) if (v(a(i)).le.vt) go to 30 a(ij)=a(i) a(i)=t t=a(ij) vt=v(t) 30 l=j if (v(a(j)).ge.vt) go to 50 a(ij)=a(j) a(j)=t t=a(ij) vt=v(t) if (v(a(i)).le.vt) go to 50 a(ij)=a(i) a(i)=t t=a(ij) vt=v(t) go to 50 40 a(l)=a(k) a(k)=tt 50 l=l-1 if (v(a(l)).gt.vt) go to 50 tt=a(l) vtt=v(tt) 60 k=k+1 if (v(a(k)).lt.vt) go to 60 if (k.le.l) go to 40 if (l-i.le.j-k) go to 70 il(m)=i iu(m)=l i=k m=m+1 go to 90 70 il(m)=k iu(m)=j j=l m=m+1 go to 90 80 m=m-1 if (m.eq.0) return i=il(m) j=iu(m) 90 if (j-i.gt.10) go to 20 if (i.eq.ii) go to 10 i=i-1 100 i=i+1 if (i.eq.j) go to 80 t=a(i+1) vt=v(t) if (v(a(i)).le.vt) go to 100 k=i 110 a(k+1)=a(k) k=k-1 if (vt.lt.v(a(k))) go to 110 a(k+1)=t go to 100 end