LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
dlasd1.f
Go to the documentation of this file.
1 *> \brief \b DLASD1 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
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd1.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd1.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd1.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLASD1( 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 * DOUBLE PRECISION ALPHA, BETA
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IDXQ( * ), IWORK( * )
30 * DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
40 *> where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
41 *>
42 *> A related subroutine DLASD7 handles the case in which the singular
43 *> values (and the singular vectors in factored form) are desired.
44 *>
45 *> DLASD1 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 occurence the dimension of the
64 *> secular equation problem is reduced by one. This stage is
65 *> performed by the routine DLASD2.
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 DLASD4 (as called
70 *> by DLASD3). 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 DOUBLE PRECISION array,
107 *> dimension (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 DOUBLE PRECISION
117 *> Contains the diagonal element associated with the added row.
118 *> \endverbatim
119 *>
120 *> \param[in,out] BETA
121 *> \verbatim
122 *> BETA is DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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[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 DOUBLE PRECISION 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 *> \date September 2012
194 *
195 *> \ingroup auxOTHERauxiliary
196 *
197 *> \par Contributors:
198 * ==================
199 *>
200 *> Ming Gu and Huan Ren, Computer Science Division, University of
201 *> California at Berkeley, USA
202 *>
203 * =====================================================================
204  SUBROUTINE dlasd1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
205  \$ idxq, iwork, work, info )
206 *
207 * -- LAPACK auxiliary routine (version 3.4.2) --
208 * -- LAPACK is a software package provided by Univ. of Tennessee, --
209 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
210 * September 2012
211 *
212 * .. Scalar Arguments ..
213  INTEGER info, ldu, ldvt, nl, nr, sqre
214  DOUBLE PRECISION alpha, beta
215 * ..
216 * .. Array Arguments ..
217  INTEGER idxq( * ), iwork( * )
218  DOUBLE PRECISION d( * ), u( ldu, * ), vt( ldvt, * ), work( * )
219 * ..
220 *
221 * =====================================================================
222 *
223 * .. Parameters ..
224 *
225  DOUBLE PRECISION one, zero
226  parameter( one = 1.0d+0, zero = 0.0d+0 )
227 * ..
228 * .. Local Scalars ..
229  INTEGER coltyp, i, idx, idxc, idxp, iq, isigma, iu2,
230  \$ ivt2, iz, k, ldq, ldu2, ldvt2, m, n, n1, n2
231  DOUBLE PRECISION orgnrm
232 * ..
233 * .. External Subroutines ..
234  EXTERNAL dlamrg, dlascl, dlasd2, dlasd3, xerbla
235 * ..
236 * .. Intrinsic Functions ..
237  INTRINSIC abs, max
238 * ..
239 * .. Executable Statements ..
240 *
241 * Test the input parameters.
242 *
243  info = 0
244 *
245  IF( nl.LT.1 ) THEN
246  info = -1
247  ELSE IF( nr.LT.1 ) THEN
248  info = -2
249  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
250  info = -3
251  END IF
252  IF( info.NE.0 ) THEN
253  CALL xerbla( 'DLASD1', -info )
254  return
255  END IF
256 *
257  n = nl + nr + 1
258  m = n + sqre
259 *
260 * The following values are for bookkeeping purposes only. They are
261 * integer pointers which indicate the portion of the workspace
262 * used by a particular array in DLASD2 and DLASD3.
263 *
264  ldu2 = n
265  ldvt2 = m
266 *
267  iz = 1
268  isigma = iz + m
269  iu2 = isigma + n
270  ivt2 = iu2 + ldu2*n
271  iq = ivt2 + ldvt2*m
272 *
273  idx = 1
274  idxc = idx + n
275  coltyp = idxc + n
276  idxp = coltyp + n
277 *
278 * Scale.
279 *
280  orgnrm = max( abs( alpha ), abs( beta ) )
281  d( nl+1 ) = zero
282  DO 10 i = 1, n
283  IF( abs( d( i ) ).GT.orgnrm ) THEN
284  orgnrm = abs( d( i ) )
285  END IF
286  10 continue
287  CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
288  alpha = alpha / orgnrm
289  beta = beta / orgnrm
290 *
291 * Deflate singular values.
292 *
293  CALL dlasd2( nl, nr, sqre, k, d, work( iz ), alpha, beta, u, ldu,
294  \$ vt, ldvt, work( isigma ), work( iu2 ), ldu2,
295  \$ work( ivt2 ), ldvt2, iwork( idxp ), iwork( idx ),
296  \$ iwork( idxc ), idxq, iwork( coltyp ), info )
297 *
298 * Solve Secular Equation and update singular vectors.
299 *
300  ldq = k
301  CALL dlasd3( nl, nr, sqre, k, d, work( iq ), ldq, work( isigma ),
302  \$ u, ldu, work( iu2 ), ldu2, vt, ldvt, work( ivt2 ),
303  \$ ldvt2, iwork( idxc ), iwork( coltyp ), work( iz ),
304  \$ info )
305  IF( info.NE.0 ) THEN
306  return
307  END IF
308 *
309 * Unscale.
310 *
311  CALL dlascl( '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 dlamrg( n1, n2, d, 1, -1, idxq )
318 *
319  return
320 *
321 * End of DLASD1
322 *
323  END