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