1      SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
2     $                   U, LDU, C, LDC, WORK, INFO )
3*
4*  -- LAPACK auxiliary routine (version 3.0) --
5*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6*     Courant Institute, Argonne National Lab, and Rice University
7*     October 31, 1999
8*
9*     .. Scalar Arguments ..
10      CHARACTER          UPLO
11      INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
12*     ..
13*     .. Array Arguments ..
14      DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
15     $                   VT( LDVT, * ), WORK( * )
16*     ..
17*
18*  Purpose
19*  =======
20*
21*  DLASDQ computes the singular value decomposition (SVD) of a real
22*  (upper or lower) bidiagonal matrix with diagonal D and offdiagonal
23*  E, accumulating the transformations if desired. Letting B denote
24*  the input bidiagonal matrix, the algorithm computes orthogonal
25*  matrices Q and P such that B = Q * S * P' (P' denotes the transpose
26*  of P). The singular values S are overwritten on D.
27*
28*  The input matrix U  is changed to U  * Q  if desired.
29*  The input matrix VT is changed to P' * VT if desired.
30*  The input matrix C  is changed to Q' * C  if desired.
31*
32*  See "Computing  Small Singular Values of Bidiagonal Matrices With
33*  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
34*  LAPACK Working Note #3, for a detailed description of the algorithm.
35*
36*  Arguments
37*  =========
38*
39*  UPLO  (input) CHARACTER*1
40*        On entry, UPLO specifies whether the input bidiagonal matrix
41*        is upper or lower bidiagonal, and wether it is square are
42*        not.
43*           UPLO = 'U' or 'u'   B is upper bidiagonal.
44*           UPLO = 'L' or 'l'   B is lower bidiagonal.
45*
46*  SQRE  (input) INTEGER
47*        = 0: then the input matrix is N-by-N.
48*        = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and
49*             (N+1)-by-N if UPLU = 'L'.
50*
51*        The bidiagonal matrix has
52*        N = NL + NR + 1 rows and
53*        M = N + SQRE >= N columns.
54*
55*  N     (input) INTEGER
56*        On entry, N specifies the number of rows and columns
57*        in the matrix. N must be at least 0.
58*
59*  NCVT  (input) INTEGER
60*        On entry, NCVT specifies the number of columns of
61*        the matrix VT. NCVT must be at least 0.
62*
63*  NRU   (input) INTEGER
64*        On entry, NRU specifies the number of rows of
65*        the matrix U. NRU must be at least 0.
66*
67*  NCC   (input) INTEGER
68*        On entry, NCC specifies the number of columns of
69*        the matrix C. NCC must be at least 0.
70*
71*  D     (input/output) DOUBLE PRECISION array, dimension (N)
72*        On entry, D contains the diagonal entries of the
73*        bidiagonal matrix whose SVD is desired. On normal exit,
74*        D contains the singular values in ascending order.
75*
76*  E     (input/output) DOUBLE PRECISION array.
77*        dimension is (N-1) if SQRE = 0 and N if SQRE = 1.
78*        On entry, the entries of E contain the offdiagonal entries
79*        of the bidiagonal matrix whose SVD is desired. On normal
80*        exit, E will contain 0. If the algorithm does not converge,
81*        D and E will contain the diagonal and superdiagonal entries
82*        of a bidiagonal matrix orthogonally equivalent to the one
83*        given as input.
84*
85*  VT    (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
86*        On entry, contains a matrix which on exit has been
87*        premultiplied by P', dimension N-by-NCVT if SQRE = 0
88*        and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
89*
90*  LDVT  (input) INTEGER
91*        On entry, LDVT specifies the leading dimension of VT as
92*        declared in the calling (sub) program. LDVT must be at
93*        least 1. If NCVT is nonzero LDVT must also be at least N.
94*
95*  U     (input/output) DOUBLE PRECISION array, dimension (LDU, N)
96*        On entry, contains a  matrix which on exit has been
97*        postmultiplied by Q, dimension NRU-by-N if SQRE = 0
98*        and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
99*
100*  LDU   (input) INTEGER
101*        On entry, LDU  specifies the leading dimension of U as
102*        declared in the calling (sub) program. LDU must be at
103*        least max( 1, NRU ) .
104*
105*  C     (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
106*        On entry, contains an N-by-NCC matrix which on exit
107*        has been premultiplied by Q'  dimension N-by-NCC if SQRE = 0
108*        and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
109*
110*  LDC   (input) INTEGER
111*        On entry, LDC  specifies the leading dimension of C as
112*        declared in the calling (sub) program. LDC must be at
113*        least 1. If NCC is nonzero, LDC must also be at least N.
114*
115*  WORK  (workspace) DOUBLE PRECISION array, dimension (4*N)
116*        Workspace. Only referenced if one of NCVT, NRU, or NCC is
117*        nonzero, and if N is at least 2.
118*
119*  INFO  (output) INTEGER
120*        On exit, a value of 0 indicates a successful exit.
121*        If INFO < 0, argument number -INFO is illegal.
122*        If INFO > 0, the algorithm did not converge, and INFO
123*        specifies how many superdiagonals did not converge.
124*
125*  Further Details
126*  ===============
127*
128*  Based on contributions by
129*     Ming Gu and Huan Ren, Computer Science Division, University of
130*     California at Berkeley, USA
131*
132*  =====================================================================
133*
134*     .. Parameters ..
135      DOUBLE PRECISION   ZERO
136      PARAMETER          ( ZERO = 0.0D+0 )
137*     ..
138*     .. Local Scalars ..
139      LOGICAL            ROTATE
140      INTEGER            I, ISUB, IUPLO, J, NP1, SQRE1
141      DOUBLE PRECISION   CS, R, SMIN, SN
142*     ..
143*     .. External Subroutines ..
144      EXTERNAL           DBDSQR, DLARTG, DLASR, DSWAP, XERBLA
145*     ..
146*     .. External Functions ..
147      LOGICAL            LSAME
148      EXTERNAL           LSAME
149*     ..
150*     .. Intrinsic Functions ..
151      INTRINSIC          MAX
152*     ..
153*     .. Executable Statements ..
154*
155*     Test the input parameters.
156*
157      INFO = 0
158      IUPLO = 0
159      IF( LSAME( UPLO, 'U' ) )
160     $   IUPLO = 1
161      IF( LSAME( UPLO, 'L' ) )
162     $   IUPLO = 2
163      IF( IUPLO.EQ.0 ) THEN
164         INFO = -1
165      ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
166         INFO = -2
167      ELSE IF( N.LT.0 ) THEN
168         INFO = -3
169      ELSE IF( NCVT.LT.0 ) THEN
170         INFO = -4
171      ELSE IF( NRU.LT.0 ) THEN
172         INFO = -5
173      ELSE IF( NCC.LT.0 ) THEN
174         INFO = -6
175      ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
176     $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
177         INFO = -10
178      ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
179         INFO = -12
180      ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
181     $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
182         INFO = -14
183      END IF
184      IF( INFO.NE.0 ) THEN
185         CALL XERBLA( 'DLASDQ', -INFO )
186         RETURN
187      END IF
188      IF( N.EQ.0 )
189     $   RETURN
190*
191*     ROTATE is true if any singular vectors desired, false otherwise
192*
193      ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
194      NP1 = N + 1
195      SQRE1 = SQRE
196*
197*     If matrix non-square upper bidiagonal, rotate to be lower
198*     bidiagonal.  The rotations are on the right.
199*
200      IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN
201         DO 10 I = 1, N - 1
202            CALL DLARTG( D( I ), E( I ), CS, SN, R )
203            D( I ) = R
204            E( I ) = SN*D( I+1 )
205            D( I+1 ) = CS*D( I+1 )
206            IF( ROTATE ) THEN
207               WORK( I ) = CS
208               WORK( N+I ) = SN
209            END IF
210   10    CONTINUE
211         CALL DLARTG( D( N ), E( N ), CS, SN, R )
212         D( N ) = R
213         E( N ) = ZERO
214         IF( ROTATE ) THEN
215            WORK( N ) = CS
216            WORK( N+N ) = SN
217         END IF
218         IUPLO = 2
219         SQRE1 = 0
220*
221*        Update singular vectors if desired.
222*
223         IF( NCVT.GT.0 )
224     $      CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ),
225     $                  WORK( NP1 ), VT, LDVT )
226      END IF
227*
228*     If matrix lower bidiagonal, rotate to be upper bidiagonal
229*     by applying Givens rotations on the left.
230*
231      IF( IUPLO.EQ.2 ) THEN
232         DO 20 I = 1, N - 1
233            CALL DLARTG( D( I ), E( I ), CS, SN, R )
234            D( I ) = R
235            E( I ) = SN*D( I+1 )
236            D( I+1 ) = CS*D( I+1 )
237            IF( ROTATE ) THEN
238               WORK( I ) = CS
239               WORK( N+I ) = SN
240            END IF
241   20    CONTINUE
242*
243*        If matrix (N+1)-by-N lower bidiagonal, one additional
244*        rotation is needed.
245*
246         IF( SQRE1.EQ.1 ) THEN
247            CALL DLARTG( D( N ), E( N ), CS, SN, R )
248            D( N ) = R
249            IF( ROTATE ) THEN
250               WORK( N ) = CS
251               WORK( N+N ) = SN
252            END IF
253         END IF
254*
255*        Update singular vectors if desired.
256*
257         IF( NRU.GT.0 ) THEN
258            IF( SQRE1.EQ.0 ) THEN
259               CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),
260     $                     WORK( NP1 ), U, LDU )
261            ELSE
262               CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ),
263     $                     WORK( NP1 ), U, LDU )
264            END IF
265         END IF
266         IF( NCC.GT.0 ) THEN
267            IF( SQRE1.EQ.0 ) THEN
268               CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),
269     $                     WORK( NP1 ), C, LDC )
270            ELSE
271               CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),
272     $                     WORK( NP1 ), C, LDC )
273            END IF
274         END IF
275      END IF
276*
277*     Call DBDSQR to compute the SVD of the reduced real
278*     N-by-N upper bidiagonal matrix.
279*
280      CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C,
281     $             LDC, WORK, INFO )
282*
283*     Sort the singular values into ascending order (insertion sort on
284*     singular values, but only one transposition per singular vector)
285*
286      DO 40 I = 1, N
287*
288*        Scan for smallest D(I).
289*
290         ISUB = I
291         SMIN = D( I )
292         DO 30 J = I + 1, N
293            IF( D( J ).LT.SMIN ) THEN
294               ISUB = J
295               SMIN = D( J )
296            END IF
297   30    CONTINUE
298         IF( ISUB.NE.I ) THEN
299*
300*           Swap singular values and vectors.
301*
302            D( ISUB ) = D( I )
303            D( I ) = SMIN
304            IF( NCVT.GT.0 )
305     $         CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )
306            IF( NRU.GT.0 )
307     $         CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )
308            IF( NCC.GT.0 )
309     $         CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )
310         END IF
311   40 CONTINUE
312*
313      RETURN
314*
315*     End of DLASDQ
316*
317      END
318