The R Project SVN R-packages

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3821 mrmanese 1 
 
2 
 
3 
      SUBROUTINE INCLUD(NP, NRBAR, WEIGHT, XROW, YELEM, D,
4 
     +      RBAR, THETAB, SSERR, IER)
5 
C
6 
C     ALGORITHM AS274  APPL. STATIST. (1992) VOL.41, NO. 2
7 
C     Modified from algorithm AS 75.1
8 
C
9 
C     Calling this routine updates d, rbar, thetab and sserr by the
10 
C     inclusion of xrow, yelem with the specified weight.   The number
11 
C     of columns (variables) may exceed the number of rows (cases).
12 
C
13 
C**** WARNING: The elements of XROW are overwritten  ****
14 
C
15 
      INTEGER NP, NRBAR, IER
16 
      DOUBLE PRECISION WEIGHT, XROW(NP), YELEM, D(NP), RBAR(*),
17 
     +    THETAB(NP), SSERR
18 
C
19 
C     Local variables
20 
C
21 
      INTEGER I, K, NEXTR
22 
      DOUBLE PRECISION ZERO, W, Y, XI, DI, WXI, DPI, CBAR, SBAR, XK
23 
C
24 
      DATA ZERO/0.D0/
25 
C
26 
C     Some checks.
27 
C
28 
      IER = 0
29 
      IF (NP .LT. 1) IER = 1
30 
      IF (NRBAR .LT. NP*(NP-1)/2) IER = IER + 2
31 
      IF (IER .NE. 0) RETURN
32 
C
33 
      W = WEIGHT
34 
      Y = YELEM
35 
      NEXTR = 1
36 
      DO 30 I = 1, NP
37 
C
38 
C     Skip unnecessary transformations.   Test on exact zeroes must be
39 
C     used or stability can be destroyed.
40 
C
41 
      IF (W .EQ. ZERO) RETURN
42 
       XI = XROW(I)
43 
      IF (XI .EQ. ZERO) THEN
44 
       NEXTR = NEXTR + NP - I
45 
      GO TO 30
46 
      END IF
47 
      DI = D(I)
48 
      WXI = W * XI
49 
      DPI = DI + WXI*XI
50 
      CBAR = DI / DPI
51 
      SBAR = WXI / DPI
52 
      W = CBAR * W
53 
      D(I) = DPI
54 
      IF (I .EQ. NP) GO TO 20
55 
      DO 10 K = I+1, NP
56 
	XK = XROW(K)
57 
	XROW(K) = XK - XI * RBAR(NEXTR)
58 
	RBAR(NEXTR) = CBAR * RBAR(NEXTR) + SBAR * XK
59 
	NEXTR = NEXTR + 1
60 
   10   CONTINUE
61 
   20   XK = Y
62 
       Y = XK - XI * THETAB(I)
63 
       THETAB(I) = CBAR * THETAB(I) + SBAR * XK
64 
   30  CONTINUE
65 
C
66 
C     Y * SQRT(W) is now equal to Brown & Durbin's recursive residual.
67 
C
68 
      SSERR = SSERR + W * Y * Y
69 
C
70 
      RETURN
71 
      END
72 
C
73 
C
74 
      SUBROUTINE TOLSET(NP, NRBAR, D, RBAR, TOL, WORK, IER)
75 
C
76 
C     ALGORITHM AS274  APPL. STATIST. (1992) VOL.41, NO. 2
77 
C
78 
C     Sets up array TOL for testing for zeroes in an orthogonal
79 
C     reduction formed using AS75.1.
80 
C
81 
      INTEGER NP, NRBAR, IER
82 
      DOUBLE PRECISION D(NP), RBAR(*), TOL(NP), WORK(NP)
83 
C
84 
C     Local variables.
85 
C
86 
      INTEGER COL, ROW, POS
87 
      DOUBLE PRECISION EPS, SUM, ZERO
88 
C
89 
C     EPS is a machine-dependent constant.   For compilers which use
90 
C     the IEEE format for floating-point numbers, recommended values
91 
C     are 1.E-06 for single precision and 1.D-12 for double precision.
92 
C
93 
c     changed EPS from 10^-12 to 5x10^-10 to try to fix a bug
94 
      DATA EPS/1.D-12/, ZERO/0.D0/
95 
C
96 
C     Some checks.
97 
C
98 
      IER = 0
99 
      IF (NP .LT. 1) IER = 1
100 
      IF (NRBAR .LT. NP*(NP-1)/2) IER = IER + 2
101 
      IF (IER .NE. 0) RETURN
102 
C
103 
C     Set TOL(I) = sum of absolute values in column I of RBAR after
104 
C     scaling each element by the square root of its row multiplier.
105 
C
106 
      DO 10 ROW = 1, NP
107 
   10 WORK(ROW) = SQRT(D(ROW))
108 
      DO 30 COL = 1, NP
109 
      POS = COL - 1
110 
      IF (COL .LE. NP) THEN
111 
      SUM = WORK(COL)
112 
      ELSE
113 
      SUM = ZERO
114 
      END IF
115 
      DO 20 ROW = 1, MIN(COL-1, NP)
116 
      SUM = SUM + ABS(RBAR(POS)) * WORK(ROW)
117 
      POS = POS + NP - ROW - 1
