001:       SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
002:      $                   U, LDU, C, LDC, WORK, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          UPLO
010:       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
014:      $                   VT( LDVT, * ), WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DLASDQ computes the singular value decomposition (SVD) of a real
021: *  (upper or lower) bidiagonal matrix with diagonal D and offdiagonal
022: *  E, accumulating the transformations if desired. Letting B denote
023: *  the input bidiagonal matrix, the algorithm computes orthogonal
024: *  matrices Q and P such that B = Q * S * P' (P' denotes the transpose
025: *  of P). The singular values S are overwritten on D.
026: *
027: *  The input matrix U  is changed to U  * Q  if desired.
028: *  The input matrix VT is changed to P' * VT if desired.
029: *  The input matrix C  is changed to Q' * C  if desired.
030: *
031: *  See "Computing  Small Singular Values of Bidiagonal Matrices With
032: *  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
033: *  LAPACK Working Note #3, for a detailed description of the algorithm.
034: *
035: *  Arguments
036: *  =========
037: *
038: *  UPLO  (input) CHARACTER*1
039: *        On entry, UPLO specifies whether the input bidiagonal matrix
040: *        is upper or lower bidiagonal, and wether it is square are
041: *        not.
042: *           UPLO = 'U' or 'u'   B is upper bidiagonal.
043: *           UPLO = 'L' or 'l'   B is lower bidiagonal.
044: *
045: *  SQRE  (input) INTEGER
046: *        = 0: then the input matrix is N-by-N.
047: *        = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and
048: *             (N+1)-by-N if UPLU = 'L'.
049: *
050: *        The bidiagonal matrix has
051: *        N = NL + NR + 1 rows and
052: *        M = N + SQRE >= N columns.
053: *
054: *  N     (input) INTEGER
055: *        On entry, N specifies the number of rows and columns
056: *        in the matrix. N must be at least 0.
057: *
058: *  NCVT  (input) INTEGER
059: *        On entry, NCVT specifies the number of columns of
060: *        the matrix VT. NCVT must be at least 0.
061: *
062: *  NRU   (input) INTEGER
063: *        On entry, NRU specifies the number of rows of
064: *        the matrix U. NRU must be at least 0.
065: *
066: *  NCC   (input) INTEGER
067: *        On entry, NCC specifies the number of columns of
068: *        the matrix C. NCC must be at least 0.
069: *
070: *  D     (input/output) DOUBLE PRECISION array, dimension (N)
071: *        On entry, D contains the diagonal entries of the
072: *        bidiagonal matrix whose SVD is desired. On normal exit,
073: *        D contains the singular values in ascending order.
074: *
075: *  E     (input/output) DOUBLE PRECISION array.
076: *        dimension is (N-1) if SQRE = 0 and N if SQRE = 1.
077: *        On entry, the entries of E contain the offdiagonal entries
078: *        of the bidiagonal matrix whose SVD is desired. On normal
079: *        exit, E will contain 0. If the algorithm does not converge,
080: *        D and E will contain the diagonal and superdiagonal entries
081: *        of a bidiagonal matrix orthogonally equivalent to the one
082: *        given as input.
083: *
084: *  VT    (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
085: *        On entry, contains a matrix which on exit has been
086: *        premultiplied by P', dimension N-by-NCVT if SQRE = 0
087: *        and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
088: *
089: *  LDVT  (input) INTEGER
090: *        On entry, LDVT specifies the leading dimension of VT as
091: *        declared in the calling (sub) program. LDVT must be at
092: *        least 1. If NCVT is nonzero LDVT must also be at least N.
093: *
094: *  U     (input/output) DOUBLE PRECISION array, dimension (LDU, N)
095: *        On entry, contains a  matrix which on exit has been
096: *        postmultiplied by Q, dimension NRU-by-N if SQRE = 0
097: *        and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
098: *
099: *  LDU   (input) INTEGER
100: *        On entry, LDU  specifies the leading dimension of U as
101: *        declared in the calling (sub) program. LDU must be at
102: *        least max( 1, NRU ) .
103: *
104: *  C     (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
105: *        On entry, contains an N-by-NCC matrix which on exit
106: *        has been premultiplied by Q'  dimension N-by-NCC if SQRE = 0
107: *        and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
108: *
109: *  LDC   (input) INTEGER
110: *        On entry, LDC  specifies the leading dimension of C as
111: *        declared in the calling (sub) program. LDC must be at
112: *        least 1. If NCC is nonzero, LDC must also be at least N.
113: *
114: *  WORK  (workspace) DOUBLE PRECISION array, dimension (4*N)
115: *        Workspace. Only referenced if one of NCVT, NRU, or NCC is
116: *        nonzero, and if N is at least 2.
117: *
118: *  INFO  (output) INTEGER
119: *        On exit, a value of 0 indicates a successful exit.
120: *        If INFO < 0, argument number -INFO is illegal.
121: *        If INFO > 0, the algorithm did not converge, and INFO
122: *        specifies how many superdiagonals did not converge.
123: *
124: *  Further Details
125: *  ===============
126: *
127: *  Based on contributions by
128: *     Ming Gu and Huan Ren, Computer Science Division, University of
129: *     California at Berkeley, USA
130: *
131: *  =====================================================================
132: *
133: *     .. Parameters ..
134:       DOUBLE PRECISION   ZERO
135:       PARAMETER          ( ZERO = 0.0D+0 )
136: *     ..
137: *     .. Local Scalars ..
138:       LOGICAL            ROTATE
139:       INTEGER            I, ISUB, IUPLO, J, NP1, SQRE1
140:       DOUBLE PRECISION   CS, R, SMIN, SN
141: *     ..
142: *     .. External Subroutines ..
143:       EXTERNAL           DBDSQR, DLARTG, DLASR, DSWAP, XERBLA
144: *     ..
145: *     .. External Functions ..
146:       LOGICAL            LSAME
147:       EXTERNAL           LSAME
148: *     ..
149: *     .. Intrinsic Functions ..
150:       INTRINSIC          MAX
151: *     ..
152: *     .. Executable Statements ..
153: *
154: *     Test the input parameters.
155: *
156:       INFO = 0
157:       IUPLO = 0
158:       IF( LSAME( UPLO, 'U' ) )
159:      $   IUPLO = 1
160:       IF( LSAME( UPLO, 'L' ) )
161:      $   IUPLO = 2
162:       IF( IUPLO.EQ.0 ) THEN
163:          INFO = -1
164:       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
165:          INFO = -2
166:       ELSE IF( N.LT.0 ) THEN
167:          INFO = -3
168:       ELSE IF( NCVT.LT.0 ) THEN
169:          INFO = -4
170:       ELSE IF( NRU.LT.0 ) THEN
171:          INFO = -5
172:       ELSE IF( NCC.LT.0 ) THEN
173:          INFO = -6
174:       ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
175:      $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
176:          INFO = -10
177:       ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
178:          INFO = -12
179:       ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
180:      $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
181:          INFO = -14
182:       END IF
183:       IF( INFO.NE.0 ) THEN
184:          CALL XERBLA( 'DLASDQ', -INFO )
185:          RETURN
186:       END IF
187:       IF( N.EQ.0 )
188:      $   RETURN
189: *
190: *     ROTATE is true if any singular vectors desired, false otherwise
191: *
192:       ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
193:       NP1 = N + 1
194:       SQRE1 = SQRE
195: *
196: *     If matrix non-square upper bidiagonal, rotate to be lower
197: *     bidiagonal.  The rotations are on the right.
198: *
199:       IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN
200:          DO 10 I = 1, N - 1
201:             CALL DLARTG( D( I ), E( I ), CS, SN, R )
202:             D( I ) = R
203:             E( I ) = SN*D( I+1 )
204:             D( I+1 ) = CS*D( I+1 )
205:             IF( ROTATE ) THEN
206:                WORK( I ) = CS
207:                WORK( N+I ) = SN
208:             END IF
209:    10    CONTINUE
210:          CALL DLARTG( D( N ), E( N ), CS, SN, R )
211:          D( N ) = R
212:          E( N ) = ZERO
213:          IF( ROTATE ) THEN
214:             WORK( N ) = CS
215:             WORK( N+N ) = SN
216:          END IF
217:          IUPLO = 2
218:          SQRE1 = 0
219: *
220: *        Update singular vectors if desired.
221: *
222:          IF( NCVT.GT.0 )
223:      $      CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ),
224:      $                  WORK( NP1 ), VT, LDVT )
225:       END IF
226: *
227: *     If matrix lower bidiagonal, rotate to be upper bidiagonal
228: *     by applying Givens rotations on the left.
229: *
230:       IF( IUPLO.EQ.2 ) THEN
231:          DO 20 I = 1, N - 1
232:             CALL DLARTG( D( I ), E( I ), CS, SN, R )
233:             D( I ) = R
234:             E( I ) = SN*D( I+1 )
235:             D( I+1 ) = CS*D( I+1 )
236:             IF( ROTATE ) THEN
237:                WORK( I ) = CS
238:                WORK( N+I ) = SN
239:             END IF
240:    20    CONTINUE
241: *
242: *        If matrix (N+1)-by-N lower bidiagonal, one additional
243: *        rotation is needed.
244: *
245:          IF( SQRE1.EQ.1 ) THEN
246:             CALL DLARTG( D( N ), E( N ), CS, SN, R )
247:             D( N ) = R
248:             IF( ROTATE ) THEN
249:                WORK( N ) = CS
250:                WORK( N+N ) = SN
251:             END IF
252:          END IF
253: *
254: *        Update singular vectors if desired.
255: *
256:          IF( NRU.GT.0 ) THEN
257:             IF( SQRE1.EQ.0 ) THEN
258:                CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),
259:      $                     WORK( NP1 ), U, LDU )
260:             ELSE
261:                CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ),
262:      $                     WORK( NP1 ), U, LDU )
263:             END IF
264:          END IF
265:          IF( NCC.GT.0 ) THEN
266:             IF( SQRE1.EQ.0 ) THEN
267:                CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),
268:      $                     WORK( NP1 ), C, LDC )
269:             ELSE
270:                CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),
271:      $                     WORK( NP1 ), C, LDC )
272:             END IF
273:          END IF
274:       END IF
275: *
276: *     Call DBDSQR to compute the SVD of the reduced real
277: *     N-by-N upper bidiagonal matrix.
278: *
279:       CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C,
280:      $             LDC, WORK, INFO )
281: *
282: *     Sort the singular values into ascending order (insertion sort on
283: *     singular values, but only one transposition per singular vector)
284: *
285:       DO 40 I = 1, N
286: *
287: *        Scan for smallest D(I).
288: *
289:          ISUB = I
290:          SMIN = D( I )
291:          DO 30 J = I + 1, N
292:             IF( D( J ).LT.SMIN ) THEN
293:                ISUB = J
294:                SMIN = D( J )
295:             END IF
296:    30    CONTINUE
297:          IF( ISUB.NE.I ) THEN
298: *
299: *           Swap singular values and vectors.
300: *
301:             D( ISUB ) = D( I )
302:             D( I ) = SMIN
303:             IF( NCVT.GT.0 )
304:      $         CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )
305:             IF( NRU.GT.0 )
306:      $         CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )
307:             IF( NCC.GT.0 )
308:      $         CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )
309:          END IF
310:    40 CONTINUE
311: *
312:       RETURN
313: *
314: *     End of DLASDQ
315: *
316:       END
317: