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