118 
  20  CONTINUE
119 
      TOL(COL) = EPS * SUM
120 
  30  CONTINUE
121 
C
122 
      RETURN
123 
      END
124 
 
125 
      SUBROUTINE SINGCHK(NP, NRBAR, D, RBAR, THETAB, SSERR, TOL,
126 
     +   LINDEP, WORK, IER)
127 
C
128 
C     ALGORITHM AS274  APPL. STATIST. (1992) VOL.41, NO. 2
129 
C
130 
C     Checks for singularities, reports, and adjusts orthogonal
131 
C     reductions produced by AS75.1.
132 
C
133 
      INTEGER NP, NRBAR, IER
134 
      DOUBLE PRECISION D(NP), RBAR(NRBAR), THETAB(NP), SSERR,
135 
     +      TOL(NP), WORK(NP)
136 
      LOGICAL LINDEP(NP)
137 
C
138 
C     Local variables
139 
C
140 
      DOUBLE PRECISION ZERO, TEMP
141 
      INTEGER COL, POS, ROW, NC2, POS2
142 
C
143 
      DATA ZERO/0.D0/
144 
C
145 
C     Check input parameters
146 
C
147 
      IER = 0
148 
      IF (NP .LT. 1) IER = 1
149 
      IF (NRBAR .LT. NP*(NP-1)/2) IER = IER + 2
150 
      IF (IER .NE. 0) RETURN
151 
C
152 
      DO 10 COL = 1, NP
153 
   10 WORK(COL) = SQRT(D(COL))
154 
C
155 
      DO 40 COL = 1, NP
156 
C
157 
C     Set elements within RBAR to zero if they are less than TOL(COL) in
158 
C     absolute value after being scaled by the square root of their row
159 
C     multiplier.
160 
C
161 
      TEMP = TOL(COL)
162 
      POS = COL - 1
163 
      DO 30 ROW = 1, COL-1
164 
      IF (ABS(RBAR(POS)) * WORK(ROW) .LT. TEMP) RBAR(POS) = ZERO
165 
      POS = POS + NP - ROW - 1
166 
   30 CONTINUE
167 
C
168 
C     If diagonal element is near zero, set it to zero, set appropriate
169 
C     element of LINDEP, and use INCLUD to augment the projections in
170 
C     the lower rows of the orthogonalization.
171 
C
172 
      LINDEP(COL) = .FALSE.
173 
      IF (WORK(COL) .LE. TEMP) THEN
174 
      LINDEP(COL) = .TRUE.
175 
      IER = IER - 1
176 
      IF (COL .LT. NP) THEN
177 
	NC2 = NP - COL
178 
	POS2 = POS + NP - COL + 1
179 
	CALL INCLUD(NC2, NC2*(NC2-1)/2, D(COL), RBAR(POS+1),
180 
     +            THETAB(COL), D(COL+1), RBAR(POS2), THETAB(COL+1),
181 
     +            SSERR, IER)
182 
      ELSE
183 
	SSERR = SSERR + D(COL) * THETAB(COL)**2
184 
      END IF
185 
      D(COL) = ZERO
186 
      WORK(COL) = ZERO
187 
      THETAB(COL) = ZERO
188 
      END IF
189 
   40 CONTINUE
190 
      RETURN
191 
      END
192 
 
193 
 
194 
      SUBROUTINE REGCF(NP, NRBAR, D, RBAR, THETAB, TOL, BETA,
195 
     +     NREQ, IER)
196 
C
197 
C     ALGORITHM AS274  APPL. STATIST. (1992) VOL 41, NO. x
198 
C
199 
C     Modified version of AS75.4 to calculate regression coefficients
200 
C     for the first NREQ variables, given an orthogonal reduction from
201 
C     AS75.1.
202 
C
203 
      INTEGER NP, NRBAR, NREQ, IER
204 
      DOUBLE PRECISION D(NP), RBAR(*), THETAB(NP), TOL(NP),
205 
     +     BETA(NP)
206 
C
207 
C     Local variables
208 
C
209 
      INTEGER I, J, NEXTR
210 
      DOUBLE PRECISION ZERO
211 
C
212 
      DATA ZERO/0.D0/
213 
C
214 
C     Some checks.
215 
C
216 
      IER = 0
217 
      IF (NP .LT. 1) IER = 1
218 
      IF (NRBAR .LT. NP*(NP-1)/2) IER = IER + 2
219 
      IF (NREQ .LT. 1 .OR. NREQ .GT. NP) IER = IER + 4
220 
      IF (IER .NE. 0) RETURN
221 
C
222 
      DO 20 I = NREQ, 1, -1
223 
      IF (SQRT(D(I)) .LT. TOL(I)) THEN
224 
      BETA(I) = ZERO
225 
      D(I) = ZERO
226 
      GO TO 20
227 
      END IF
228 
      BETA(I) = THETAB(I)
229 
      NEXTR = (I-1) * (NP+NP-I)/2 + 1
230 
      DO 10 J = I+1, NREQ
231 
      BETA(I) = BETA(I) - RBAR(NEXTR) * BETA(J)
232 
      NEXTR = NEXTR + 1
233 
   10 CONTINUE
234 
   20 CONTINUE
235 
C
236 
      RETURN
237 
      END