Annotation of rpl/lapack/lapack/dlasq1.f, revision 1.9
1.8 bertrand 1: *> \brief \b DLASQ1
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
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">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, N
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION D( * ), E( * ), WORK( * )
28: * ..
29: *
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
92: *> = 3, termination criterion of outer while loop not met
93: *> (program created more than N unreduced blocks)
94: *> \endverbatim
95: *
96: * Authors:
97: * ========
98: *
99: *> \author Univ. of Tennessee
100: *> \author Univ. of California Berkeley
101: *> \author Univ. of Colorado Denver
102: *> \author NAG Ltd.
103: *
104: *> \date November 2011
1.1 bertrand 105: *
1.8 bertrand 106: *> \ingroup auxOTHERcomputational
1.1 bertrand 107: *
1.8 bertrand 108: * =====================================================================
109: SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
1.1 bertrand 110: *
1.8 bertrand 111: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 112: * -- LAPACK is a software package provided by Univ. of Tennessee, --
113: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 bertrand 114: * November 2011
1.1 bertrand 115: *
116: * .. Scalar Arguments ..
117: INTEGER INFO, N
118: * ..
119: * .. Array Arguments ..
120: DOUBLE PRECISION D( * ), E( * ), WORK( * )
121: * ..
122: *
123: * =====================================================================
124: *
125: * .. Parameters ..
126: DOUBLE PRECISION ZERO
127: PARAMETER ( ZERO = 0.0D0 )
128: * ..
129: * .. Local Scalars ..
130: INTEGER I, IINFO
131: DOUBLE PRECISION EPS, SCALE, SAFMIN, SIGMN, SIGMX
132: * ..
133: * .. External Subroutines ..
134: EXTERNAL DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA
135: * ..
136: * .. External Functions ..
137: DOUBLE PRECISION DLAMCH
138: EXTERNAL DLAMCH
139: * ..
140: * .. Intrinsic Functions ..
141: INTRINSIC ABS, MAX, SQRT
142: * ..
143: * .. Executable Statements ..
144: *
145: INFO = 0
146: IF( N.LT.0 ) THEN
147: INFO = -2
148: CALL XERBLA( 'DLASQ1', -INFO )
149: RETURN
150: ELSE IF( N.EQ.0 ) THEN
151: RETURN
152: ELSE IF( N.EQ.1 ) THEN
153: D( 1 ) = ABS( D( 1 ) )
154: RETURN
155: ELSE IF( N.EQ.2 ) THEN
156: CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
157: D( 1 ) = SIGMX
158: D( 2 ) = SIGMN
159: RETURN
160: END IF
161: *
162: * Estimate the largest singular value.
163: *
164: SIGMX = ZERO
165: DO 10 I = 1, N - 1
166: D( I ) = ABS( D( I ) )
167: SIGMX = MAX( SIGMX, ABS( E( I ) ) )
168: 10 CONTINUE
169: D( N ) = ABS( D( N ) )
170: *
171: * Early return if SIGMX is zero (matrix is already diagonal).
172: *
173: IF( SIGMX.EQ.ZERO ) THEN
174: CALL DLASRT( 'D', N, D, IINFO )
175: RETURN
176: END IF
177: *
178: DO 20 I = 1, N
179: SIGMX = MAX( SIGMX, D( I ) )
180: 20 CONTINUE
181: *
182: * Copy D and E into WORK (in the Z format) and scale (squaring the
183: * input data makes scaling by a power of the radix pointless).
184: *
185: EPS = DLAMCH( 'Precision' )
186: SAFMIN = DLAMCH( 'Safe minimum' )
187: SCALE = SQRT( EPS / SAFMIN )
188: CALL DCOPY( N, D, 1, WORK( 1 ), 2 )
189: CALL DCOPY( N-1, E, 1, WORK( 2 ), 2 )
190: CALL DLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1,
191: $ IINFO )
192: *
193: * Compute the q's and e's.
194: *
195: DO 30 I = 1, 2*N - 1
196: WORK( I ) = WORK( I )**2
197: 30 CONTINUE
198: WORK( 2*N ) = ZERO
199: *
200: CALL DLASQ2( N, WORK, INFO )
201: *
202: IF( INFO.EQ.0 ) THEN
203: DO 40 I = 1, N
204: D( I ) = SQRT( WORK( I ) )
205: 40 CONTINUE
206: CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
1.8 bertrand 207: ELSE IF( INFO.EQ.2 ) THEN
208: *
209: * Maximum number of iterations exceeded. Move data from WORK
210: * into D and E so the calling subroutine can try to finish
211: *
212: DO I = 1, N
213: D( I ) = SQRT( WORK( 2*I-1 ) )
214: E( I ) = SQRT( WORK( 2*I ) )
215: END DO
216: CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
217: CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, E, N, IINFO )
1.1 bertrand 218: END IF
219: *
220: RETURN
221: *
222: * End of DLASQ1
223: *
224: 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