1*> \brief \b SLASDQ 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 SLASDQ + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasdq.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasdq.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasdq.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLASDQ( 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*       REAL               C( LDC, * ), D( * ), E( * ), U( LDU, * ),
30*      $                   VT( LDVT, * ), WORK( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> SLASDQ 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 whether 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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*> \ingroup OTHERauxiliary
201*
202*> \par Contributors:
203*  ==================
204*>
205*>     Ming Gu and Huan Ren, Computer Science Division, University of
206*>     California at Berkeley, USA
207*>
208*  =====================================================================
209      SUBROUTINE SLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
210     $                   U, LDU, C, LDC, WORK, INFO )
211*
212*  -- LAPACK auxiliary routine --
213*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
214*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215*
216*     .. Scalar Arguments ..
217      CHARACTER          UPLO
218      INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
219*     ..
220*     .. Array Arguments ..
221      REAL               C( LDC, * ), D( * ), E( * ), U( LDU, * ),
222     $                   VT( LDVT, * ), WORK( * )
223*     ..
224*
225*  =====================================================================
226*
227*     .. Parameters ..
228      REAL               ZERO
229      PARAMETER          ( ZERO = 0.0E+0 )
230*     ..
231*     .. Local Scalars ..
232      LOGICAL            ROTATE
233      INTEGER            I, ISUB, IUPLO, J, NP1, SQRE1
234      REAL               CS, R, SMIN, SN
235*     ..
236*     .. External Subroutines ..
237      EXTERNAL           SBDSQR, SLARTG, SLASR, SSWAP, XERBLA
238*     ..
239*     .. External Functions ..
240      LOGICAL            LSAME
241      EXTERNAL           LSAME
242*     ..
243*     .. Intrinsic Functions ..
244      INTRINSIC          MAX
245*     ..
246*     .. Executable Statements ..
247*
248*     Test the input parameters.
249*
250      INFO = 0
251      IUPLO = 0
252      IF( LSAME( UPLO, 'U' ) )
253     $   IUPLO = 1
254      IF( LSAME( UPLO, 'L' ) )
255     $   IUPLO = 2
256      IF( IUPLO.EQ.0 ) THEN
257         INFO = -1
258      ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
259         INFO = -2
260      ELSE IF( N.LT.0 ) THEN
261         INFO = -3
262      ELSE IF( NCVT.LT.0 ) THEN
263         INFO = -4
264      ELSE IF( NRU.LT.0 ) THEN
265         INFO = -5
266      ELSE IF( NCC.LT.0 ) THEN
267         INFO = -6
268      ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
269     $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
270         INFO = -10
271      ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
272         INFO = -12
273      ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
274     $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
275         INFO = -14
276      END IF
277      IF( INFO.NE.0 ) THEN
278         CALL XERBLA( 'SLASDQ', -INFO )
279         RETURN
280      END IF
281      IF( N.EQ.0 )
282     $   RETURN
283*
284*     ROTATE is true if any singular vectors desired, false otherwise
285*
286      ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
287      NP1 = N + 1
288      SQRE1 = SQRE
289*
290*     If matrix non-square upper bidiagonal, rotate to be lower
291*     bidiagonal.  The rotations are on the right.
292*
293      IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN
294         DO 10 I = 1, N - 1
295            CALL SLARTG( D( I ), E( I ), CS, SN, R )
296            D( I ) = R
297            E( I ) = SN*D( I+1 )
298            D( I+1 ) = CS*D( I+1 )
299            IF( ROTATE ) THEN
300               WORK( I ) = CS
301               WORK( N+I ) = SN
302            END IF
303   10    CONTINUE
304         CALL SLARTG( D( N ), E( N ), CS, SN, R )
305         D( N ) = R
306         E( N ) = ZERO
307         IF( ROTATE ) THEN
308            WORK( N ) = CS
309            WORK( N+N ) = SN
310         END IF
311         IUPLO = 2
312         SQRE1 = 0
313*
314*        Update singular vectors if desired.
315*
316         IF( NCVT.GT.0 )
317     $      CALL SLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ),
318     $                  WORK( NP1 ), VT, LDVT )
319      END IF
320*
321*     If matrix lower bidiagonal, rotate to be upper bidiagonal
322*     by applying Givens rotations on the left.
323*
324      IF( IUPLO.EQ.2 ) THEN
325         DO 20 I = 1, N - 1
326            CALL SLARTG( D( I ), E( I ), CS, SN, R )
327            D( I ) = R
328            E( I ) = SN*D( I+1 )
329            D( I+1 ) = CS*D( I+1 )
330            IF( ROTATE ) THEN
331               WORK( I ) = CS
332               WORK( N+I ) = SN
333            END IF
334   20    CONTINUE
335*
336*        If matrix (N+1)-by-N lower bidiagonal, one additional
337*        rotation is needed.
338*
339         IF( SQRE1.EQ.1 ) THEN
340            CALL SLARTG( D( N ), E( N ), CS, SN, R )
341            D( N ) = R
342            IF( ROTATE ) THEN
343               WORK( N ) = CS
344               WORK( N+N ) = SN
345            END IF
346         END IF
347*
348*        Update singular vectors if desired.
349*
350         IF( NRU.GT.0 ) THEN
351            IF( SQRE1.EQ.0 ) THEN
352               CALL SLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),
353     $                     WORK( NP1 ), U, LDU )
354            ELSE
355               CALL SLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ),
356     $                     WORK( NP1 ), U, LDU )
357            END IF
358         END IF
359         IF( NCC.GT.0 ) THEN
360            IF( SQRE1.EQ.0 ) THEN
361               CALL SLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),
362     $                     WORK( NP1 ), C, LDC )
363            ELSE
364               CALL SLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),
365     $                     WORK( NP1 ), C, LDC )
366            END IF
367         END IF
368      END IF
369*
370*     Call SBDSQR to compute the SVD of the reduced real
371*     N-by-N upper bidiagonal matrix.
372*
373      CALL SBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C,
374     $             LDC, WORK, INFO )
375*
376*     Sort the singular values into ascending order (insertion sort on
377*     singular values, but only one transposition per singular vector)
378*
379      DO 40 I = 1, N
380*
381*        Scan for smallest D(I).
382*
383         ISUB = I
384         SMIN = D( I )
385         DO 30 J = I + 1, N
386            IF( D( J ).LT.SMIN ) THEN
387               ISUB = J
388               SMIN = D( J )
389            END IF
390   30    CONTINUE
391         IF( ISUB.NE.I ) THEN
392*
393*           Swap singular values and vectors.
394*
395            D( ISUB ) = D( I )
396            D( I ) = SMIN
397            IF( NCVT.GT.0 )
398     $         CALL SSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )
399            IF( NRU.GT.0 )
400     $         CALL SSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )
401            IF( NCC.GT.0 )
402     $         CALL SSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )
403         END IF
404   40 CONTINUE
405*
406      RETURN
407*
408*     End of SLASDQ
409*
410      END
411