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
225