001:       SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            J, JOB
010:       DOUBLE PRECISION   C, GAMMA, S, SEST, SESTPR
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   W( J ), X( J )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DLAIC1 applies one step of incremental condition estimation in
020: *  its simplest version:
021: *
022: *  Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
023: *  lower triangular matrix L, such that
024: *           twonorm(L*x) = sest
025: *  Then DLAIC1 computes sestpr, s, c such that
026: *  the vector
027: *                  [ s*x ]
028: *           xhat = [  c  ]
029: *  is an approximate singular vector of
030: *                  [ L     0  ]
031: *           Lhat = [ w' gamma ]
032: *  in the sense that
033: *           twonorm(Lhat*xhat) = sestpr.
034: *
035: *  Depending on JOB, an estimate for the largest or smallest singular
036: *  value is computed.
037: *
038: *  Note that [s c]' and sestpr**2 is an eigenpair of the system
039: *
040: *      diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]
041: *                                            [ gamma ]
042: *
043: *  where  alpha =  x'*w.
044: *
045: *  Arguments
046: *  =========
047: *
048: *  JOB     (input) INTEGER
049: *          = 1: an estimate for the largest singular value is computed.
050: *          = 2: an estimate for the smallest singular value is computed.
051: *
052: *  J       (input) INTEGER
053: *          Length of X and W
054: *
055: *  X       (input) DOUBLE PRECISION array, dimension (J)
056: *          The j-vector x.
057: *
058: *  SEST    (input) DOUBLE PRECISION
059: *          Estimated singular value of j by j matrix L
060: *
061: *  W       (input) DOUBLE PRECISION array, dimension (J)
062: *          The j-vector w.
063: *
064: *  GAMMA   (input) DOUBLE PRECISION
065: *          The diagonal element gamma.
066: *
067: *  SESTPR  (output) DOUBLE PRECISION
068: *          Estimated singular value of (j+1) by (j+1) matrix Lhat.
069: *
070: *  S       (output) DOUBLE PRECISION
071: *          Sine needed in forming xhat.
072: *
073: *  C       (output) DOUBLE PRECISION
074: *          Cosine needed in forming xhat.
075: *
076: *  =====================================================================
077: *
078: *     .. Parameters ..
079:       DOUBLE PRECISION   ZERO, ONE, TWO
080:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
081:       DOUBLE PRECISION   HALF, FOUR
082:       PARAMETER          ( HALF = 0.5D0, FOUR = 4.0D0 )
083: *     ..
084: *     .. Local Scalars ..
085:       DOUBLE PRECISION   ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS,
086:      $                   NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2
087: *     ..
088: *     .. Intrinsic Functions ..
089:       INTRINSIC          ABS, MAX, SIGN, SQRT
090: *     ..
091: *     .. External Functions ..
092:       DOUBLE PRECISION   DDOT, DLAMCH
093:       EXTERNAL           DDOT, DLAMCH
094: *     ..
095: *     .. Executable Statements ..
096: *
097:       EPS = DLAMCH( 'Epsilon' )
098:       ALPHA = DDOT( J, X, 1, W, 1 )
099: *
100:       ABSALP = ABS( ALPHA )
101:       ABSGAM = ABS( GAMMA )
102:       ABSEST = ABS( SEST )
103: *
104:       IF( JOB.EQ.1 ) THEN
105: *
106: *        Estimating largest singular value
107: *
108: *        special cases
109: *
110:          IF( SEST.EQ.ZERO ) THEN
111:             S1 = MAX( ABSGAM, ABSALP )
112:             IF( S1.EQ.ZERO ) THEN
113:                S = ZERO
114:                C = ONE
115:                SESTPR = ZERO
116:             ELSE
117:                S = ALPHA / S1
118:                C = GAMMA / S1
119:                TMP = SQRT( S*S+C*C )
120:                S = S / TMP
121:                C = C / TMP
122:                SESTPR = S1*TMP
123:             END IF
124:             RETURN
125:          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
126:             S = ONE
127:             C = ZERO
128:             TMP = MAX( ABSEST, ABSALP )
129:             S1 = ABSEST / TMP
130:             S2 = ABSALP / TMP
131:             SESTPR = TMP*SQRT( S1*S1+S2*S2 )
132:             RETURN
133:          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
134:             S1 = ABSGAM
135:             S2 = ABSEST
136:             IF( S1.LE.S2 ) THEN
137:                S = ONE
138:                C = ZERO
139:                SESTPR = S2
140:             ELSE
141:                S = ZERO
142:                C = ONE
143:                SESTPR = S1
144:             END IF
145:             RETURN
146:          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
147:             S1 = ABSGAM
148:             S2 = ABSALP
149:             IF( S1.LE.S2 ) THEN
150:                TMP = S1 / S2
151:                S = SQRT( ONE+TMP*TMP )
152:                SESTPR = S2*S
153:                C = ( GAMMA / S2 ) / S
154:                S = SIGN( ONE, ALPHA ) / S
155:             ELSE
156:                TMP = S2 / S1
157:                C = SQRT( ONE+TMP*TMP )
158:                SESTPR = S1*C
159:                S = ( ALPHA / S1 ) / C
160:                C = SIGN( ONE, GAMMA ) / C
161:             END IF
162:             RETURN
163:          ELSE
164: *
165: *           normal case
166: *
167:             ZETA1 = ALPHA / ABSEST
168:             ZETA2 = GAMMA / ABSEST
169: *
170:             B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
171:             C = ZETA1*ZETA1
172:             IF( B.GT.ZERO ) THEN
173:                T = C / ( B+SQRT( B*B+C ) )
174:             ELSE
175:                T = SQRT( B*B+C ) - B
176:             END IF
177: *
178:             SINE = -ZETA1 / T
179:             COSINE = -ZETA2 / ( ONE+T )
180:             TMP = SQRT( SINE*SINE+COSINE*COSINE )
181:             S = SINE / TMP
182:             C = COSINE / TMP
183:             SESTPR = SQRT( T+ONE )*ABSEST
184:             RETURN
185:          END IF
186: *
187:       ELSE IF( JOB.EQ.2 ) THEN
188: *
189: *        Estimating smallest singular value
190: *
191: *        special cases
192: *
193:          IF( SEST.EQ.ZERO ) THEN
194:             SESTPR = ZERO
195:             IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
196:                SINE = ONE
197:                COSINE = ZERO
198:             ELSE
199:                SINE = -GAMMA
200:                COSINE = ALPHA
201:             END IF
202:             S1 = MAX( ABS( SINE ), ABS( COSINE ) )
203:             S = SINE / S1
204:             C = COSINE / S1
205:             TMP = SQRT( S*S+C*C )
206:             S = S / TMP
207:             C = C / TMP
208:             RETURN
209:          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
210:             S = ZERO
211:             C = ONE
212:             SESTPR = ABSGAM
213:             RETURN
214:          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
215:             S1 = ABSGAM
216:             S2 = ABSEST
217:             IF( S1.LE.S2 ) THEN
218:                S = ZERO
219:                C = ONE
220:                SESTPR = S1
221:             ELSE
222:                S = ONE
223:                C = ZERO
224:                SESTPR = S2
225:             END IF
226:             RETURN
227:          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
228:             S1 = ABSGAM
229:             S2 = ABSALP
230:             IF( S1.LE.S2 ) THEN
231:                TMP = S1 / S2
232:                C = SQRT( ONE+TMP*TMP )
233:                SESTPR = ABSEST*( TMP / C )
234:                S = -( GAMMA / S2 ) / C
235:                C = SIGN( ONE, ALPHA ) / C
236:             ELSE
237:                TMP = S2 / S1
238:                S = SQRT( ONE+TMP*TMP )
239:                SESTPR = ABSEST / S
240:                C = ( ALPHA / S1 ) / S
241:                S = -SIGN( ONE, GAMMA ) / S
242:             END IF
243:             RETURN
244:          ELSE
245: *
246: *           normal case
247: *
248:             ZETA1 = ALPHA / ABSEST
249:             ZETA2 = GAMMA / ABSEST
250: *
251:             NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ),
252:      $              ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 )
253: *
254: *           See if root is closer to zero or to ONE
255: *
256:             TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
257:             IF( TEST.GE.ZERO ) THEN
258: *
259: *              root is close to zero, compute directly
260: *
261:                B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
262:                C = ZETA2*ZETA2
263:                T = C / ( B+SQRT( ABS( B*B-C ) ) )
264:                SINE = ZETA1 / ( ONE-T )
265:                COSINE = -ZETA2 / T
266:                SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
267:             ELSE
268: *
269: *              root is closer to ONE, shift by that amount
270: *
271:                B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
272:                C = ZETA1*ZETA1
273:                IF( B.GE.ZERO ) THEN
274:                   T = -C / ( B+SQRT( B*B+C ) )
275:                ELSE
276:                   T = B - SQRT( B*B+C )
277:                END IF
278:                SINE = -ZETA1 / T
279:                COSINE = -ZETA2 / ( ONE+T )
280:                SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
281:             END IF
282:             TMP = SQRT( SINE*SINE+COSINE*COSINE )
283:             S = SINE / TMP
284:             C = COSINE / TMP
285:             RETURN
286: *
287:          END IF
288:       END IF
289:       RETURN
290: *
291: *     End of DLAIC1
292: *
293:       END
294: