1      SUBROUTINE SBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
2     $                   WORK, IWORK, INFO )
3*
4*  -- LAPACK 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          COMPQ, UPLO
11      INTEGER            INFO, LDU, LDVT, N
12*     ..
13*     .. Array Arguments ..
14      INTEGER            IQ( * ), IWORK( * )
15      REAL               D( * ), E( * ), Q( * ), U( LDU, * ),
16     $                   VT( LDVT, * ), WORK( * )
17*     ..
18*     .. Common blocks ..
19      COMMON             / LATIME / OPS, ITCNT
20*     ..
21*     .. Scalars in Common ..
22      REAL               ITCNT, OPS
23*     ..
24*
25*  Purpose
26*  =======
27*
28*  SBDSDC computes the singular value decomposition (SVD) of a real
29*  N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,
30*  using a divide and conquer method, where S is a diagonal matrix
31*  with non-negative diagonal elements (the singular values of B), and
32*  U and VT are orthogonal matrices of left and right singular vectors,
33*  respectively. SBDSDC can be used to compute all singular values,
34*  and optionally, singular vectors or singular vectors in compact form.
35*
36*  This code makes very mild assumptions about floating point
37*  arithmetic. It will work on machines with a guard digit in
38*  add/subtract, or on those binary machines without guard digits
39*  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
40*  It could conceivably fail on hexadecimal or decimal machines
41*  without guard digits, but we know of none.  See SLASD3 for details.
42*
43*  The code currently call SLASDQ if singular values only are desired.
44*  However, it can be slightly modified to compute singular values
45*  using the divide and conquer method.
46*
47*  Arguments
48*  =========
49*
50*  UPLO    (input) CHARACTER*1
51*          = 'U':  B is upper bidiagonal.
52*          = 'L':  B is lower bidiagonal.
53*
54*  COMPQ   (input) CHARACTER*1
55*          Specifies whether singular vectors are to be computed
56*          as follows:
57*          = 'N':  Compute singular values only;
58*          = 'P':  Compute singular values and compute singular
59*                  vectors in compact form;
60*          = 'I':  Compute singular values and singular vectors.
61*
62*  N       (input) INTEGER
63*          The order of the matrix B.  N >= 0.
64*
65*  D       (input/output) REAL array, dimension (N)
66*          On entry, the n diagonal elements of the bidiagonal matrix B.
67*          On exit, if INFO=0, the singular values of B.
68*
69*  E       (input/output) REAL array, dimension (N)
70*          On entry, the elements of E contain the offdiagonal
71*          elements of the bidiagonal matrix whose SVD is desired.
72*          On exit, E has been destroyed.
73*
74*  U       (output) REAL array, dimension (LDU,N)
75*          If  COMPQ = 'I', then:
76*             On exit, if INFO = 0, U contains the left singular vectors
77*             of the bidiagonal matrix.
78*          For other values of COMPQ, U is not referenced.
79*
80*  LDU     (input) INTEGER
81*          The leading dimension of the array U.  LDU >= 1.
82*          If singular vectors are desired, then LDU >= max( 1, N ).
83*
84*  VT      (output) REAL array, dimension (LDVT,N)
85*          If  COMPQ = 'I', then:
86*             On exit, if INFO = 0, VT' contains the right singular
87*             vectors of the bidiagonal matrix.
88*          For other values of COMPQ, VT is not referenced.
89*
90*  LDVT    (input) INTEGER
91*          The leading dimension of the array VT.  LDVT >= 1.
92*          If singular vectors are desired, then LDVT >= max( 1, N ).
93*
94*  Q       (output) REAL array, dimension (LDQ)
95*          If  COMPQ = 'P', then:
96*             On exit, if INFO = 0, Q and IQ contain the left
97*             and right singular vectors in a compact form,
98*             requiring O(N log N) space instead of 2*N**2.
99*             In particular, Q contains all the REAL data in
100*             LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
101*             words of memory, where SMLSIZ is returned by ILAENV and
102*             is equal to the maximum size of the subproblems at the
103*             bottom of the computation tree (usually about 25).
104*          For other values of COMPQ, Q is not referenced.
105*
106*  IQ      (output) INTEGER array, dimension (LDIQ)
107*          If  COMPQ = 'P', then:
108*             On exit, if INFO = 0, Q and IQ contain the left
109*             and right singular vectors in a compact form,
110*             requiring O(N log N) space instead of 2*N**2.
111*             In particular, IQ contains all INTEGER data in
112*             LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
113*             words of memory, where SMLSIZ is returned by ILAENV and
114*             is equal to the maximum size of the subproblems at the
115*             bottom of the computation tree (usually about 25).
116*          For other values of COMPQ, IQ is not referenced.
117*
118*  WORK    (workspace) REAL array, dimension (LWORK)
119*          If COMPQ = 'N' then LWORK >= (4 * N).
120*          If COMPQ = 'P' then LWORK >= (6 * N).
121*          If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
122*
123*  IWORK   (workspace) INTEGER array, dimension (7*N)
124*
125*  INFO    (output) INTEGER
126*          = 0:  successful exit.
127*          < 0:  if INFO = -i, the i-th argument had an illegal value.
128*          > 0:  The algorithm failed to compute an singular value.
129*                The update process of divide and conquer failed.
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      REAL               ZERO, ONE, TWO
142      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
143*     ..
144*     .. Local Scalars ..
145      INTEGER            DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
146     $                   ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
147     $                   MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
148     $                   SMLSZP, SQRE, START, WSTART, Z
149      REAL               CS, EPS, ORGNRM, P, R, SN
150*     ..
151*     .. External Functions ..
152      LOGICAL            LSAME
153      INTEGER            ILAENV
154      REAL               SLAMCH, SLANST
155      EXTERNAL           SLAMCH, SLANST, ILAENV, LSAME
156*     ..
157*     .. External Subroutines ..
158      EXTERNAL           SCOPY, SLARTG, SLASCL, SLASD0, SLASDA, SLASDQ,
159     $                   SLASET, SLASR, SSWAP, XERBLA
160*     ..
161*     .. Intrinsic Functions ..
162      INTRINSIC          REAL, ABS, INT, LOG, SIGN
163*     ..
164*     .. Executable Statements ..
165*
166*     Test the input parameters.
167*
168      INFO = 0
169*
170      IUPLO = 0
171      IF( LSAME( UPLO, 'U' ) )
172     $   IUPLO = 1
173      IF( LSAME( UPLO, 'L' ) )
174     $   IUPLO = 2
175      IF( LSAME( COMPQ, 'N' ) ) THEN
176         ICOMPQ = 0
177      ELSE IF( LSAME( COMPQ, 'P' ) ) THEN
178         ICOMPQ = 1
179      ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
180         ICOMPQ = 2
181      ELSE
182         ICOMPQ = -1
183      END IF
184      IF( IUPLO.EQ.0 ) THEN
185         INFO = -1
186      ELSE IF( ICOMPQ.LT.0 ) THEN
187         INFO = -2
188      ELSE IF( N.LT.0 ) THEN
189         INFO = -3
190      ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
191     $         N ) ) ) THEN
192         INFO = -7
193      ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
194     $         N ) ) ) THEN
195         INFO = -9
196      END IF
197      IF( INFO.NE.0 ) THEN
198         CALL XERBLA( 'SBDSDC', -INFO )
199         RETURN
200      END IF
201*
202*     Quick return if possible
203*
204      IF( N.EQ.0 )
205     $   RETURN
206      SMLSIZ = ILAENV( 9, 'SBDSDC', ' ', 0, 0, 0, 0 )
207      IF( N.EQ.1 ) THEN
208         IF( ICOMPQ.EQ.1 ) THEN
209            Q( 1 ) = SIGN( ONE, D( 1 ) )
210            Q( 1+SMLSIZ*N ) = ONE
211         ELSE IF( ICOMPQ.EQ.2 ) THEN
212            U( 1, 1 ) = SIGN( ONE, D( 1 ) )
213            VT( 1, 1 ) = ONE
214         END IF
215         D( 1 ) = ABS( D( 1 ) )
216         RETURN
217      END IF
218      NM1 = N - 1
219*
220*     If matrix lower bidiagonal, rotate to be upper bidiagonal
221*     by applying Givens rotations on the left
222*
223      WSTART = 1
224      QSTART = 3
225      IF( ICOMPQ.EQ.1 ) THEN
226         CALL SCOPY( N, D, 1, Q( 1 ), 1 )
227         CALL SCOPY( N-1, E, 1, Q( N+1 ), 1 )
228      END IF
229      IF( IUPLO.EQ.2 ) THEN
230         QSTART = 5
231         WSTART = 2*N - 1
232         OPS = OPS + REAL( 8*( N-1 ) )
233         DO 10 I = 1, N - 1
234            CALL SLARTG( D( I ), E( I ), CS, SN, R )
235            D( I ) = R
236            E( I ) = SN*D( I+1 )
237            D( I+1 ) = CS*D( I+1 )
238            IF( ICOMPQ.EQ.1 ) THEN
239               Q( I+2*N ) = CS
240               Q( I+3*N ) = SN
241            ELSE IF( ICOMPQ.EQ.2 ) THEN
242               WORK( I ) = CS
243               WORK( NM1+I ) = -SN
244            END IF
245   10    CONTINUE
246      END IF
247*
248*     If ICOMPQ = 0, use SLASDQ to compute the singular values.
249*
250      IF( ICOMPQ.EQ.0 ) THEN
251         CALL SLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
252     $                LDU, WORK( WSTART ), INFO )
253         GO TO 40
254      END IF
255*
256*     If N is smaller than the minimum divide size SMLSIZ, then solve
257*     the problem with another solver.
258*
259      IF( N.LE.SMLSIZ ) THEN
260         IF( ICOMPQ.EQ.2 ) THEN
261            CALL SLASET( 'A', N, N, ZERO, ONE, U, LDU )
262            CALL SLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
263            CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
264     $                   LDU, WORK( WSTART ), INFO )
265         ELSE IF( ICOMPQ.EQ.1 ) THEN
266            IU = 1
267            IVT = IU + N
268            CALL SLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
269     $                   N )
270            CALL SLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
271     $                   N )
272            CALL SLASDQ( 'U', 0, N, N, N, 0, D, E,
273     $                   Q( IVT+( QSTART-1 )*N ), N,
274     $                   Q( IU+( QSTART-1 )*N ), N,
275     $                   Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
276     $                   INFO )
277         END IF
278         GO TO 40
279      END IF
280*
281      IF( ICOMPQ.EQ.2 ) THEN
282         CALL SLASET( 'A', N, N, ZERO, ONE, U, LDU )
283         CALL SLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
284      END IF
285*
286*     Scale.
287*
288      ORGNRM = SLANST( 'M', N, D, E )
289      IF( ORGNRM.EQ.ZERO )
290     $   RETURN
291      OPS = OPS + REAL( N + NM1 )
292      CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
293      CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
294*
295      EPS = SLAMCH( 'Epsilon' )
296*
297      MLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
298      SMLSZP = SMLSIZ + 1
299*
300      IF( ICOMPQ.EQ.1 ) THEN
301         IU = 1
302         IVT = 1 + SMLSIZ
303         DIFL = IVT + SMLSZP
304         DIFR = DIFL + MLVL
305         Z = DIFR + MLVL*2
306         IC = Z + MLVL
307         IS = IC + 1
308         POLES = IS + 1
309         GIVNUM = POLES + 2*MLVL
310*
311         K = 1
312         GIVPTR = 2
313         PERM = 3
314         GIVCOL = PERM + MLVL
315      END IF
316*
317      DO 20 I = 1, N
318         IF( ABS( D( I ) ).LT.EPS ) THEN
319            D( I ) = SIGN( EPS, D( I ) )
320         END IF
321   20 CONTINUE
322*
323      START = 1
324      SQRE = 0
325*
326      DO 30 I = 1, NM1
327         IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
328*
329*        Subproblem found. First determine its size and then
330*        apply divide and conquer on it.
331*
332            IF( I.LT.NM1 ) THEN
333*
334*        A subproblem with E(I) small for I < NM1.
335*
336               NSIZE = I - START + 1
337            ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
338*
339*        A subproblem with E(NM1) not too small but I = NM1.
340*
341               NSIZE = N - START + 1
342            ELSE
343*
344*        A subproblem with E(NM1) small. This implies an
345*        1-by-1 subproblem at D(N). Solve this 1-by-1 problem
346*        first.
347*
348               NSIZE = I - START + 1
349               IF( ICOMPQ.EQ.2 ) THEN
350                  U( N, N ) = SIGN( ONE, D( N ) )
351                  VT( N, N ) = ONE
352               ELSE IF( ICOMPQ.EQ.1 ) THEN
353                  Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
354                  Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
355               END IF
356               D( N ) = ABS( D( N ) )
357            END IF
358            IF( ICOMPQ.EQ.2 ) THEN
359               CALL SLASD0( NSIZE, SQRE, D( START ), E( START ),
360     $                      U( START, START ), LDU, VT( START, START ),
361     $                      LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
362            ELSE
363               CALL SLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
364     $                      E( START ), Q( START+( IU+QSTART-2 )*N ), N,
365     $                      Q( START+( IVT+QSTART-2 )*N ),
366     $                      IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
367     $                      N ), Q( START+( DIFR+QSTART-2 )*N ),
368     $                      Q( START+( Z+QSTART-2 )*N ),
369     $                      Q( START+( POLES+QSTART-2 )*N ),
370     $                      IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
371     $                      N, IQ( START+PERM*N ),
372     $                      Q( START+( GIVNUM+QSTART-2 )*N ),
373     $                      Q( START+( IC+QSTART-2 )*N ),
374     $                      Q( START+( IS+QSTART-2 )*N ),
375     $                      WORK( WSTART ), IWORK, INFO )
376               IF( INFO.NE.0 ) THEN
377                  RETURN
378               END IF
379            END IF
380            START = I + 1
381         END IF
382   30 CONTINUE
383*
384*     Unscale
385*
386      OPS = OPS + REAL( N )
387      CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
388   40 CONTINUE
389*
390*     Use Selection Sort to minimize swaps of singular vectors
391*
392      DO 60 II = 2, N
393         I = II - 1
394         KK = I
395         P = D( I )
396         DO 50 J = II, N
397            IF( D( J ).GT.P ) THEN
398               KK = J
399               P = D( J )
400            END IF
401   50    CONTINUE
402         IF( KK.NE.I ) THEN
403            D( KK ) = D( I )
404            D( I ) = P
405            IF( ICOMPQ.EQ.1 ) THEN
406               IQ( I ) = KK
407            ELSE IF( ICOMPQ.EQ.2 ) THEN
408               CALL SSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
409               CALL SSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
410            END IF
411         ELSE IF( ICOMPQ.EQ.1 ) THEN
412            IQ( I ) = I
413         END IF
414   60 CONTINUE
415*
416*     If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
417*
418      IF( ICOMPQ.EQ.1 ) THEN
419         IF( IUPLO.EQ.1 ) THEN
420            IQ( N ) = 1
421         ELSE
422            IQ( N ) = 0
423         END IF
424      END IF
425*
426*     If B is lower bidiagonal, update U by those Givens rotations
427*     which rotated B to be upper bidiagonal
428*
429      IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) ) THEN
430         OPS = OPS + REAL( 6*( N-1 )*N )
431         CALL SLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU )
432      END IF
433*
434      RETURN
435*
436*     End of SBDSDC
437*
438      END
439