1c--- For lmrob.lar() in ../R/lmrob.M.S.R 2c--- ~~~~~~~~~~~ 3C======================================================================= 4 SUBROUTINE rlSTORm2(Y,N,J,YJ) 5C....................................................................... 6 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 7 DIMENSION Y(N) 8C----------------------------------------------------------------------- 9C rlSTORm2 SEARCHES THE J-TH VALUE IN ORDER OF MAGNITUDE IN 10C A VECTOR OF LENGTH N. 11C----------------------------------------------------------------------- 12C--- copied from robust package: src/lmrobmm.f ------------------------- 13 L=1 14 LR=N 15 20 IF (L.GE.LR) GOTO 90 16 AX=Y(J) 17 JNC=L 18 JJ=LR 19 30 IF(JNC.GT.JJ) GOTO 80 20 40 IF (Y(JNC).GE.AX) GOTO 50 21 JNC=JNC+1 22 GOTO 40 23 50 IF(Y(JJ).LE.AX) GOTO 60 24 JJ=JJ-1 25 GOTO 50 26 60 IF(JNC.GT.JJ) GOTO 70 27 WA=Y(JNC) 28 Y(JNC)=Y(JJ) 29 Y(JJ)=WA 30 JNC=JNC+1 31 JJ=JJ-1 32 70 GOTO 30 33 80 IF(JJ.LT.J) L=JNC 34 IF(J.LT.JNC) LR=JJ 35 GOTO 20 36 90 YJ=Y(J) 37 RETURN 38 END 39C======================================================================= 40 SUBROUTINE rlCOLbi(V1,V2,MLT,M,IOUT) 41C....................................................................... 42 DOUBLE PRECISION V1(M),V2(M),MLT 43C----------------------------------------------------------------------- 44C AUXILIARY ROUTINE FOR rlLARSbi 45C----------------------------------------------------------------------- 46C--- copied from robust package: src/lmrobbi.f ------------------------- 47 DO 220 I=1,M 48 IF (I .EQ. IOUT) GOTO 220 49 V1(I)=V1(I)-V2(I)*MLT 50 220 CONTINUE 51 RETURN 52 END 53C======================================================================= 54 SUBROUTINE rlICHGbi(A,B) 55C....................................................................... 56C AUXILIARY ROUTINE FOR rlLARSbi 57C----------------------------------------------------------------------- 58C--- copied from robust package: src/lmrobbi.f ------------------------- 59 DOUBLE PRECISION A,B,C 60 C=A 61 A=B 62 B=C 63 RETURN 64 END 65C======================================================================= 66 SUBROUTINE rlLARSbi(X,Y,N,NP,MDX,MDT,TOL,NIT,K, 67 + KODE,SIGMA,THETA,RS,SC1,SC2,SC3,SC4,BET0) 68C....................................................................... 69 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 70 DIMENSION X(MDX,NP),Y(N),THETA(MDT),RS(N),SC1(N),SC2(NP), 71 + SC3(NP),SC4(NP) 72 INTEGER OUT 73 LOGICAL STAGE,TEST 74 DATA ZERO,TWO,EPS,BIG/0.D0,2.D0,1.0D-10,3.401D38/ 75cMM would think rather this: double precision --- but it breaks our checks 76C DATA ZERO,TWO,EPS,BIG/0.D0,2.D0,2.22D-16,1.796D308/ 77 78C----------------------------------------------------------------------- 79C LEAST ABSOLUTE RESIDUALS -- aka L_1 - Regression 80C --> Result in THETA[1:NP] 81C----------------------------------------------------------------------- 82C--- copied from robust package: src/lmrobbi.f ------------------------- 83 DO J=1,NP 84 SC4(J)=DBLE(J) 85 SC2(J)=ZERO 86 end do 87 SUM=ZERO 88 DO I=1,N 89 SC1(I)=DBLE(NP+I) 90 THETA(I)=Y(I) 91 IF (Y(I) .lt. ZERO) then 92 DO J=1,NP 93 X(I,J)=-X(I,J) 94 end do 95 THETA(I)=-THETA(I) 96 SC1(I)=-SC1(I) 97 endif 98 SUM=SUM+THETA(I) 99 end do 100C----------------------------------------------------------------------- 101C COMPUTE THE MARGINAL COSTS. 102C----------------------------------------------------------------------- 103 SUMIN=SUM 104 DO J=1,NP 105 SUM=ZERO 106 DO I=1,N 107 SUM=SUM+X(I,J) 108 end do 109 SC3(J)=SUM 110 end do 111C----------------------------------------------------------------------- 112C STAGE I. DETERMINE THE VECTOR TO ENTER THE BASIS. 113C----------------------------------------------------------------------- 114 TEST=.FALSE. ! -Wall 115 STAGE=.TRUE. 116 KOUNT=0 117 KR=1 118 KL=1 119 IN=1 ! -Wall 120 121c-- ---------------- LOOP (Stage I) ------------------------------------ 122 70 VMAX=-1.D0 123 DNP=DBLE(NP) 124 DO J=KR,NP 125 IF (DABS(SC4(J)) .GT. DNP) cycle ! = continue 126 D=DABS(SC3(J)) 127 IF (D-VMAX .LE. ZERO) cycle 128 IF (D-VMAX .LE. EPS) cycle 129 VMAX=D 130 IN=J 131 end do 132 IF (SC3(IN) .lt. ZERO) then ! swap signs 133 do I=1,N 134 X(I,IN)=-X(I,IN) 135 end do 136 SC3(IN)=-SC3(IN) 137 SC4(IN)=-SC4(IN) 138 endif 139 140C----------------------------------------------------------------------- 141C DETERMINE THE VECTOR TO LEAVE THE BASIS. 142C----------------------------------------------------------------------- 143 144cvvv ------------ 2nd-level loop --------------------------------- 145 100 K=0 146 DO I=KL,N 147 D=X(I,IN) 148 IF (D .LE. TOL) cycle 149 K=K+1 150 Y(K)=THETA(I)/D 151 RS(K)=DBLE(I) 152 TEST=.TRUE. 153 end do 154 155C--- -------------- 3rd-level loop ------------------ 156 120 IF (K .le. 0) then 157 TEST=.FALSE. ! and GOTO 150 158 else ! 130 159 VMIN=BIG 160 DO I=1,K 161 IF (Y(I)-VMIN .GE. ZERO) cycle 162 IF (VMIN-Y(I) .LE. EPS) cycle 163 J=I 164 VMIN=Y(I) 165 OUT=INT(RS(I)) 166 end do 167 Y(J)=Y(K) 168 RS(J)=RS(K) 169 K=K-1 170 endif 171 172C----------------------------------------------------------------------- 173C CHECK FOR LINEAR DEPENDENCE IN STAGE I. 174C----------------------------------------------------------------------- 175c 150 176 IF (.not.TEST .and. STAGE) then 177 DO I=1,N 178 CALL rlICHGbi(X(I,KR),X(I,IN)) 179 end do 180 CALL rlICHGbi(SC3(KR),SC3(IN)) 181 CALL rlICHGbi(SC4(KR),SC4(IN)) 182 KR=KR+1 183c GOTO 260 184 else 185c 170 186 IF (.not. TEST) then 187 KODE=2 188 GOTO 350 189 endif 190c 180 191 PIVOT=X(OUT,IN) 192 IF (SC3(IN)-PIVOT-PIVOT .gt. TOL) then ! not converged 193 DO J=KR,NP 194 D=X(OUT,J) 195 SC3(J)=SC3(J)-D-D 196 X(OUT,J)=-D 197 end do 198 D=THETA(OUT) 199 SUMIN=SUMIN-D-D 200 THETA(OUT)=-D 201 SC1(OUT)=-SC1(OUT) 202 GOTO 120 203c -----------end{ 3rd-level loop } ----------------- 204 endif 205 206C----------------------------------------------------------------------- 207C 200 PIVOT ON X(OUT,IN). 208C----------------------------------------------------------------------- 209 DO J=KR,NP 210 IF (J.EQ.IN) cycle ! = continue 211 X(OUT,J)=X(OUT,J)/PIVOT 212 end do 213 THETA(OUT)=THETA(OUT)/PIVOT 214 DO J=KR,NP 215 IF (J .EQ. IN) cycle 216 D=X(OUT,J) 217 SC3(J)=SC3(J)-D*SC3(IN) 218 CALL rlCOLbi(X(1,J),X(1,IN),D,N,OUT) 219 end do 220 SUMIN=SUMIN-SC3(IN)*THETA(OUT) 221 DO I=1,N 222 IF (I .EQ. OUT) cycle 223 D=X(I,IN) 224 THETA(I)=THETA(I)-D*THETA(OUT) 225 X(I,IN)=-D/PIVOT 226 end do 227 SC3(IN)=-SC3(IN)/PIVOT 228 X(OUT,IN)=1.D0/PIVOT 229 CALL rlICHGbi(SC1(OUT),SC4(IN)) 230 KOUNT=KOUNT+1 231 IF (.NOT. STAGE) GOTO 270 232C----------------------------------------------------------------------- 233C INTERCHANGE ROWS IN STAGE I. 234C----------------------------------------------------------------------- 235 KL=KL+1 236 DO J=KR,NP 237 CALL rlICHGbi(X(OUT,J),X(KOUNT,J)) 238 enddo 239 CALL rlICHGbi(THETA(OUT),THETA(KOUNT)) 240 CALL rlICHGbi(SC1(OUT),SC1(KOUNT)) 241 endif 242 243 IF (KOUNT+KR .NE. NP+1) GOTO 70 244c ======= 245 246C----------------------------------------------------------------------- 247C STAGE II. DETERMINE THE VECTOR TO ENTER THE BASIS. 248C----------------------------------------------------------------------- 249 STAGE=.FALSE. 250cvvv 251 270 VMAX=-BIG 252 DO J=KR,NP 253 D=SC3(J) 254 IF (D .lt. ZERO) then 255 IF (D+TWO .GT. ZERO) cycle 256 D=-D-TWO 257 endif 258 IF (D-VMAX .LE. ZERO) cycle 259 IF (D-VMAX .LE. EPS) cycle 260 VMAX=D 261 IN=J 262 end do 263 IF (VMAX .gt. TOL) then ! not converged 264 IF (SC3(IN) .le. ZERO) then 265 DO I=1,N 266 X(I,IN)=-X(I,IN) 267 end do 268 SC3(IN)=-SC3(IN)-2.D0 269 SC4(IN)=-SC4(IN) 270 endif 271 GOTO 100 272c ======== 273 endif 274C----------------------------------------------------------------------- 275C 310 PREPARE OUTPUT 276C----------------------------------------------------------------------- 277 L=KL-1 278 DO I=1,N 279 RS(I)=ZERO 280 IF (I .GT. L .OR. THETA(I) .GE. ZERO) cycle 281 do J=KR,NP 282 X(I,J)=-X(I,J) 283 end do 284 THETA(I)=-THETA(I) 285 SC1(I)=-SC1(I) 286 end do 287 KODE=0 288 IF (KR .eq. 1) then ! first time only 289 do J=1,NP 290 D=DABS(SC3(J)) 291 IF (D .LE. TOL .OR. TWO-D .LE. TOL) GOTO 350 292 end do 293 KODE=1 294 endif 295c--- 296 350 DO I=1,N 297 K=INT(SC1(I)) 298 D=THETA(I) 299 IF (K .le. 0) then 300 K=-K 301 D=-D 302 endif 303 IF (I .lt. KL) then 304 SC2(K)=D 305 else 306 K=K-NP 307 RS(K)=D 308 endif 309 end do 310 K=NP+1-KR 311 SUM=ZERO 312 DO I=KL,N 313 SUM=SUM+THETA(I) 314 end do 315 SUMIN=SUM 316 NIT=KOUNT 317 DO J=1,NP 318 THETA(J)=SC2(J) 319 end do 320 DO I=1,N 321 Y(I)=DABS(RS(I)) 322 end do 323 N2=N/2+1 324 CALL RLSTORM2(Y,N,N2,SIGMA) 325 SIGMA=SIGMA/BET0 326 RETURN 327 END 328