001:       SUBROUTINE CLAIC1( 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               SEST, SESTPR
010:       COMPLEX            C, GAMMA, S
011: *     ..
012: *     .. Array Arguments ..
013:       COMPLEX            W( J ), X( J )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  CLAIC1 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 CLAIC1 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] * [ conjg(alpha) ]
041: *                                            [ conjg(gamma) ]
042: *
043: *  where  alpha =  conjg(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) COMPLEX array, dimension (J)
056: *          The j-vector x.
057: *
058: *  SEST    (input) REAL
059: *          Estimated singular value of j by j matrix L
060: *
061: *  W       (input) COMPLEX array, dimension (J)
062: *          The j-vector w.
063: *
064: *  GAMMA   (input) COMPLEX
065: *          The diagonal element gamma.
066: *
067: *  SESTPR  (output) REAL
068: *          Estimated singular value of (j+1) by (j+1) matrix Lhat.
069: *
070: *  S       (output) COMPLEX
071: *          Sine needed in forming xhat.
072: *
073: *  C       (output) COMPLEX
074: *          Cosine needed in forming xhat.
075: *
076: *  =====================================================================
077: *
078: *     .. Parameters ..
079:       REAL               ZERO, ONE, TWO
080:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
081:       REAL               HALF, FOUR
082:       PARAMETER          ( HALF = 0.5E0, FOUR = 4.0E0 )
083: *     ..
084: *     .. Local Scalars ..
085:       REAL               ABSALP, ABSEST, ABSGAM, B, EPS, NORMA, S1, S2,
086:      $                   SCL, T, TEST, TMP, ZETA1, ZETA2
087:       COMPLEX            ALPHA, COSINE, SINE
088: *     ..
089: *     .. Intrinsic Functions ..
090:       INTRINSIC          ABS, CONJG, MAX, SQRT
091: *     ..
092: *     .. External Functions ..
093:       REAL               SLAMCH
094:       COMPLEX            CDOTC
095:       EXTERNAL           SLAMCH, CDOTC
096: *     ..
097: *     .. Executable Statements ..
098: *
099:       EPS = SLAMCH( 'Epsilon' )
100:       ALPHA = CDOTC( J, X, 1, W, 1 )
101: *
102:       ABSALP = ABS( ALPHA )
103:       ABSGAM = ABS( GAMMA )
104:       ABSEST = ABS( SEST )
105: *
106:       IF( JOB.EQ.1 ) THEN
107: *
108: *        Estimating largest singular value
109: *
110: *        special cases
111: *
112:          IF( SEST.EQ.ZERO ) THEN
113:             S1 = MAX( ABSGAM, ABSALP )
114:             IF( S1.EQ.ZERO ) THEN
115:                S = ZERO
116:                C = ONE
117:                SESTPR = ZERO
118:             ELSE
119:                S = ALPHA / S1
120:                C = GAMMA / S1
121:                TMP = SQRT( S*CONJG( S )+C*CONJG( C ) )
122:                S = S / TMP
123:                C = C / TMP
124:                SESTPR = S1*TMP
125:             END IF
126:             RETURN
127:          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
128:             S = ONE
129:             C = ZERO
130:             TMP = MAX( ABSEST, ABSALP )
131:             S1 = ABSEST / TMP
132:             S2 = ABSALP / TMP
133:             SESTPR = TMP*SQRT( S1*S1+S2*S2 )
134:             RETURN
135:          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
136:             S1 = ABSGAM
137:             S2 = ABSEST
138:             IF( S1.LE.S2 ) THEN
139:                S = ONE
140:                C = ZERO
141:                SESTPR = S2
142:             ELSE
143:                S = ZERO
144:                C = ONE
145:                SESTPR = S1
146:             END IF
147:             RETURN
148:          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
149:             S1 = ABSGAM
150:             S2 = ABSALP
151:             IF( S1.LE.S2 ) THEN
152:                TMP = S1 / S2
153:                SCL = SQRT( ONE+TMP*TMP )
154:                SESTPR = S2*SCL
155:                S = ( ALPHA / S2 ) / SCL
156:                C = ( GAMMA / S2 ) / SCL
157:             ELSE
158:                TMP = S2 / S1
159:                SCL = SQRT( ONE+TMP*TMP )
160:                SESTPR = S1*SCL
161:                S = ( ALPHA / S1 ) / SCL
162:                C = ( GAMMA / S1 ) / SCL
163:             END IF
164:             RETURN
165:          ELSE
166: *
167: *           normal case
168: *
169:             ZETA1 = ABSALP / ABSEST
170:             ZETA2 = ABSGAM / ABSEST
171: *
172:             B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
173:             C = ZETA1*ZETA1
174:             IF( B.GT.ZERO ) THEN
175:                T = C / ( B+SQRT( B*B+C ) )
176:             ELSE
177:                T = SQRT( B*B+C ) - B
178:             END IF
179: *
180:             SINE = -( ALPHA / ABSEST ) / T
181:             COSINE = -( GAMMA / ABSEST ) / ( ONE+T )
182:             TMP = SQRT( SINE*CONJG( SINE )+COSINE*CONJG( COSINE ) )
183:             S = SINE / TMP
184:             C = COSINE / TMP
185:             SESTPR = SQRT( T+ONE )*ABSEST
186:             RETURN
187:          END IF
188: *
189:       ELSE IF( JOB.EQ.2 ) THEN
190: *
191: *        Estimating smallest singular value
192: *
193: *        special cases
194: *
195:          IF( SEST.EQ.ZERO ) THEN
196:             SESTPR = ZERO
197:             IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
198:                SINE = ONE
199:                COSINE = ZERO
200:             ELSE
201:                SINE = -CONJG( GAMMA )
202:                COSINE = CONJG( ALPHA )
203:             END IF
204:             S1 = MAX( ABS( SINE ), ABS( COSINE ) )
205:             S = SINE / S1
206:             C = COSINE / S1
207:             TMP = SQRT( S*CONJG( S )+C*CONJG( C ) )
208:             S = S / TMP
209:             C = C / TMP
210:             RETURN
211:          ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
212:             S = ZERO
213:             C = ONE
214:             SESTPR = ABSGAM
215:             RETURN
216:          ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
217:             S1 = ABSGAM
218:             S2 = ABSEST
219:             IF( S1.LE.S2 ) THEN
220:                S = ZERO
221:                C = ONE
222:                SESTPR = S1
223:             ELSE
224:                S = ONE
225:                C = ZERO
226:                SESTPR = S2
227:             END IF
228:             RETURN
229:          ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
230:             S1 = ABSGAM
231:             S2 = ABSALP
232:             IF( S1.LE.S2 ) THEN
233:                TMP = S1 / S2
234:                SCL = SQRT( ONE+TMP*TMP )
235:                SESTPR = ABSEST*( TMP / SCL )
236:                S = -( CONJG( GAMMA ) / S2 ) / SCL
237:                C = ( CONJG( ALPHA ) / S2 ) / SCL
238:             ELSE
239:                TMP = S2 / S1
240:                SCL = SQRT( ONE+TMP*TMP )
241:                SESTPR = ABSEST / SCL
242:                S = -( CONJG( GAMMA ) / S1 ) / SCL
243:                C = ( CONJG( ALPHA ) / S1 ) / SCL
244:             END IF
245:             RETURN
246:          ELSE
247: *
248: *           normal case
249: *
250:             ZETA1 = ABSALP / ABSEST
251:             ZETA2 = ABSGAM / ABSEST
252: *
253:             NORMA = MAX( ONE+ZETA1*ZETA1+ZETA1*ZETA2,
254:      $              ZETA1*ZETA2+ZETA2*ZETA2 )
255: *
256: *           See if root is closer to zero or to ONE
257: *
258:             TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
259:             IF( TEST.GE.ZERO ) THEN
260: *
261: *              root is close to zero, compute directly
262: *
263:                B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
264:                C = ZETA2*ZETA2
265:                T = C / ( B+SQRT( ABS( B*B-C ) ) )
266:                SINE = ( ALPHA / ABSEST ) / ( ONE-T )
267:                COSINE = -( GAMMA / ABSEST ) / T
268:                SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
269:             ELSE
270: *
271: *              root is closer to ONE, shift by that amount
272: *
273:                B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
274:                C = ZETA1*ZETA1
275:                IF( B.GE.ZERO ) THEN
276:                   T = -C / ( B+SQRT( B*B+C ) )
277:                ELSE
278:                   T = B - SQRT( B*B+C )
279:                END IF
280:                SINE = -( ALPHA / ABSEST ) / T
281:                COSINE = -( GAMMA / ABSEST ) / ( ONE+T )
282:                SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
283:             END IF
284:             TMP = SQRT( SINE*CONJG( SINE )+COSINE*CONJG( COSINE ) )
285:             S = SINE / TMP
286:             C = COSINE / TMP
287:             RETURN
288: *
289:          END IF
290:       END IF
291:       RETURN
292: *
293: *     End of CLAIC1
294: *
295:       END
296: