Annotation of rpl/lapack/lapack/dlasq1.f, revision 1.19
1.11 bertrand 1: *> \brief \b DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
1.8 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.16 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.8 bertrand 7: *
8: *> \htmlonly
1.16 bertrand 9: *> Download DLASQ1 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq1.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq1.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq1.f">
1.8 bertrand 15: *> [TXT]</a>
1.16 bertrand 16: *> \endhtmlonly
1.8 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
1.16 bertrand 22: *
1.8 bertrand 23: * .. Scalar Arguments ..
24: * INTEGER INFO, N
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION D( * ), E( * ), WORK( * )
28: * ..
1.16 bertrand 29: *
1.8 bertrand 30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> DLASQ1 computes the singular values of a real N-by-N bidiagonal
37: *> matrix with diagonal D and off-diagonal E. The singular values
38: *> are computed to high relative accuracy, in the absence of
39: *> denormalization, underflow and overflow. The algorithm was first
40: *> presented in
41: *>
42: *> "Accurate singular values and differential qd algorithms" by K. V.
43: *> Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
44: *> 1994,
45: *>
46: *> and the present implementation is described in "An implementation of
47: *> the dqds Algorithm (Positive Case)", LAPACK Working Note.
48: *> \endverbatim
49: *
50: * Arguments:
51: * ==========
52: *
53: *> \param[in] N
54: *> \verbatim
55: *> N is INTEGER
56: *> The number of rows and columns in the matrix. N >= 0.
57: *> \endverbatim
58: *>
59: *> \param[in,out] D
60: *> \verbatim
61: *> D is DOUBLE PRECISION array, dimension (N)
62: *> On entry, D contains the diagonal elements of the
63: *> bidiagonal matrix whose SVD is desired. On normal exit,
64: *> D contains the singular values in decreasing order.
65: *> \endverbatim
66: *>
67: *> \param[in,out] E
68: *> \verbatim
69: *> E is DOUBLE PRECISION array, dimension (N)
70: *> On entry, elements E(1:N-1) contain the off-diagonal elements
71: *> of the bidiagonal matrix whose SVD is desired.
72: *> On exit, E is overwritten.
73: *> \endverbatim
74: *>
75: *> \param[out] WORK
76: *> \verbatim
77: *> WORK is DOUBLE PRECISION array, dimension (4*N)
78: *> \endverbatim
79: *>
80: *> \param[out] INFO
81: *> \verbatim
82: *> INFO is INTEGER
83: *> = 0: successful exit
84: *> < 0: if INFO = -i, the i-th argument had an illegal value
85: *> > 0: the algorithm failed
86: *> = 1, a split was marked by a positive value in E
87: *> = 2, current block of Z not diagonalized after 100*N
88: *> iterations (in inner while loop) On exit D and E
89: *> represent a matrix with the same singular values
90: *> which the calling subroutine could use to finish the
91: *> computation, or even feed back into DLASQ1
1.16 bertrand 92: *> = 3, termination criterion of outer while loop not met
1.8 bertrand 93: *> (program created more than N unreduced blocks)
94: *> \endverbatim
95: *
96: * Authors:
97: * ========
98: *
1.16 bertrand 99: *> \author Univ. of Tennessee
100: *> \author Univ. of California Berkeley
101: *> \author Univ. of Colorado Denver
102: *> \author NAG Ltd.
1.8 bertrand 103: *
104: *> \ingroup auxOTHERcomputational
1.1 bertrand 105: *
1.8 bertrand 106: * =====================================================================
107: SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
1.1 bertrand 108: *
1.19 ! bertrand 109: * -- LAPACK computational routine --
1.1 bertrand 110: * -- LAPACK is a software package provided by Univ. of Tennessee, --
111: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
112: *
113: * .. Scalar Arguments ..
114: INTEGER INFO, N
115: * ..
116: * .. Array Arguments ..
117: DOUBLE PRECISION D( * ), E( * ), WORK( * )
118: * ..
119: *
120: * =====================================================================
121: *
122: * .. Parameters ..
123: DOUBLE PRECISION ZERO
124: PARAMETER ( ZERO = 0.0D0 )
125: * ..
126: * .. Local Scalars ..
127: INTEGER I, IINFO
128: DOUBLE PRECISION EPS, SCALE, SAFMIN, SIGMN, SIGMX
129: * ..
130: * .. External Subroutines ..
131: EXTERNAL DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA
132: * ..
133: * .. External Functions ..
134: DOUBLE PRECISION DLAMCH
135: EXTERNAL DLAMCH
136: * ..
137: * .. Intrinsic Functions ..
138: INTRINSIC ABS, MAX, SQRT
139: * ..
140: * .. Executable Statements ..
141: *
142: INFO = 0
143: IF( N.LT.0 ) THEN
1.14 bertrand 144: INFO = -1
1.1 bertrand 145: CALL XERBLA( 'DLASQ1', -INFO )
146: RETURN
147: ELSE IF( N.EQ.0 ) THEN
148: RETURN
149: ELSE IF( N.EQ.1 ) THEN
150: D( 1 ) = ABS( D( 1 ) )
151: RETURN
152: ELSE IF( N.EQ.2 ) THEN
153: CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
154: D( 1 ) = SIGMX
155: D( 2 ) = SIGMN
156: RETURN
157: END IF
158: *
159: * Estimate the largest singular value.
160: *
161: SIGMX = ZERO
162: DO 10 I = 1, N - 1
163: D( I ) = ABS( D( I ) )
164: SIGMX = MAX( SIGMX, ABS( E( I ) ) )
165: 10 CONTINUE
166: D( N ) = ABS( D( N ) )
167: *
168: * Early return if SIGMX is zero (matrix is already diagonal).
169: *
170: IF( SIGMX.EQ.ZERO ) THEN
171: CALL DLASRT( 'D', N, D, IINFO )
172: RETURN
173: END IF
174: *
175: DO 20 I = 1, N
176: SIGMX = MAX( SIGMX, D( I ) )
177: 20 CONTINUE
178: *
179: * Copy D and E into WORK (in the Z format) and scale (squaring the
180: * input data makes scaling by a power of the radix pointless).
181: *
182: EPS = DLAMCH( 'Precision' )
183: SAFMIN = DLAMCH( 'Safe minimum' )
184: SCALE = SQRT( EPS / SAFMIN )
185: CALL DCOPY( N, D, 1, WORK( 1 ), 2 )
186: CALL DCOPY( N-1, E, 1, WORK( 2 ), 2 )
187: CALL DLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1,
188: $ IINFO )
1.16 bertrand 189: *
1.1 bertrand 190: * Compute the q's and e's.
191: *
192: DO 30 I = 1, 2*N - 1
193: WORK( I ) = WORK( I )**2
194: 30 CONTINUE
195: WORK( 2*N ) = ZERO
196: *
197: CALL DLASQ2( N, WORK, INFO )
198: *
199: IF( INFO.EQ.0 ) THEN
200: DO 40 I = 1, N
201: D( I ) = SQRT( WORK( I ) )
202: 40 CONTINUE
203: CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
1.8 bertrand 204: ELSE IF( INFO.EQ.2 ) THEN
205: *
206: * Maximum number of iterations exceeded. Move data from WORK
207: * into D and E so the calling subroutine can try to finish
208: *
209: DO I = 1, N
210: D( I ) = SQRT( WORK( 2*I-1 ) )
211: E( I ) = SQRT( WORK( 2*I ) )
212: END DO
213: CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
214: CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, E, N, IINFO )
1.1 bertrand 215: END IF
216: *
217: RETURN
218: *
219: * End of DLASQ1
220: *
221: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack