LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slasd6.f
Go to the documentation of this file.
1 *> \brief \b SLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by appending a row. 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 SLASD6 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd6.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd6.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd6.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
22 * IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
23 * LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
24 * IWORK, INFO )
25 *
26 * .. Scalar Arguments ..
27 * INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
28 * $ NR, SQRE
29 * REAL ALPHA, BETA, C, S
30 * ..
31 * .. Array Arguments ..
32 * INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ),
33 * $ PERM( * )
34 * REAL D( * ), DIFL( * ), DIFR( * ),
35 * $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
36 * $ VF( * ), VL( * ), WORK( * ), Z( * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> SLASD6 computes the SVD of an updated upper bidiagonal matrix B
46 *> obtained by merging two smaller ones by appending a row. This
47 *> routine is used only for the problem which requires all singular
48 *> values and optionally singular vector matrices in factored form.
49 *> B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE.
50 *> A related subroutine, SLASD1, handles the case in which all singular
51 *> values and singular vectors of the bidiagonal matrix are desired.
52 *>
53 *> SLASD6 computes the SVD as follows:
54 *>
55 *> ( D1(in) 0 0 0 )
56 *> B = U(in) * ( Z1**T a Z2**T b ) * VT(in)
57 *> ( 0 0 D2(in) 0 )
58 *>
59 *> = U(out) * ( D(out) 0) * VT(out)
60 *>
61 *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
62 *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
63 *> elsewhere; and the entry b is empty if SQRE = 0.
64 *>
65 *> The singular values of B can be computed using D1, D2, the first
66 *> components of all the right singular vectors of the lower block, and
67 *> the last components of all the right singular vectors of the upper
68 *> block. These components are stored and updated in VF and VL,
69 *> respectively, in SLASD6. Hence U and VT are not explicitly
70 *> referenced.
71 *>
72 *> The singular values are stored in D. The algorithm consists of two
73 *> stages:
74 *>
75 *> The first stage consists of deflating the size of the problem
76 *> when there are multiple singular values or if there is a zero
77 *> in the Z vector. For each such occurence the dimension of the
78 *> secular equation problem is reduced by one. This stage is
79 *> performed by the routine SLASD7.
80 *>
81 *> The second stage consists of calculating the updated
82 *> singular values. This is done by finding the roots of the
83 *> secular equation via the routine SLASD4 (as called by SLASD8).
84 *> This routine also updates VF and VL and computes the distances
85 *> between the updated singular values and the old singular
86 *> values.
87 *>
88 *> SLASD6 is called from SLASDA.
89 *> \endverbatim
90 *
91 * Arguments:
92 * ==========
93 *
94 *> \param[in] ICOMPQ
95 *> \verbatim
96 *> ICOMPQ is INTEGER
97 *> Specifies whether singular vectors are to be computed in
98 *> factored form:
99 *> = 0: Compute singular values only.
100 *> = 1: Compute singular vectors in factored form as well.
101 *> \endverbatim
102 *>
103 *> \param[in] NL
104 *> \verbatim
105 *> NL is INTEGER
106 *> The row dimension of the upper block. NL >= 1.
107 *> \endverbatim
108 *>
109 *> \param[in] NR
110 *> \verbatim
111 *> NR is INTEGER
112 *> The row dimension of the lower block. NR >= 1.
113 *> \endverbatim
114 *>
115 *> \param[in] SQRE
116 *> \verbatim
117 *> SQRE is INTEGER
118 *> = 0: the lower block is an NR-by-NR square matrix.
119 *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
120 *>
121 *> The bidiagonal matrix has row dimension N = NL + NR + 1,
122 *> and column dimension M = N + SQRE.
123 *> \endverbatim
124 *>
125 *> \param[in,out] D
126 *> \verbatim
127 *> D is REAL array, dimension (NL+NR+1).
128 *> On entry D(1:NL,1:NL) contains the singular values of the
129 *> upper block, and D(NL+2:N) contains the singular values
130 *> of the lower block. On exit D(1:N) contains the singular
131 *> values of the modified matrix.
132 *> \endverbatim
133 *>
134 *> \param[in,out] VF
135 *> \verbatim
136 *> VF is REAL array, dimension (M)
137 *> On entry, VF(1:NL+1) contains the first components of all
138 *> right singular vectors of the upper block; and VF(NL+2:M)
139 *> contains the first components of all right singular vectors
140 *> of the lower block. On exit, VF contains the first components
141 *> of all right singular vectors of the bidiagonal matrix.
142 *> \endverbatim
143 *>
144 *> \param[in,out] VL
145 *> \verbatim
146 *> VL is REAL array, dimension (M)
147 *> On entry, VL(1:NL+1) contains the last components of all
148 *> right singular vectors of the upper block; and VL(NL+2:M)
149 *> contains the last components of all right singular vectors of
150 *> the lower block. On exit, VL contains the last components of
151 *> all right singular vectors of the bidiagonal matrix.
152 *> \endverbatim
153 *>
154 *> \param[in,out] ALPHA
155 *> \verbatim
156 *> ALPHA is REAL
157 *> Contains the diagonal element associated with the added row.
158 *> \endverbatim
159 *>
160 *> \param[in,out] BETA
161 *> \verbatim
162 *> BETA is REAL
163 *> Contains the off-diagonal element associated with the added
164 *> row.
165 *> \endverbatim
166 *>
167 *> \param[out] IDXQ
168 *> \verbatim
169 *> IDXQ is INTEGER array, dimension (N)
170 *> This contains the permutation which will reintegrate the
171 *> subproblem just solved back into sorted order, i.e.
172 *> D( IDXQ( I = 1, N ) ) will be in ascending order.
173 *> \endverbatim
174 *>
175 *> \param[out] PERM
176 *> \verbatim
177 *> PERM is INTEGER array, dimension ( N )
178 *> The permutations (from deflation and sorting) to be applied
179 *> to each block. Not referenced if ICOMPQ = 0.
180 *> \endverbatim
181 *>
182 *> \param[out] GIVPTR
183 *> \verbatim
184 *> GIVPTR is INTEGER
185 *> The number of Givens rotations which took place in this
186 *> subproblem. Not referenced if ICOMPQ = 0.
187 *> \endverbatim
188 *>
189 *> \param[out] GIVCOL
190 *> \verbatim
191 *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
192 *> Each pair of numbers indicates a pair of columns to take place
193 *> in a Givens rotation. Not referenced if ICOMPQ = 0.
194 *> \endverbatim
195 *>
196 *> \param[in] LDGCOL
197 *> \verbatim
198 *> LDGCOL is INTEGER
199 *> leading dimension of GIVCOL, must be at least N.
200 *> \endverbatim
201 *>
202 *> \param[out] GIVNUM
203 *> \verbatim
204 *> GIVNUM is REAL array, dimension ( LDGNUM, 2 )
205 *> Each number indicates the C or S value to be used in the
206 *> corresponding Givens rotation. Not referenced if ICOMPQ = 0.
207 *> \endverbatim
208 *>
209 *> \param[in] LDGNUM
210 *> \verbatim
211 *> LDGNUM is INTEGER
212 *> The leading dimension of GIVNUM and POLES, must be at least N.
213 *> \endverbatim
214 *>
215 *> \param[out] POLES
216 *> \verbatim
217 *> POLES is REAL array, dimension ( LDGNUM, 2 )
218 *> On exit, POLES(1,*) is an array containing the new singular
219 *> values obtained from solving the secular equation, and
220 *> POLES(2,*) is an array containing the poles in the secular
221 *> equation. Not referenced if ICOMPQ = 0.
222 *> \endverbatim
223 *>
224 *> \param[out] DIFL
225 *> \verbatim
226 *> DIFL is REAL array, dimension ( N )
227 *> On exit, DIFL(I) is the distance between I-th updated
228 *> (undeflated) singular value and the I-th (undeflated) old
229 *> singular value.
230 *> \endverbatim
231 *>
232 *> \param[out] DIFR
233 *> \verbatim
234 *> DIFR is REAL array,
235 *> dimension ( LDGNUM, 2 ) if ICOMPQ = 1 and
236 *> dimension ( N ) if ICOMPQ = 0.
237 *> On exit, DIFR(I, 1) is the distance between I-th updated
238 *> (undeflated) singular value and the I+1-th (undeflated) old
239 *> singular value.
240 *>
241 *> If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
242 *> normalizing factors for the right singular vector matrix.
243 *>
244 *> See SLASD8 for details on DIFL and DIFR.
245 *> \endverbatim
246 *>
247 *> \param[out] Z
248 *> \verbatim
249 *> Z is REAL array, dimension ( M )
250 *> The first elements of this array contain the components
251 *> of the deflation-adjusted updating row vector.
252 *> \endverbatim
253 *>
254 *> \param[out] K
255 *> \verbatim
256 *> K is INTEGER
257 *> Contains the dimension of the non-deflated matrix,
258 *> This is the order of the related secular equation. 1 <= K <=N.
259 *> \endverbatim
260 *>
261 *> \param[out] C
262 *> \verbatim
263 *> C is REAL
264 *> C contains garbage if SQRE =0 and the C-value of a Givens
265 *> rotation related to the right null space if SQRE = 1.
266 *> \endverbatim
267 *>
268 *> \param[out] S
269 *> \verbatim
270 *> S is REAL
271 *> S contains garbage if SQRE =0 and the S-value of a Givens
272 *> rotation related to the right null space if SQRE = 1.
273 *> \endverbatim
274 *>
275 *> \param[out] WORK
276 *> \verbatim
277 *> WORK is REAL array, dimension ( 4 * M )
278 *> \endverbatim
279 *>
280 *> \param[out] IWORK
281 *> \verbatim
282 *> IWORK is INTEGER array, dimension ( 3 * N )
283 *> \endverbatim
284 *>
285 *> \param[out] INFO
286 *> \verbatim
287 *> INFO is INTEGER
288 *> = 0: successful exit.
289 *> < 0: if INFO = -i, the i-th argument had an illegal value.
290 *> > 0: if INFO = 1, a singular value did not converge
291 *> \endverbatim
292 *
293 * Authors:
294 * ========
295 *
296 *> \author Univ. of Tennessee
297 *> \author Univ. of California Berkeley
298 *> \author Univ. of Colorado Denver
299 *> \author NAG Ltd.
300 *
301 *> \date September 2012
302 *
303 *> \ingroup auxOTHERauxiliary
304 *
305 *> \par Contributors:
306 * ==================
307 *>
308 *> Ming Gu and Huan Ren, Computer Science Division, University of
309 *> California at Berkeley, USA
310 *>
311 * =====================================================================
312  SUBROUTINE slasd6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
313  $ idxq, perm, givptr, givcol, ldgcol, givnum,
314  $ ldgnum, poles, difl, difr, z, k, c, s, work,
315  $ iwork, info )
316 *
317 * -- LAPACK auxiliary routine (version 3.4.2) --
318 * -- LAPACK is a software package provided by Univ. of Tennessee, --
319 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
320 * September 2012
321 *
322 * .. Scalar Arguments ..
323  INTEGER givptr, icompq, info, k, ldgcol, ldgnum, nl,
324  $ nr, sqre
325  REAL alpha, beta, c, s
326 * ..
327 * .. Array Arguments ..
328  INTEGER givcol( ldgcol, * ), idxq( * ), iwork( * ),
329  $ perm( * )
330  REAL d( * ), difl( * ), difr( * ),
331  $ givnum( ldgnum, * ), poles( ldgnum, * ),
332  $ vf( * ), vl( * ), work( * ), z( * )
333 * ..
334 *
335 * =====================================================================
336 *
337 * .. Parameters ..
338  REAL one, zero
339  parameter( one = 1.0e+0, zero = 0.0e+0 )
340 * ..
341 * .. Local Scalars ..
342  INTEGER i, idx, idxc, idxp, isigma, ivfw, ivlw, iw, m,
343  $ n, n1, n2
344  REAL orgnrm
345 * ..
346 * .. External Subroutines ..
347  EXTERNAL scopy, slamrg, slascl, slasd7, slasd8, xerbla
348 * ..
349 * .. Intrinsic Functions ..
350  INTRINSIC abs, max
351 * ..
352 * .. Executable Statements ..
353 *
354 * Test the input parameters.
355 *
356  info = 0
357  n = nl + nr + 1
358  m = n + sqre
359 *
360  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
361  info = -1
362  ELSE IF( nl.LT.1 ) THEN
363  info = -2
364  ELSE IF( nr.LT.1 ) THEN
365  info = -3
366  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
367  info = -4
368  ELSE IF( ldgcol.LT.n ) THEN
369  info = -14
370  ELSE IF( ldgnum.LT.n ) THEN
371  info = -16
372  END IF
373  IF( info.NE.0 ) THEN
374  CALL xerbla( 'SLASD6', -info )
375  return
376  END IF
377 *
378 * The following values are for bookkeeping purposes only. They are
379 * integer pointers which indicate the portion of the workspace
380 * used by a particular array in SLASD7 and SLASD8.
381 *
382  isigma = 1
383  iw = isigma + n
384  ivfw = iw + m
385  ivlw = ivfw + m
386 *
387  idx = 1
388  idxc = idx + n
389  idxp = idxc + n
390 *
391 * Scale.
392 *
393  orgnrm = max( abs( alpha ), abs( beta ) )
394  d( nl+1 ) = zero
395  DO 10 i = 1, n
396  IF( abs( d( i ) ).GT.orgnrm ) THEN
397  orgnrm = abs( d( i ) )
398  END IF
399  10 continue
400  CALL slascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
401  alpha = alpha / orgnrm
402  beta = beta / orgnrm
403 *
404 * Sort and Deflate singular values.
405 *
406  CALL slasd7( icompq, nl, nr, sqre, k, d, z, work( iw ), vf,
407  $ work( ivfw ), vl, work( ivlw ), alpha, beta,
408  $ work( isigma ), iwork( idx ), iwork( idxp ), idxq,
409  $ perm, givptr, givcol, ldgcol, givnum, ldgnum, c, s,
410  $ info )
411 *
412 * Solve Secular Equation, compute DIFL, DIFR, and update VF, VL.
413 *
414  CALL slasd8( icompq, k, d, z, vf, vl, difl, difr, ldgnum,
415  $ work( isigma ), work( iw ), info )
416 *
417 * Handle error returned
418 *
419  IF( info.NE.0 ) THEN
420  CALL xerbla( 'SLASD8', -info )
421  return
422  END IF
423 *
424 * Save the poles if ICOMPQ = 1.
425 *
426  IF( icompq.EQ.1 ) THEN
427  CALL scopy( k, d, 1, poles( 1, 1 ), 1 )
428  CALL scopy( k, work( isigma ), 1, poles( 1, 2 ), 1 )
429  END IF
430 *
431 * Unscale.
432 *
433  CALL slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
434 *
435 * Prepare the IDXQ sorting permutation.
436 *
437  n1 = k
438  n2 = n - k
439  CALL slamrg( n1, n2, d, 1, -1, idxq )
440 *
441  return
442 *
443 * End of SLASD6
444 *
445  END