1*> \brief \b SLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. 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 SLASD1 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd1.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd1.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd1.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
22*                          IDXQ, IWORK, WORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       INTEGER            INFO, LDU, LDVT, NL, NR, SQRE
26*       REAL               ALPHA, BETA
27*       ..
28*       .. Array Arguments ..
29*       INTEGER            IDXQ( * ), IWORK( * )
30*       REAL               D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> SLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
40*> where N = NL + NR + 1 and M = N + SQRE. SLASD1 is called from SLASD0.
41*>
42*> A related subroutine SLASD7 handles the case in which the singular
43*> values (and the singular vectors in factored form) are desired.
44*>
45*> SLASD1 computes the SVD as follows:
46*>
47*>               ( D1(in)    0    0       0 )
48*>   B = U(in) * (   Z1**T   a   Z2**T    b ) * VT(in)
49*>               (   0       0   D2(in)   0 )
50*>
51*>     = U(out) * ( D(out) 0) * VT(out)
52*>
53*> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
54*> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
55*> elsewhere; and the entry b is empty if SQRE = 0.
56*>
57*> The left singular vectors of the original matrix are stored in U, and
58*> the transpose of the right singular vectors are stored in VT, and the
59*> singular values are in D.  The algorithm consists of three stages:
60*>
61*>    The first stage consists of deflating the size of the problem
62*>    when there are multiple singular values or when there are zeros in
63*>    the Z vector.  For each such occurrence the dimension of the
64*>    secular equation problem is reduced by one.  This stage is
65*>    performed by the routine SLASD2.
66*>
67*>    The second stage consists of calculating the updated
68*>    singular values. This is done by finding the square roots of the
69*>    roots of the secular equation via the routine SLASD4 (as called
70*>    by SLASD3). This routine also calculates the singular vectors of
71*>    the current problem.
72*>
73*>    The final stage consists of computing the updated singular vectors
74*>    directly using the updated singular values.  The singular vectors
75*>    for the current problem are multiplied with the singular vectors
76*>    from the overall problem.
77*> \endverbatim
78*
79*  Arguments:
80*  ==========
81*
82*> \param[in] NL
83*> \verbatim
84*>          NL is INTEGER
85*>         The row dimension of the upper block.  NL >= 1.
86*> \endverbatim
87*>
88*> \param[in] NR
89*> \verbatim
90*>          NR is INTEGER
91*>         The row dimension of the lower block.  NR >= 1.
92*> \endverbatim
93*>
94*> \param[in] SQRE
95*> \verbatim
96*>          SQRE is INTEGER
97*>         = 0: the lower block is an NR-by-NR square matrix.
98*>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
99*>
100*>         The bidiagonal matrix has row dimension N = NL + NR + 1,
101*>         and column dimension M = N + SQRE.
102*> \endverbatim
103*>
104*> \param[in,out] D
105*> \verbatim
106*>          D is REAL array, dimension (NL+NR+1).
107*>         N = NL+NR+1
108*>         On entry D(1:NL,1:NL) contains the singular values of the
109*>         upper block; and D(NL+2:N) contains the singular values of
110*>         the lower block. On exit D(1:N) contains the singular values
111*>         of the modified matrix.
112*> \endverbatim
113*>
114*> \param[in,out] ALPHA
115*> \verbatim
116*>          ALPHA is REAL
117*>         Contains the diagonal element associated with the added row.
118*> \endverbatim
119*>
120*> \param[in,out] BETA
121*> \verbatim
122*>          BETA is REAL
123*>         Contains the off-diagonal element associated with the added
124*>         row.
125*> \endverbatim
126*>
127*> \param[in,out] U
128*> \verbatim
129*>          U is REAL array, dimension (LDU,N)
130*>         On entry U(1:NL, 1:NL) contains the left singular vectors of
131*>         the upper block; U(NL+2:N, NL+2:N) contains the left singular
132*>         vectors of the lower block. On exit U contains the left
133*>         singular vectors of the bidiagonal matrix.
134*> \endverbatim
135*>
136*> \param[in] LDU
137*> \verbatim
138*>          LDU is INTEGER
139*>         The leading dimension of the array U.  LDU >= max( 1, N ).
140*> \endverbatim
141*>
142*> \param[in,out] VT
143*> \verbatim
144*>          VT is REAL array, dimension (LDVT,M)
145*>         where M = N + SQRE.
146*>         On entry VT(1:NL+1, 1:NL+1)**T contains the right singular
147*>         vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains
148*>         the right singular vectors of the lower block. On exit
149*>         VT**T contains the right singular vectors of the
150*>         bidiagonal matrix.
151*> \endverbatim
152*>
153*> \param[in] LDVT
154*> \verbatim
155*>          LDVT is INTEGER
156*>         The leading dimension of the array VT.  LDVT >= max( 1, M ).
157*> \endverbatim
158*>
159*> \param[in,out] IDXQ
160*> \verbatim
161*>          IDXQ is INTEGER array, dimension (N)
162*>         This contains the permutation which will reintegrate the
163*>         subproblem just solved back into sorted order, i.e.
164*>         D( IDXQ( I = 1, N ) ) will be in ascending order.
165*> \endverbatim
166*>
167*> \param[out] IWORK
168*> \verbatim
169*>          IWORK is INTEGER array, dimension (4*N)
170*> \endverbatim
171*>
172*> \param[out] WORK
173*> \verbatim
174*>          WORK is REAL array, dimension (3*M**2+2*M)
175*> \endverbatim
176*>
177*> \param[out] INFO
178*> \verbatim
179*>          INFO is INTEGER
180*>          = 0:  successful exit.
181*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
182*>          > 0:  if INFO = 1, a singular value did not converge
183*> \endverbatim
184*
185*  Authors:
186*  ========
187*
188*> \author Univ. of Tennessee
189*> \author Univ. of California Berkeley
190*> \author Univ. of Colorado Denver
191*> \author NAG Ltd.
192*
193*> \ingroup OTHERauxiliary
194*
195*> \par Contributors:
196*  ==================
197*>
198*>     Ming Gu and Huan Ren, Computer Science Division, University of
199*>     California at Berkeley, USA
200*>
201*  =====================================================================
202      SUBROUTINE SLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
203     $                   IDXQ, IWORK, WORK, INFO )
204*
205*  -- LAPACK auxiliary routine --
206*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
207*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
208*
209*     .. Scalar Arguments ..
210      INTEGER            INFO, LDU, LDVT, NL, NR, SQRE
211      REAL               ALPHA, BETA
212*     ..
213*     .. Array Arguments ..
214      INTEGER            IDXQ( * ), IWORK( * )
215      REAL               D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
216*     ..
217*
218*  =====================================================================
219*
220*     .. Parameters ..
221*
222      REAL               ONE, ZERO
223      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
224*     ..
225*     .. Local Scalars ..
226      INTEGER            COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
227     $                   IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
228      REAL               ORGNRM
229*     ..
230*     .. External Subroutines ..
231      EXTERNAL           SLAMRG, SLASCL, SLASD2, SLASD3, XERBLA
232*     ..
233*     .. Intrinsic Functions ..
234      INTRINSIC          ABS, MAX
235*     ..
236*     .. Executable Statements ..
237*
238*     Test the input parameters.
239*
240      INFO = 0
241*
242      IF( NL.LT.1 ) THEN
243         INFO = -1
244      ELSE IF( NR.LT.1 ) THEN
245         INFO = -2
246      ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
247         INFO = -3
248      END IF
249      IF( INFO.NE.0 ) THEN
250         CALL XERBLA( 'SLASD1', -INFO )
251         RETURN
252      END IF
253*
254      N = NL + NR + 1
255      M = N + SQRE
256*
257*     The following values are for bookkeeping purposes only.  They are
258*     integer pointers which indicate the portion of the workspace
259*     used by a particular array in SLASD2 and SLASD3.
260*
261      LDU2 = N
262      LDVT2 = M
263*
264      IZ = 1
265      ISIGMA = IZ + M
266      IU2 = ISIGMA + N
267      IVT2 = IU2 + LDU2*N
268      IQ = IVT2 + LDVT2*M
269*
270      IDX = 1
271      IDXC = IDX + N
272      COLTYP = IDXC + N
273      IDXP = COLTYP + N
274*
275*     Scale.
276*
277      ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
278      D( NL+1 ) = ZERO
279      DO 10 I = 1, N
280         IF( ABS( D( I ) ).GT.ORGNRM ) THEN
281            ORGNRM = ABS( D( I ) )
282         END IF
283   10 CONTINUE
284      CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
285      ALPHA = ALPHA / ORGNRM
286      BETA = BETA / ORGNRM
287*
288*     Deflate singular values.
289*
290      CALL SLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU,
291     $             VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2,
292     $             WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ),
293     $             IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO )
294*
295*     Solve Secular Equation and update singular vectors.
296*
297      LDQ = K
298      CALL SLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ),
299     $             U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ),
300     $             LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ),
301     $             INFO )
302*
303*     Report the possible convergence failure.
304*
305      IF( INFO.NE.0 ) THEN
306         RETURN
307      END IF
308*
309*     Unscale.
310*
311      CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
312*
313*     Prepare the IDXQ sorting permutation.
314*
315      N1 = K
316      N2 = N - K
317      CALL SLAMRG( N1, N2, D, 1, -1, IDXQ )
318*
319      RETURN
320*
321*     End of SLASD1
322*
323      END
324