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