LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
dlasq1.f
Go to the documentation of this file.
1 *> \brief \b DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
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 December 2016
105 *
106 *> \ingroup auxOTHERcomputational
107 *
108 * =====================================================================
109  SUBROUTINE dlasq1( N, D, E, WORK, INFO )
110 *
111 * -- LAPACK computational routine (version 3.7.0) --
112 * -- LAPACK is a software package provided by Univ. of Tennessee, --
113 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
114 * December 2016
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 = -1
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 )
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 )
218  END IF
219 *
220  RETURN
221 *
222 * End of DLASQ1
223 *
224  END
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:145
subroutine dlasq1(N, D, E, WORK, INFO)
DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
Definition: dlasq1.f:110
subroutine dlasq2(N, Z, INFO)
DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
Definition: dlasq2.f:114
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:90
subroutine dlas2(F, G, H, SSMIN, SSMAX)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: dlas2.f:109