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