001:       SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
002:      $                   WORK, IWORK, INFO )
003: *
004: *  -- LAPACK 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          COMPQ, UPLO
011:       INTEGER            INFO, LDU, LDVT, N
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IQ( * ), IWORK( * )
015:       DOUBLE PRECISION   D( * ), E( * ), Q( * ), U( LDU, * ),
016:      $                   VT( LDVT, * ), WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  DBDSDC computes the singular value decomposition (SVD) of a real
023: *  N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,
024: *  using a divide and conquer method, where S is a diagonal matrix
025: *  with non-negative diagonal elements (the singular values of B), and
026: *  U and VT are orthogonal matrices of left and right singular vectors,
027: *  respectively. DBDSDC can be used to compute all singular values,
028: *  and optionally, singular vectors or singular vectors in compact form.
029: *
030: *  This code makes very mild assumptions about floating point
031: *  arithmetic. It will work on machines with a guard digit in
032: *  add/subtract, or on those binary machines without guard digits
033: *  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
034: *  It could conceivably fail on hexadecimal or decimal machines
035: *  without guard digits, but we know of none.  See DLASD3 for details.
036: *
037: *  The code currently calls DLASDQ if singular values only are desired.
038: *  However, it can be slightly modified to compute singular values
039: *  using the divide and conquer method.
040: *
041: *  Arguments
042: *  =========
043: *
044: *  UPLO    (input) CHARACTER*1
045: *          = 'U':  B is upper bidiagonal.
046: *          = 'L':  B is lower bidiagonal.
047: *
048: *  COMPQ   (input) CHARACTER*1
049: *          Specifies whether singular vectors are to be computed
050: *          as follows:
051: *          = 'N':  Compute singular values only;
052: *          = 'P':  Compute singular values and compute singular
053: *                  vectors in compact form;
054: *          = 'I':  Compute singular values and singular vectors.
055: *
056: *  N       (input) INTEGER
057: *          The order of the matrix B.  N >= 0.
058: *
059: *  D       (input/output) DOUBLE PRECISION array, dimension (N)
060: *          On entry, the n diagonal elements of the bidiagonal matrix B.
061: *          On exit, if INFO=0, the singular values of B.
062: *
063: *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
064: *          On entry, the elements of E contain the offdiagonal
065: *          elements of the bidiagonal matrix whose SVD is desired.
066: *          On exit, E has been destroyed.
067: *
068: *  U       (output) DOUBLE PRECISION array, dimension (LDU,N)
069: *          If  COMPQ = 'I', then:
070: *             On exit, if INFO = 0, U contains the left singular vectors
071: *             of the bidiagonal matrix.
072: *          For other values of COMPQ, U is not referenced.
073: *
074: *  LDU     (input) INTEGER
075: *          The leading dimension of the array U.  LDU >= 1.
076: *          If singular vectors are desired, then LDU >= max( 1, N ).
077: *
078: *  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)
079: *          If  COMPQ = 'I', then:
080: *             On exit, if INFO = 0, VT' contains the right singular
081: *             vectors of the bidiagonal matrix.
082: *          For other values of COMPQ, VT is not referenced.
083: *
084: *  LDVT    (input) INTEGER
085: *          The leading dimension of the array VT.  LDVT >= 1.
086: *          If singular vectors are desired, then LDVT >= max( 1, N ).
087: *
088: *  Q       (output) DOUBLE PRECISION array, dimension (LDQ)
089: *          If  COMPQ = 'P', then:
090: *             On exit, if INFO = 0, Q and IQ contain the left
091: *             and right singular vectors in a compact form,
092: *             requiring O(N log N) space instead of 2*N**2.
093: *             In particular, Q contains all the DOUBLE PRECISION data in
094: *             LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
095: *             words of memory, where SMLSIZ is returned by ILAENV and
096: *             is equal to the maximum size of the subproblems at the
097: *             bottom of the computation tree (usually about 25).
098: *          For other values of COMPQ, Q is not referenced.
099: *
100: *  IQ      (output) INTEGER array, dimension (LDIQ)
101: *          If  COMPQ = 'P', then:
102: *             On exit, if INFO = 0, Q and IQ contain the left
103: *             and right singular vectors in a compact form,
104: *             requiring O(N log N) space instead of 2*N**2.
105: *             In particular, IQ contains all INTEGER data in
106: *             LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
107: *             words of memory, where SMLSIZ is returned by ILAENV and
108: *             is equal to the maximum size of the subproblems at the
109: *             bottom of the computation tree (usually about 25).
110: *          For other values of COMPQ, IQ is not referenced.
111: *
112: *  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
113: *          If COMPQ = 'N' then LWORK >= (4 * N).
114: *          If COMPQ = 'P' then LWORK >= (6 * N).
115: *          If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
116: *
117: *  IWORK   (workspace) INTEGER array, dimension (8*N)
118: *
119: *  INFO    (output) INTEGER
120: *          = 0:  successful exit.
121: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
122: *          > 0:  The algorithm failed to compute an singular value.
123: *                The update process of divide and conquer failed.
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: *  Changed dimension statement in comment describing E from (N) to
134: *  (N-1).  Sven, 17 Feb 05.
135: *  =====================================================================
136: *
137: *     .. Parameters ..
138:       DOUBLE PRECISION   ZERO, ONE, TWO
139:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
140: *     ..
141: *     .. Local Scalars ..
142:       INTEGER            DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
143:      $                   ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
144:      $                   MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
145:      $                   SMLSZP, SQRE, START, WSTART, Z
146:       DOUBLE PRECISION   CS, EPS, ORGNRM, P, R, SN
147: *     ..
148: *     .. External Functions ..
149:       LOGICAL            LSAME
150:       INTEGER            ILAENV
151:       DOUBLE PRECISION   DLAMCH, DLANST
152:       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANST
153: *     ..
154: *     .. External Subroutines ..
155:       EXTERNAL           DCOPY, DLARTG, DLASCL, DLASD0, DLASDA, DLASDQ,
156:      $                   DLASET, DLASR, DSWAP, XERBLA
157: *     ..
158: *     .. Intrinsic Functions ..
159:       INTRINSIC          ABS, DBLE, INT, LOG, SIGN
160: *     ..
161: *     .. Executable Statements ..
162: *
163: *     Test the input parameters.
164: *
165:       INFO = 0
166: *
167:       IUPLO = 0
168:       IF( LSAME( UPLO, 'U' ) )
169:      $   IUPLO = 1
170:       IF( LSAME( UPLO, 'L' ) )
171:      $   IUPLO = 2
172:       IF( LSAME( COMPQ, 'N' ) ) THEN
173:          ICOMPQ = 0
174:       ELSE IF( LSAME( COMPQ, 'P' ) ) THEN
175:          ICOMPQ = 1
176:       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
177:          ICOMPQ = 2
178:       ELSE
179:          ICOMPQ = -1
180:       END IF
181:       IF( IUPLO.EQ.0 ) THEN
182:          INFO = -1
183:       ELSE IF( ICOMPQ.LT.0 ) THEN
184:          INFO = -2
185:       ELSE IF( N.LT.0 ) THEN
186:          INFO = -3
187:       ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
188:      $         N ) ) ) THEN
189:          INFO = -7
190:       ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
191:      $         N ) ) ) THEN
192:          INFO = -9
193:       END IF
194:       IF( INFO.NE.0 ) THEN
195:          CALL XERBLA( 'DBDSDC', -INFO )
196:          RETURN
197:       END IF
198: *
199: *     Quick return if possible
200: *
201:       IF( N.EQ.0 )
202:      $   RETURN
203:       SMLSIZ = ILAENV( 9, 'DBDSDC', ' ', 0, 0, 0, 0 )
204:       IF( N.EQ.1 ) THEN
205:          IF( ICOMPQ.EQ.1 ) THEN
206:             Q( 1 ) = SIGN( ONE, D( 1 ) )
207:             Q( 1+SMLSIZ*N ) = ONE
208:          ELSE IF( ICOMPQ.EQ.2 ) THEN
209:             U( 1, 1 ) = SIGN( ONE, D( 1 ) )
210:             VT( 1, 1 ) = ONE
211:          END IF
212:          D( 1 ) = ABS( D( 1 ) )
213:          RETURN
214:       END IF
215:       NM1 = N - 1
216: *
217: *     If matrix lower bidiagonal, rotate to be upper bidiagonal
218: *     by applying Givens rotations on the left
219: *
220:       WSTART = 1
221:       QSTART = 3
222:       IF( ICOMPQ.EQ.1 ) THEN
223:          CALL DCOPY( N, D, 1, Q( 1 ), 1 )
224:          CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 )
225:       END IF
226:       IF( IUPLO.EQ.2 ) THEN
227:          QSTART = 5
228:          WSTART = 2*N - 1
229:          DO 10 I = 1, N - 1
230:             CALL DLARTG( D( I ), E( I ), CS, SN, R )
231:             D( I ) = R
232:             E( I ) = SN*D( I+1 )
233:             D( I+1 ) = CS*D( I+1 )
234:             IF( ICOMPQ.EQ.1 ) THEN
235:                Q( I+2*N ) = CS
236:                Q( I+3*N ) = SN
237:             ELSE IF( ICOMPQ.EQ.2 ) THEN
238:                WORK( I ) = CS
239:                WORK( NM1+I ) = -SN
240:             END IF
241:    10    CONTINUE
242:       END IF
243: *
244: *     If ICOMPQ = 0, use DLASDQ to compute the singular values.
245: *
246:       IF( ICOMPQ.EQ.0 ) THEN
247:          CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
248:      $                LDU, WORK( WSTART ), INFO )
249:          GO TO 40
250:       END IF
251: *
252: *     If N is smaller than the minimum divide size SMLSIZ, then solve
253: *     the problem with another solver.
254: *
255:       IF( N.LE.SMLSIZ ) THEN
256:          IF( ICOMPQ.EQ.2 ) THEN
257:             CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
258:             CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
259:             CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
260:      $                   LDU, WORK( WSTART ), INFO )
261:          ELSE IF( ICOMPQ.EQ.1 ) THEN
262:             IU = 1
263:             IVT = IU + N
264:             CALL DLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
265:      $                   N )
266:             CALL DLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
267:      $                   N )
268:             CALL DLASDQ( 'U', 0, N, N, N, 0, D, E,
269:      $                   Q( IVT+( QSTART-1 )*N ), N,
270:      $                   Q( IU+( QSTART-1 )*N ), N,
271:      $                   Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
272:      $                   INFO )
273:          END IF
274:          GO TO 40
275:       END IF
276: *
277:       IF( ICOMPQ.EQ.2 ) THEN
278:          CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
279:          CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
280:       END IF
281: *
282: *     Scale.
283: *
284:       ORGNRM = DLANST( 'M', N, D, E )
285:       IF( ORGNRM.EQ.ZERO )
286:      $   RETURN
287:       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
288:       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
289: *
290:       EPS = DLAMCH( 'Epsilon' )
291: *
292:       MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
293:       SMLSZP = SMLSIZ + 1
294: *
295:       IF( ICOMPQ.EQ.1 ) THEN
296:          IU = 1
297:          IVT = 1 + SMLSIZ
298:          DIFL = IVT + SMLSZP
299:          DIFR = DIFL + MLVL
300:          Z = DIFR + MLVL*2
301:          IC = Z + MLVL
302:          IS = IC + 1
303:          POLES = IS + 1
304:          GIVNUM = POLES + 2*MLVL
305: *
306:          K = 1
307:          GIVPTR = 2
308:          PERM = 3
309:          GIVCOL = PERM + MLVL
310:       END IF
311: *
312:       DO 20 I = 1, N
313:          IF( ABS( D( I ) ).LT.EPS ) THEN
314:             D( I ) = SIGN( EPS, D( I ) )
315:          END IF
316:    20 CONTINUE
317: *
318:       START = 1
319:       SQRE = 0
320: *
321:       DO 30 I = 1, NM1
322:          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
323: *
324: *        Subproblem found. First determine its size and then
325: *        apply divide and conquer on it.
326: *
327:             IF( I.LT.NM1 ) THEN
328: *
329: *        A subproblem with E(I) small for I < NM1.
330: *
331:                NSIZE = I - START + 1
332:             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
333: *
334: *        A subproblem with E(NM1) not too small but I = NM1.
335: *
336:                NSIZE = N - START + 1
337:             ELSE
338: *
339: *        A subproblem with E(NM1) small. This implies an
340: *        1-by-1 subproblem at D(N). Solve this 1-by-1 problem
341: *        first.
342: *
343:                NSIZE = I - START + 1
344:                IF( ICOMPQ.EQ.2 ) THEN
345:                   U( N, N ) = SIGN( ONE, D( N ) )
346:                   VT( N, N ) = ONE
347:                ELSE IF( ICOMPQ.EQ.1 ) THEN
348:                   Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
349:                   Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
350:                END IF
351:                D( N ) = ABS( D( N ) )
352:             END IF
353:             IF( ICOMPQ.EQ.2 ) THEN
354:                CALL DLASD0( NSIZE, SQRE, D( START ), E( START ),
355:      $                      U( START, START ), LDU, VT( START, START ),
356:      $                      LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
357:             ELSE
358:                CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
359:      $                      E( START ), Q( START+( IU+QSTART-2 )*N ), N,
360:      $                      Q( START+( IVT+QSTART-2 )*N ),
361:      $                      IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
362:      $                      N ), Q( START+( DIFR+QSTART-2 )*N ),
363:      $                      Q( START+( Z+QSTART-2 )*N ),
364:      $                      Q( START+( POLES+QSTART-2 )*N ),
365:      $                      IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
366:      $                      N, IQ( START+PERM*N ),
367:      $                      Q( START+( GIVNUM+QSTART-2 )*N ),
368:      $                      Q( START+( IC+QSTART-2 )*N ),
369:      $                      Q( START+( IS+QSTART-2 )*N ),
370:      $                      WORK( WSTART ), IWORK, INFO )
371:                IF( INFO.NE.0 ) THEN
372:                   RETURN
373:                END IF
374:             END IF
375:             START = I + 1
376:          END IF
377:    30 CONTINUE
378: *
379: *     Unscale
380: *
381:       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
382:    40 CONTINUE
383: *
384: *     Use Selection Sort to minimize swaps of singular vectors
385: *
386:       DO 60 II = 2, N
387:          I = II - 1
388:          KK = I
389:          P = D( I )
390:          DO 50 J = II, N
391:             IF( D( J ).GT.P ) THEN
392:                KK = J
393:                P = D( J )
394:             END IF
395:    50    CONTINUE
396:          IF( KK.NE.I ) THEN
397:             D( KK ) = D( I )
398:             D( I ) = P
399:             IF( ICOMPQ.EQ.1 ) THEN
400:                IQ( I ) = KK
401:             ELSE IF( ICOMPQ.EQ.2 ) THEN
402:                CALL DSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
403:                CALL DSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
404:             END IF
405:          ELSE IF( ICOMPQ.EQ.1 ) THEN
406:             IQ( I ) = I
407:          END IF
408:    60 CONTINUE
409: *
410: *     If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
411: *
412:       IF( ICOMPQ.EQ.1 ) THEN
413:          IF( IUPLO.EQ.1 ) THEN
414:             IQ( N ) = 1
415:          ELSE
416:             IQ( N ) = 0
417:          END IF
418:       END IF
419: *
420: *     If B is lower bidiagonal, update U by those Givens rotations
421: *     which rotated B to be upper bidiagonal
422: *
423:       IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) )
424:      $   CALL DLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU )
425: *
426:       RETURN
427: *
428: *     End of DBDSDC
429: *
430:       END
431: