1*> \brief \b DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLASDQ + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasdq.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasdq.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasdq.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
22*                          U, LDU, C, LDC, WORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          UPLO
26*       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
27*       ..
28*       .. Array Arguments ..
29*       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
30*      $                   VT( LDVT, * ), WORK( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> DLASDQ computes the singular value decomposition (SVD) of a real
40*> (upper or lower) bidiagonal matrix with diagonal D and offdiagonal
41*> E, accumulating the transformations if desired. Letting B denote
42*> the input bidiagonal matrix, the algorithm computes orthogonal
43*> matrices Q and P such that B = Q * S * P**T (P**T denotes the transpose
44*> of P). The singular values S are overwritten on D.
45*>
46*> The input matrix U  is changed to U  * Q  if desired.
47*> The input matrix VT is changed to P**T * VT if desired.
48*> The input matrix C  is changed to Q**T * C  if desired.
49*>
50*> See "Computing  Small Singular Values of Bidiagonal Matrices With
51*> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
52*> LAPACK Working Note #3, for a detailed description of the algorithm.
53*> \endverbatim
54*
55*  Arguments:
56*  ==========
57*
58*> \param[in] UPLO
59*> \verbatim
60*>          UPLO is CHARACTER*1
61*>        On entry, UPLO specifies whether the input bidiagonal matrix
62*>        is upper or lower bidiagonal, and wether it is square are
63*>        not.
64*>           UPLO = 'U' or 'u'   B is upper bidiagonal.
65*>           UPLO = 'L' or 'l'   B is lower bidiagonal.
66*> \endverbatim
67*>
68*> \param[in] SQRE
69*> \verbatim
70*>          SQRE is INTEGER
71*>        = 0: then the input matrix is N-by-N.
72*>        = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and
73*>             (N+1)-by-N if UPLU = 'L'.
74*>
75*>        The bidiagonal matrix has
76*>        N = NL + NR + 1 rows and
77*>        M = N + SQRE >= N columns.
78*> \endverbatim
79*>
80*> \param[in] N
81*> \verbatim
82*>          N is INTEGER
83*>        On entry, N specifies the number of rows and columns
84*>        in the matrix. N must be at least 0.
85*> \endverbatim
86*>
87*> \param[in] NCVT
88*> \verbatim
89*>          NCVT is INTEGER
90*>        On entry, NCVT specifies the number of columns of
91*>        the matrix VT. NCVT must be at least 0.
92*> \endverbatim
93*>
94*> \param[in] NRU
95*> \verbatim
96*>          NRU is INTEGER
97*>        On entry, NRU specifies the number of rows of
98*>        the matrix U. NRU must be at least 0.
99*> \endverbatim
100*>
101*> \param[in] NCC
102*> \verbatim
103*>          NCC is INTEGER
104*>        On entry, NCC specifies the number of columns of
105*>        the matrix C. NCC must be at least 0.
106*> \endverbatim
107*>
108*> \param[in,out] D
109*> \verbatim
110*>          D is DOUBLE PRECISION array, dimension (N)
111*>        On entry, D contains the diagonal entries of the
112*>        bidiagonal matrix whose SVD is desired. On normal exit,
113*>        D contains the singular values in ascending order.
114*> \endverbatim
115*>
116*> \param[in,out] E
117*> \verbatim
118*>          E is DOUBLE PRECISION array.
119*>        dimension is (N-1) if SQRE = 0 and N if SQRE = 1.
120*>        On entry, the entries of E contain the offdiagonal entries
121*>        of the bidiagonal matrix whose SVD is desired. On normal
122*>        exit, E will contain 0. If the algorithm does not converge,
123*>        D and E will contain the diagonal and superdiagonal entries
124*>        of a bidiagonal matrix orthogonally equivalent to the one
125*>        given as input.
126*> \endverbatim
127*>
128*> \param[in,out] VT
129*> \verbatim
130*>          VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)
131*>        On entry, contains a matrix which on exit has been
132*>        premultiplied by P**T, dimension N-by-NCVT if SQRE = 0
133*>        and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
134*> \endverbatim
135*>
136*> \param[in] LDVT
137*> \verbatim
138*>          LDVT is INTEGER
139*>        On entry, LDVT specifies the leading dimension of VT as
140*>        declared in the calling (sub) program. LDVT must be at
141*>        least 1. If NCVT is nonzero LDVT must also be at least N.
142*> \endverbatim
143*>
144*> \param[in,out] U
145*> \verbatim
146*>          U is DOUBLE PRECISION array, dimension (LDU, N)
147*>        On entry, contains a  matrix which on exit has been
148*>        postmultiplied by Q, dimension NRU-by-N if SQRE = 0
149*>        and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
150*> \endverbatim
151*>
152*> \param[in] LDU
153*> \verbatim
154*>          LDU is INTEGER
155*>        On entry, LDU  specifies the leading dimension of U as
156*>        declared in the calling (sub) program. LDU must be at
157*>        least max( 1, NRU ) .
158*> \endverbatim
159*>
160*> \param[in,out] C
161*> \verbatim
162*>          C is DOUBLE PRECISION array, dimension (LDC, NCC)
163*>        On entry, contains an N-by-NCC matrix which on exit
164*>        has been premultiplied by Q**T  dimension N-by-NCC if SQRE = 0
165*>        and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
166*> \endverbatim
167*>
168*> \param[in] LDC
169*> \verbatim
170*>          LDC is INTEGER
171*>        On entry, LDC  specifies the leading dimension of C as
172*>        declared in the calling (sub) program. LDC must be at
173*>        least 1. If NCC is nonzero, LDC must also be at least N.
174*> \endverbatim
175*>
176*> \param[out] WORK
177*> \verbatim
178*>          WORK is DOUBLE PRECISION array, dimension (4*N)
179*>        Workspace. Only referenced if one of NCVT, NRU, or NCC is
180*>        nonzero, and if N is at least 2.
181*> \endverbatim
182*>
183*> \param[out] INFO
184*> \verbatim
185*>          INFO is INTEGER
186*>        On exit, a value of 0 indicates a successful exit.
187*>        If INFO < 0, argument number -INFO is illegal.
188*>        If INFO > 0, the algorithm did not converge, and INFO
189*>        specifies how many superdiagonals did not converge.
190*> \endverbatim
191*
192*  Authors:
193*  ========
194*
195*> \author Univ. of Tennessee
196*> \author Univ. of California Berkeley
197*> \author Univ. of Colorado Denver
198*> \author NAG Ltd.
199*
200*> \date September 2012
201*
202*> \ingroup auxOTHERauxiliary
203*
204*> \par Contributors:
205*  ==================
206*>
207*>     Ming Gu and Huan Ren, Computer Science Division, University of
208*>     California at Berkeley, USA
209*>
210*  =====================================================================
211      SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
212     $                   U, LDU, C, LDC, WORK, INFO )
213*
214*  -- LAPACK auxiliary routine (version 3.4.2) --
215*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
216*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217*     September 2012
218*
219*     .. Scalar Arguments ..
220      CHARACTER          UPLO
221      INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
222*     ..
223*     .. Array Arguments ..
224      DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
225     $                   VT( LDVT, * ), WORK( * )
226*     ..
227*
228*  =====================================================================
229*
230*     .. Parameters ..
231      DOUBLE PRECISION   ZERO
232      PARAMETER          ( ZERO = 0.0D+0 )
233*     ..
234*     .. Local Scalars ..
235      LOGICAL            ROTATE
236      INTEGER            I, ISUB, IUPLO, J, NP1, SQRE1
237      DOUBLE PRECISION   CS, R, SMIN, SN
238*     ..
239*     .. External Subroutines ..
240      EXTERNAL           DBDSQR, DLARTG, DLASR, DSWAP, XERBLA
241*     ..
242*     .. External Functions ..
243      LOGICAL            LSAME
244      EXTERNAL           LSAME
245*     ..
246*     .. Intrinsic Functions ..
247      INTRINSIC          MAX
248*     ..
249*     .. Executable Statements ..
250*
251*     Test the input parameters.
252*
253      INFO = 0
254      IUPLO = 0
255      IF( LSAME( UPLO, 'U' ) )
256     $   IUPLO = 1
257      IF( LSAME( UPLO, 'L' ) )
258     $   IUPLO = 2
259      IF( IUPLO.EQ.0 ) THEN
260         INFO = -1
261      ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
262         INFO = -2
263      ELSE IF( N.LT.0 ) THEN
264         INFO = -3
265      ELSE IF( NCVT.LT.0 ) THEN
266         INFO = -4
267      ELSE IF( NRU.LT.0 ) THEN
268         INFO = -5
269      ELSE IF( NCC.LT.0 ) THEN
270         INFO = -6
271      ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
272     $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
273         INFO = -10
274      ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
275         INFO = -12
276      ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
277     $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
278         INFO = -14
279      END IF
280      IF( INFO.NE.0 ) THEN
281         CALL XERBLA( 'DLASDQ', -INFO )
282         RETURN
283      END IF
284      IF( N.EQ.0 )
285     $   RETURN
286*
287*     ROTATE is true if any singular vectors desired, false otherwise
288*
289      ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
290      NP1 = N + 1
291      SQRE1 = SQRE
292*
293*     If matrix non-square upper bidiagonal, rotate to be lower
294*     bidiagonal.  The rotations are on the right.
295*
296      IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN
297         DO 10 I = 1, N - 1
298            CALL DLARTG( D( I ), E( I ), CS, SN, R )
299            D( I ) = R
300            E( I ) = SN*D( I+1 )
301            D( I+1 ) = CS*D( I+1 )
302            IF( ROTATE ) THEN
303               WORK( I ) = CS
304               WORK( N+I ) = SN
305            END IF
306   10    CONTINUE
307         CALL DLARTG( D( N ), E( N ), CS, SN, R )
308         D( N ) = R
309         E( N ) = ZERO
310         IF( ROTATE ) THEN
311            WORK( N ) = CS
312            WORK( N+N ) = SN
313         END IF
314         IUPLO = 2
315         SQRE1 = 0
316*
317*        Update singular vectors if desired.
318*
319         IF( NCVT.GT.0 )
320     $      CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ),
321     $                  WORK( NP1 ), VT, LDVT )
322      END IF
323*
324*     If matrix lower bidiagonal, rotate to be upper bidiagonal
325*     by applying Givens rotations on the left.
326*
327      IF( IUPLO.EQ.2 ) THEN
328         DO 20 I = 1, N - 1
329            CALL DLARTG( D( I ), E( I ), CS, SN, R )
330            D( I ) = R
331            E( I ) = SN*D( I+1 )
332            D( I+1 ) = CS*D( I+1 )
333            IF( ROTATE ) THEN
334               WORK( I ) = CS
335               WORK( N+I ) = SN
336            END IF
337   20    CONTINUE
338*
339*        If matrix (N+1)-by-N lower bidiagonal, one additional
340*        rotation is needed.
341*
342         IF( SQRE1.EQ.1 ) THEN
343            CALL DLARTG( D( N ), E( N ), CS, SN, R )
344            D( N ) = R
345            IF( ROTATE ) THEN
346               WORK( N ) = CS
347               WORK( N+N ) = SN
348            END IF
349         END IF
350*
351*        Update singular vectors if desired.
352*
353         IF( NRU.GT.0 ) THEN
354            IF( SQRE1.EQ.0 ) THEN
355               CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),
356     $                     WORK( NP1 ), U, LDU )
357            ELSE
358               CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ),
359     $                     WORK( NP1 ), U, LDU )
360            END IF
361         END IF
362         IF( NCC.GT.0 ) THEN
363            IF( SQRE1.EQ.0 ) THEN
364               CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),
365     $                     WORK( NP1 ), C, LDC )
366            ELSE
367               CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),
368     $                     WORK( NP1 ), C, LDC )
369            END IF
370         END IF
371      END IF
372*
373*     Call DBDSQR to compute the SVD of the reduced real
374*     N-by-N upper bidiagonal matrix.
375*
376      CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C,
377     $             LDC, WORK, INFO )
378*
379*     Sort the singular values into ascending order (insertion sort on
380*     singular values, but only one transposition per singular vector)
381*
382      DO 40 I = 1, N
383*
384*        Scan for smallest D(I).
385*
386         ISUB = I
387         SMIN = D( I )
388         DO 30 J = I + 1, N
389            IF( D( J ).LT.SMIN ) THEN
390               ISUB = J
391               SMIN = D( J )
392            END IF
393   30    CONTINUE
394         IF( ISUB.NE.I ) THEN
395*
396*           Swap singular values and vectors.
397*
398            D( ISUB ) = D( I )
399            D( I ) = SMIN
400            IF( NCVT.GT.0 )
401     $         CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )
402            IF( NRU.GT.0 )
403     $         CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )
404            IF( NCC.GT.0 )
405     $         CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )
406         END IF
407   40 CONTINUE
408*
409      RETURN
410*
411*     End of DLASDQ
412*
413      END
414