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