1*> \brief \b DLASD0 computes the singular values of a real upper bidiagonal n-by-m matrix B with diagonal d and off-diagonal e. Used by sbdsdc.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLASD0 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd0.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd0.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd0.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLASD0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK,
22*                          WORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       INTEGER            INFO, LDU, LDVT, N, SMLSIZ, SQRE
26*       ..
27*       .. Array Arguments ..
28*       INTEGER            IWORK( * )
29*       DOUBLE PRECISION   D( * ), E( * ), U( LDU, * ), VT( LDVT, * ),
30*      $                   WORK( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> Using a divide and conquer approach, DLASD0 computes the singular
40*> value decomposition (SVD) of a real upper bidiagonal N-by-M
41*> matrix B with diagonal D and offdiagonal E, where M = N + SQRE.
42*> The algorithm computes orthogonal matrices U and VT such that
43*> B = U * S * VT. The singular values S are overwritten on D.
44*>
45*> A related subroutine, DLASDA, computes only the singular values,
46*> and optionally, the singular vectors in compact form.
47*> \endverbatim
48*
49*  Arguments:
50*  ==========
51*
52*> \param[in] N
53*> \verbatim
54*>          N is INTEGER
55*>         On entry, the row dimension of the upper bidiagonal matrix.
56*>         This is also the dimension of the main diagonal array D.
57*> \endverbatim
58*>
59*> \param[in] SQRE
60*> \verbatim
61*>          SQRE is INTEGER
62*>         Specifies the column dimension of the bidiagonal matrix.
63*>         = 0: The bidiagonal matrix has column dimension M = N;
64*>         = 1: The bidiagonal matrix has column dimension M = N+1;
65*> \endverbatim
66*>
67*> \param[in,out] D
68*> \verbatim
69*>          D is DOUBLE PRECISION array, dimension (N)
70*>         On entry D contains the main diagonal of the bidiagonal
71*>         matrix.
72*>         On exit D, if INFO = 0, contains its singular values.
73*> \endverbatim
74*>
75*> \param[in,out] E
76*> \verbatim
77*>          E is DOUBLE PRECISION array, dimension (M-1)
78*>         Contains the subdiagonal entries of the bidiagonal matrix.
79*>         On exit, E has been destroyed.
80*> \endverbatim
81*>
82*> \param[out] U
83*> \verbatim
84*>          U is DOUBLE PRECISION array, dimension (LDU, N)
85*>         On exit, U contains the left singular vectors.
86*> \endverbatim
87*>
88*> \param[in] LDU
89*> \verbatim
90*>          LDU is INTEGER
91*>         On entry, leading dimension of U.
92*> \endverbatim
93*>
94*> \param[out] VT
95*> \verbatim
96*>          VT is DOUBLE PRECISION array, dimension (LDVT, M)
97*>         On exit, VT**T contains the right singular vectors.
98*> \endverbatim
99*>
100*> \param[in] LDVT
101*> \verbatim
102*>          LDVT is INTEGER
103*>         On entry, leading dimension of VT.
104*> \endverbatim
105*>
106*> \param[in] SMLSIZ
107*> \verbatim
108*>          SMLSIZ is INTEGER
109*>         On entry, maximum size of the subproblems at the
110*>         bottom of the computation tree.
111*> \endverbatim
112*>
113*> \param[out] IWORK
114*> \verbatim
115*>          IWORK is INTEGER array, dimension (8*N)
116*> \endverbatim
117*>
118*> \param[out] WORK
119*> \verbatim
120*>          WORK is DOUBLE PRECISION array, dimension (3*M**2+2*M)
121*> \endverbatim
122*>
123*> \param[out] INFO
124*> \verbatim
125*>          INFO is INTEGER
126*>          = 0:  successful exit.
127*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
128*>          > 0:  if INFO = 1, a singular value did not converge
129*> \endverbatim
130*
131*  Authors:
132*  ========
133*
134*> \author Univ. of Tennessee
135*> \author Univ. of California Berkeley
136*> \author Univ. of Colorado Denver
137*> \author NAG Ltd.
138*
139*> \ingroup OTHERauxiliary
140*
141*> \par Contributors:
142*  ==================
143*>
144*>     Ming Gu and Huan Ren, Computer Science Division, University of
145*>     California at Berkeley, USA
146*>
147*  =====================================================================
148      SUBROUTINE DLASD0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK,
149     $                   WORK, INFO )
150*
151*  -- LAPACK auxiliary routine --
152*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
153*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
154*
155*     .. Scalar Arguments ..
156      INTEGER            INFO, LDU, LDVT, N, SMLSIZ, SQRE
157*     ..
158*     .. Array Arguments ..
159      INTEGER            IWORK( * )
160      DOUBLE PRECISION   D( * ), E( * ), U( LDU, * ), VT( LDVT, * ),
161     $                   WORK( * )
162*     ..
163*
164*  =====================================================================
165*
166*     .. Local Scalars ..
167      INTEGER            I, I1, IC, IDXQ, IDXQC, IM1, INODE, ITEMP, IWK,
168     $                   J, LF, LL, LVL, M, NCC, ND, NDB1, NDIML, NDIMR,
169     $                   NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQREI
170      DOUBLE PRECISION   ALPHA, BETA
171*     ..
172*     .. External Subroutines ..
173      EXTERNAL           DLASD1, DLASDQ, DLASDT, XERBLA
174*     ..
175*     .. Executable Statements ..
176*
177*     Test the input parameters.
178*
179      INFO = 0
180*
181      IF( N.LT.0 ) THEN
182         INFO = -1
183      ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
184         INFO = -2
185      END IF
186*
187      M = N + SQRE
188*
189      IF( LDU.LT.N ) THEN
190         INFO = -6
191      ELSE IF( LDVT.LT.M ) THEN
192         INFO = -8
193      ELSE IF( SMLSIZ.LT.3 ) THEN
194         INFO = -9
195      END IF
196      IF( INFO.NE.0 ) THEN
197         CALL XERBLA( 'DLASD0', -INFO )
198         RETURN
199      END IF
200*
201*     If the input matrix is too small, call DLASDQ to find the SVD.
202*
203      IF( N.LE.SMLSIZ ) THEN
204         CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDVT, U, LDU, U,
205     $                LDU, WORK, INFO )
206         RETURN
207      END IF
208*
209*     Set up the computation tree.
210*
211      INODE = 1
212      NDIML = INODE + N
213      NDIMR = NDIML + N
214      IDXQ = NDIMR + N
215      IWK = IDXQ + N
216      CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
217     $             IWORK( NDIMR ), SMLSIZ )
218*
219*     For the nodes on bottom level of the tree, solve
220*     their subproblems by DLASDQ.
221*
222      NDB1 = ( ND+1 ) / 2
223      NCC = 0
224      DO 30 I = NDB1, ND
225*
226*     IC : center row of each node
227*     NL : number of rows of left  subproblem
228*     NR : number of rows of right subproblem
229*     NLF: starting row of the left   subproblem
230*     NRF: starting row of the right  subproblem
231*
232         I1 = I - 1
233         IC = IWORK( INODE+I1 )
234         NL = IWORK( NDIML+I1 )
235         NLP1 = NL + 1
236         NR = IWORK( NDIMR+I1 )
237         NRP1 = NR + 1
238         NLF = IC - NL
239         NRF = IC + 1
240         SQREI = 1
241         CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), E( NLF ),
242     $                VT( NLF, NLF ), LDVT, U( NLF, NLF ), LDU,
243     $                U( NLF, NLF ), LDU, WORK, INFO )
244         IF( INFO.NE.0 ) THEN
245            RETURN
246         END IF
247         ITEMP = IDXQ + NLF - 2
248         DO 10 J = 1, NL
249            IWORK( ITEMP+J ) = J
250   10    CONTINUE
251         IF( I.EQ.ND ) THEN
252            SQREI = SQRE
253         ELSE
254            SQREI = 1
255         END IF
256         NRP1 = NR + SQREI
257         CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), E( NRF ),
258     $                VT( NRF, NRF ), LDVT, U( NRF, NRF ), LDU,
259     $                U( NRF, NRF ), LDU, WORK, INFO )
260         IF( INFO.NE.0 ) THEN
261            RETURN
262         END IF
263         ITEMP = IDXQ + IC
264         DO 20 J = 1, NR
265            IWORK( ITEMP+J-1 ) = J
266   20    CONTINUE
267   30 CONTINUE
268*
269*     Now conquer each subproblem bottom-up.
270*
271      DO 50 LVL = NLVL, 1, -1
272*
273*        Find the first node LF and last node LL on the
274*        current level LVL.
275*
276         IF( LVL.EQ.1 ) THEN
277            LF = 1
278            LL = 1
279         ELSE
280            LF = 2**( LVL-1 )
281            LL = 2*LF - 1
282         END IF
283         DO 40 I = LF, LL
284            IM1 = I - 1
285            IC = IWORK( INODE+IM1 )
286            NL = IWORK( NDIML+IM1 )
287            NR = IWORK( NDIMR+IM1 )
288            NLF = IC - NL
289            IF( ( SQRE.EQ.0 ) .AND. ( I.EQ.LL ) ) THEN
290               SQREI = SQRE
291            ELSE
292               SQREI = 1
293            END IF
294            IDXQC = IDXQ + NLF - 1
295            ALPHA = D( IC )
296            BETA = E( IC )
297            CALL DLASD1( NL, NR, SQREI, D( NLF ), ALPHA, BETA,
298     $                   U( NLF, NLF ), LDU, VT( NLF, NLF ), LDVT,
299     $                   IWORK( IDXQC ), IWORK( IWK ), WORK, INFO )
300*
301*        Report the possible convergence failure.
302*
303            IF( INFO.NE.0 ) THEN
304               RETURN
305            END IF
306   40    CONTINUE
307   50 CONTINUE
308*
309      RETURN
310*
311*     End of DLASD0
312*
313      END
314