001:       SUBROUTINE SLASDQ( 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: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          UPLO
011:       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               C( LDC, * ), D( * ), E( * ), U( LDU, * ),
015:      $                   VT( LDVT, * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  SLASDQ computes the singular value decomposition (SVD) of a real
022: *  (upper or lower) bidiagonal matrix with diagonal D and offdiagonal
023: *  E, accumulating the transformations if desired. Letting B denote
024: *  the input bidiagonal matrix, the algorithm computes orthogonal
025: *  matrices Q and P such that B = Q * S * P' (P' denotes the transpose
026: *  of P). The singular values S are overwritten on D.
027: *
028: *  The input matrix U  is changed to U  * Q  if desired.
029: *  The input matrix VT is changed to P' * VT if desired.
030: *  The input matrix C  is changed to Q' * C  if desired.
031: *
032: *  See "Computing  Small Singular Values of Bidiagonal Matrices With
033: *  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
034: *  LAPACK Working Note #3, for a detailed description of the algorithm.
035: *
036: *  Arguments
037: *  =========
038: *
039: *  UPLO  (input) CHARACTER*1
040: *        On entry, UPLO specifies whether the input bidiagonal matrix
041: *        is upper or lower bidiagonal, and wether it is square are
042: *        not.
043: *           UPLO = 'U' or 'u'   B is upper bidiagonal.
044: *           UPLO = 'L' or 'l'   B is lower bidiagonal.
045: *
046: *  SQRE  (input) INTEGER
047: *        = 0: then the input matrix is N-by-N.
048: *        = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and
049: *             (N+1)-by-N if UPLU = 'L'.
050: *
051: *        The bidiagonal matrix has
052: *        N = NL + NR + 1 rows and
053: *        M = N + SQRE >= N columns.
054: *
055: *  N     (input) INTEGER
056: *        On entry, N specifies the number of rows and columns
057: *        in the matrix. N must be at least 0.
058: *
059: *  NCVT  (input) INTEGER
060: *        On entry, NCVT specifies the number of columns of
061: *        the matrix VT. NCVT must be at least 0.
062: *
063: *  NRU   (input) INTEGER
064: *        On entry, NRU specifies the number of rows of
065: *        the matrix U. NRU must be at least 0.
066: *
067: *  NCC   (input) INTEGER
068: *        On entry, NCC specifies the number of columns of
069: *        the matrix C. NCC must be at least 0.
070: *
071: *  D     (input/output) REAL array, dimension (N)
072: *        On entry, D contains the diagonal entries of the
073: *        bidiagonal matrix whose SVD is desired. On normal exit,
074: *        D contains the singular values in ascending order.
075: *
076: *  E     (input/output) REAL array.
077: *        dimension is (N-1) if SQRE = 0 and N if SQRE = 1.
078: *        On entry, the entries of E contain the offdiagonal entries
079: *        of the bidiagonal matrix whose SVD is desired. On normal
080: *        exit, E will contain 0. If the algorithm does not converge,
081: *        D and E will contain the diagonal and superdiagonal entries
082: *        of a bidiagonal matrix orthogonally equivalent to the one
083: *        given as input.
084: *
085: *  VT    (input/output) REAL array, dimension (LDVT, NCVT)
086: *        On entry, contains a matrix which on exit has been
087: *        premultiplied by P', dimension N-by-NCVT if SQRE = 0
088: *        and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
089: *
090: *  LDVT  (input) INTEGER
091: *        On entry, LDVT specifies the leading dimension of VT as
092: *        declared in the calling (sub) program. LDVT must be at
093: *        least 1. If NCVT is nonzero LDVT must also be at least N.
094: *
095: *  U     (input/output) REAL array, dimension (LDU, N)
096: *        On entry, contains a  matrix which on exit has been
097: *        postmultiplied by Q, dimension NRU-by-N if SQRE = 0
098: *        and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
099: *
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) REAL 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) REAL 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:       REAL               ZERO
136:       PARAMETER          ( ZERO = 0.0E+0 )
137: *     ..
138: *     .. Local Scalars ..
139:       LOGICAL            ROTATE
140:       INTEGER            I, ISUB, IUPLO, J, NP1, SQRE1
141:       REAL               CS, R, SMIN, SN
142: *     ..
143: *     .. External Subroutines ..
144:       EXTERNAL           SBDSQR, SLARTG, SLASR, SSWAP, 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( 'SLASDQ', -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 SLARTG( 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 SLARTG( 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 SLASR( '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 SLARTG( 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 SLARTG( 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 SLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),
260:      $                     WORK( NP1 ), U, LDU )
261:             ELSE
262:                CALL SLASR( '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 SLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),
269:      $                     WORK( NP1 ), C, LDC )
270:             ELSE
271:                CALL SLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),
272:      $                     WORK( NP1 ), C, LDC )
273:             END IF
274:          END IF
275:       END IF
276: *
277: *     Call SBDSQR to compute the SVD of the reduced real
278: *     N-by-N upper bidiagonal matrix.
279: *
280:       CALL SBDSQR( '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 SSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )
306:             IF( NRU.GT.0 )
307:      $         CALL SSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )
308:             IF( NCC.GT.0 )
309:      $         CALL SSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )
310:          END IF
311:    40 CONTINUE
312: *
313:       RETURN
314: *
315: *     End of SLASDQ
316: *
317:       END
318: