Annotation of rpl/lapack/lapack/dlaic1.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLAIC1
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLAIC1 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaic1.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaic1.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaic1.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * INTEGER J, JOB
! 25: * DOUBLE PRECISION C, GAMMA, S, SEST, SESTPR
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * DOUBLE PRECISION W( J ), X( J )
! 29: * ..
! 30: *
! 31: *
! 32: *> \par Purpose:
! 33: * =============
! 34: *>
! 35: *> \verbatim
! 36: *>
! 37: *> DLAIC1 applies one step of incremental condition estimation in
! 38: *> its simplest version:
! 39: *>
! 40: *> Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
! 41: *> lower triangular matrix L, such that
! 42: *> twonorm(L*x) = sest
! 43: *> Then DLAIC1 computes sestpr, s, c such that
! 44: *> the vector
! 45: *> [ s*x ]
! 46: *> xhat = [ c ]
! 47: *> is an approximate singular vector of
! 48: *> [ L 0 ]
! 49: *> Lhat = [ w**T gamma ]
! 50: *> in the sense that
! 51: *> twonorm(Lhat*xhat) = sestpr.
! 52: *>
! 53: *> Depending on JOB, an estimate for the largest or smallest singular
! 54: *> value is computed.
! 55: *>
! 56: *> Note that [s c]**T and sestpr**2 is an eigenpair of the system
! 57: *>
! 58: *> diag(sest*sest, 0) + [alpha gamma] * [ alpha ]
! 59: *> [ gamma ]
! 60: *>
! 61: *> where alpha = x**T*w.
! 62: *> \endverbatim
! 63: *
! 64: * Arguments:
! 65: * ==========
! 66: *
! 67: *> \param[in] JOB
! 68: *> \verbatim
! 69: *> JOB is INTEGER
! 70: *> = 1: an estimate for the largest singular value is computed.
! 71: *> = 2: an estimate for the smallest singular value is computed.
! 72: *> \endverbatim
! 73: *>
! 74: *> \param[in] J
! 75: *> \verbatim
! 76: *> J is INTEGER
! 77: *> Length of X and W
! 78: *> \endverbatim
! 79: *>
! 80: *> \param[in] X
! 81: *> \verbatim
! 82: *> X is DOUBLE PRECISION array, dimension (J)
! 83: *> The j-vector x.
! 84: *> \endverbatim
! 85: *>
! 86: *> \param[in] SEST
! 87: *> \verbatim
! 88: *> SEST is DOUBLE PRECISION
! 89: *> Estimated singular value of j by j matrix L
! 90: *> \endverbatim
! 91: *>
! 92: *> \param[in] W
! 93: *> \verbatim
! 94: *> W is DOUBLE PRECISION array, dimension (J)
! 95: *> The j-vector w.
! 96: *> \endverbatim
! 97: *>
! 98: *> \param[in] GAMMA
! 99: *> \verbatim
! 100: *> GAMMA is DOUBLE PRECISION
! 101: *> The diagonal element gamma.
! 102: *> \endverbatim
! 103: *>
! 104: *> \param[out] SESTPR
! 105: *> \verbatim
! 106: *> SESTPR is DOUBLE PRECISION
! 107: *> Estimated singular value of (j+1) by (j+1) matrix Lhat.
! 108: *> \endverbatim
! 109: *>
! 110: *> \param[out] S
! 111: *> \verbatim
! 112: *> S is DOUBLE PRECISION
! 113: *> Sine needed in forming xhat.
! 114: *> \endverbatim
! 115: *>
! 116: *> \param[out] C
! 117: *> \verbatim
! 118: *> C is DOUBLE PRECISION
! 119: *> Cosine needed in forming xhat.
! 120: *> \endverbatim
! 121: *
! 122: * Authors:
! 123: * ========
! 124: *
! 125: *> \author Univ. of Tennessee
! 126: *> \author Univ. of California Berkeley
! 127: *> \author Univ. of Colorado Denver
! 128: *> \author NAG Ltd.
! 129: *
! 130: *> \date November 2011
! 131: *
! 132: *> \ingroup doubleOTHERauxiliary
! 133: *
! 134: * =====================================================================
1.1 bertrand 135: SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
136: *
1.9 ! bertrand 137: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 138: * -- LAPACK is a software package provided by Univ. of Tennessee, --
139: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 140: * November 2011
1.1 bertrand 141: *
142: * .. Scalar Arguments ..
143: INTEGER J, JOB
144: DOUBLE PRECISION C, GAMMA, S, SEST, SESTPR
145: * ..
146: * .. Array Arguments ..
147: DOUBLE PRECISION W( J ), X( J )
148: * ..
149: *
150: * =====================================================================
151: *
152: * .. Parameters ..
153: DOUBLE PRECISION ZERO, ONE, TWO
154: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
155: DOUBLE PRECISION HALF, FOUR
156: PARAMETER ( HALF = 0.5D0, FOUR = 4.0D0 )
157: * ..
158: * .. Local Scalars ..
159: DOUBLE PRECISION ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS,
160: $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2
161: * ..
162: * .. Intrinsic Functions ..
163: INTRINSIC ABS, MAX, SIGN, SQRT
164: * ..
165: * .. External Functions ..
166: DOUBLE PRECISION DDOT, DLAMCH
167: EXTERNAL DDOT, DLAMCH
168: * ..
169: * .. Executable Statements ..
170: *
171: EPS = DLAMCH( 'Epsilon' )
172: ALPHA = DDOT( J, X, 1, W, 1 )
173: *
174: ABSALP = ABS( ALPHA )
175: ABSGAM = ABS( GAMMA )
176: ABSEST = ABS( SEST )
177: *
178: IF( JOB.EQ.1 ) THEN
179: *
180: * Estimating largest singular value
181: *
182: * special cases
183: *
184: IF( SEST.EQ.ZERO ) THEN
185: S1 = MAX( ABSGAM, ABSALP )
186: IF( S1.EQ.ZERO ) THEN
187: S = ZERO
188: C = ONE
189: SESTPR = ZERO
190: ELSE
191: S = ALPHA / S1
192: C = GAMMA / S1
193: TMP = SQRT( S*S+C*C )
194: S = S / TMP
195: C = C / TMP
196: SESTPR = S1*TMP
197: END IF
198: RETURN
199: ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
200: S = ONE
201: C = ZERO
202: TMP = MAX( ABSEST, ABSALP )
203: S1 = ABSEST / TMP
204: S2 = ABSALP / TMP
205: SESTPR = TMP*SQRT( S1*S1+S2*S2 )
206: RETURN
207: ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
208: S1 = ABSGAM
209: S2 = ABSEST
210: IF( S1.LE.S2 ) THEN
211: S = ONE
212: C = ZERO
213: SESTPR = S2
214: ELSE
215: S = ZERO
216: C = ONE
217: SESTPR = S1
218: END IF
219: RETURN
220: ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
221: S1 = ABSGAM
222: S2 = ABSALP
223: IF( S1.LE.S2 ) THEN
224: TMP = S1 / S2
225: S = SQRT( ONE+TMP*TMP )
226: SESTPR = S2*S
227: C = ( GAMMA / S2 ) / S
228: S = SIGN( ONE, ALPHA ) / S
229: ELSE
230: TMP = S2 / S1
231: C = SQRT( ONE+TMP*TMP )
232: SESTPR = S1*C
233: S = ( ALPHA / S1 ) / C
234: C = SIGN( ONE, GAMMA ) / C
235: END IF
236: RETURN
237: ELSE
238: *
239: * normal case
240: *
241: ZETA1 = ALPHA / ABSEST
242: ZETA2 = GAMMA / ABSEST
243: *
244: B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
245: C = ZETA1*ZETA1
246: IF( B.GT.ZERO ) THEN
247: T = C / ( B+SQRT( B*B+C ) )
248: ELSE
249: T = SQRT( B*B+C ) - B
250: END IF
251: *
252: SINE = -ZETA1 / T
253: COSINE = -ZETA2 / ( ONE+T )
254: TMP = SQRT( SINE*SINE+COSINE*COSINE )
255: S = SINE / TMP
256: C = COSINE / TMP
257: SESTPR = SQRT( T+ONE )*ABSEST
258: RETURN
259: END IF
260: *
261: ELSE IF( JOB.EQ.2 ) THEN
262: *
263: * Estimating smallest singular value
264: *
265: * special cases
266: *
267: IF( SEST.EQ.ZERO ) THEN
268: SESTPR = ZERO
269: IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
270: SINE = ONE
271: COSINE = ZERO
272: ELSE
273: SINE = -GAMMA
274: COSINE = ALPHA
275: END IF
276: S1 = MAX( ABS( SINE ), ABS( COSINE ) )
277: S = SINE / S1
278: C = COSINE / S1
279: TMP = SQRT( S*S+C*C )
280: S = S / TMP
281: C = C / TMP
282: RETURN
283: ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
284: S = ZERO
285: C = ONE
286: SESTPR = ABSGAM
287: RETURN
288: ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
289: S1 = ABSGAM
290: S2 = ABSEST
291: IF( S1.LE.S2 ) THEN
292: S = ZERO
293: C = ONE
294: SESTPR = S1
295: ELSE
296: S = ONE
297: C = ZERO
298: SESTPR = S2
299: END IF
300: RETURN
301: ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
302: S1 = ABSGAM
303: S2 = ABSALP
304: IF( S1.LE.S2 ) THEN
305: TMP = S1 / S2
306: C = SQRT( ONE+TMP*TMP )
307: SESTPR = ABSEST*( TMP / C )
308: S = -( GAMMA / S2 ) / C
309: C = SIGN( ONE, ALPHA ) / C
310: ELSE
311: TMP = S2 / S1
312: S = SQRT( ONE+TMP*TMP )
313: SESTPR = ABSEST / S
314: C = ( ALPHA / S1 ) / S
315: S = -SIGN( ONE, GAMMA ) / S
316: END IF
317: RETURN
318: ELSE
319: *
320: * normal case
321: *
322: ZETA1 = ALPHA / ABSEST
323: ZETA2 = GAMMA / ABSEST
324: *
325: NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ),
326: $ ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 )
327: *
328: * See if root is closer to zero or to ONE
329: *
330: TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
331: IF( TEST.GE.ZERO ) THEN
332: *
333: * root is close to zero, compute directly
334: *
335: B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
336: C = ZETA2*ZETA2
337: T = C / ( B+SQRT( ABS( B*B-C ) ) )
338: SINE = ZETA1 / ( ONE-T )
339: COSINE = -ZETA2 / T
340: SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
341: ELSE
342: *
343: * root is closer to ONE, shift by that amount
344: *
345: B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
346: C = ZETA1*ZETA1
347: IF( B.GE.ZERO ) THEN
348: T = -C / ( B+SQRT( B*B+C ) )
349: ELSE
350: T = B - SQRT( B*B+C )
351: END IF
352: SINE = -ZETA1 / T
353: COSINE = -ZETA2 / ( ONE+T )
354: SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
355: END IF
356: TMP = SQRT( SINE*SINE+COSINE*COSINE )
357: S = SINE / TMP
358: C = COSINE / TMP
359: RETURN
360: *
361: END IF
362: END IF
363: RETURN
364: *
365: * End of DLAIC1
366: *
367: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaic1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack