LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
slasd7.f
Go to the documentation of this file.
1 *> \brief \b SLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem. 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 SLASD7 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd7.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd7.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd7.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
22 * VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
23 * PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
24 * C, S, 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, * ), IDX( * ), IDXP( * ),
33 * $ IDXQ( * ), PERM( * )
34 * REAL D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
35 * $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
36 * $ ZW( * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> SLASD7 merges the two sets of singular values together into a single
46 *> sorted set. Then it tries to deflate the size of the problem. There
47 *> are two ways in which deflation can occur: when two or more singular
48 *> values are close together or if there is a tiny entry in the Z
49 *> vector. For each such occurrence the order of the related
50 *> secular equation problem is reduced by one.
51 *>
52 *> SLASD7 is called from SLASD6.
53 *> \endverbatim
54 *
55 * Arguments:
56 * ==========
57 *
58 *> \param[in] ICOMPQ
59 *> \verbatim
60 *> ICOMPQ is INTEGER
61 *> Specifies whether singular vectors are to be computed
62 *> in compact form, as follows:
63 *> = 0: Compute singular values only.
64 *> = 1: Compute singular vectors of upper
65 *> bidiagonal matrix in compact form.
66 *> \endverbatim
67 *>
68 *> \param[in] NL
69 *> \verbatim
70 *> NL is INTEGER
71 *> The row dimension of the upper block. NL >= 1.
72 *> \endverbatim
73 *>
74 *> \param[in] NR
75 *> \verbatim
76 *> NR is INTEGER
77 *> The row dimension of the lower block. NR >= 1.
78 *> \endverbatim
79 *>
80 *> \param[in] SQRE
81 *> \verbatim
82 *> SQRE is INTEGER
83 *> = 0: the lower block is an NR-by-NR square matrix.
84 *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
85 *>
86 *> The bidiagonal matrix has
87 *> N = NL + NR + 1 rows and
88 *> M = N + SQRE >= N columns.
89 *> \endverbatim
90 *>
91 *> \param[out] K
92 *> \verbatim
93 *> K is INTEGER
94 *> Contains the dimension of the non-deflated matrix, this is
95 *> the order of the related secular equation. 1 <= K <=N.
96 *> \endverbatim
97 *>
98 *> \param[in,out] D
99 *> \verbatim
100 *> D is REAL array, dimension ( N )
101 *> On entry D contains the singular values of the two submatrices
102 *> to be combined. On exit D contains the trailing (N-K) updated
103 *> singular values (those which were deflated) sorted into
104 *> increasing order.
105 *> \endverbatim
106 *>
107 *> \param[out] Z
108 *> \verbatim
109 *> Z is REAL array, dimension ( M )
110 *> On exit Z contains the updating row vector in the secular
111 *> equation.
112 *> \endverbatim
113 *>
114 *> \param[out] ZW
115 *> \verbatim
116 *> ZW is REAL array, dimension ( M )
117 *> Workspace for Z.
118 *> \endverbatim
119 *>
120 *> \param[in,out] VF
121 *> \verbatim
122 *> VF is REAL array, dimension ( M )
123 *> On entry, VF(1:NL+1) contains the first components of all
124 *> right singular vectors of the upper block; and VF(NL+2:M)
125 *> contains the first components of all right singular vectors
126 *> of the lower block. On exit, VF contains the first components
127 *> of all right singular vectors of the bidiagonal matrix.
128 *> \endverbatim
129 *>
130 *> \param[out] VFW
131 *> \verbatim
132 *> VFW is REAL array, dimension ( M )
133 *> Workspace for VF.
134 *> \endverbatim
135 *>
136 *> \param[in,out] VL
137 *> \verbatim
138 *> VL is REAL array, dimension ( M )
139 *> On entry, VL(1:NL+1) contains the last components of all
140 *> right singular vectors of the upper block; and VL(NL+2:M)
141 *> contains the last components of all right singular vectors
142 *> of the lower block. On exit, VL contains the last components
143 *> of all right singular vectors of the bidiagonal matrix.
144 *> \endverbatim
145 *>
146 *> \param[out] VLW
147 *> \verbatim
148 *> VLW is REAL array, dimension ( M )
149 *> Workspace for VL.
150 *> \endverbatim
151 *>
152 *> \param[in] ALPHA
153 *> \verbatim
154 *> ALPHA is REAL
155 *> Contains the diagonal element associated with the added row.
156 *> \endverbatim
157 *>
158 *> \param[in] BETA
159 *> \verbatim
160 *> BETA is REAL
161 *> Contains the off-diagonal element associated with the added
162 *> row.
163 *> \endverbatim
164 *>
165 *> \param[out] DSIGMA
166 *> \verbatim
167 *> DSIGMA is REAL array, dimension ( N )
168 *> Contains a copy of the diagonal elements (K-1 singular values
169 *> and one zero) in the secular equation.
170 *> \endverbatim
171 *>
172 *> \param[out] IDX
173 *> \verbatim
174 *> IDX is INTEGER array, dimension ( N )
175 *> This will contain the permutation used to sort the contents of
176 *> D into ascending order.
177 *> \endverbatim
178 *>
179 *> \param[out] IDXP
180 *> \verbatim
181 *> IDXP is INTEGER array, dimension ( N )
182 *> This will contain the permutation used to place deflated
183 *> values of D at the end of the array. On output IDXP(2:K)
184 *> points to the nondeflated D-values and IDXP(K+1:N)
185 *> points to the deflated singular values.
186 *> \endverbatim
187 *>
188 *> \param[in] IDXQ
189 *> \verbatim
190 *> IDXQ is INTEGER array, dimension ( N )
191 *> This contains the permutation which separately sorts the two
192 *> sub-problems in D into ascending order. Note that entries in
193 *> the first half of this permutation must first be moved one
194 *> position backward; and entries in the second half
195 *> must first have NL+1 added to their values.
196 *> \endverbatim
197 *>
198 *> \param[out] PERM
199 *> \verbatim
200 *> PERM is INTEGER array, dimension ( N )
201 *> The permutations (from deflation and sorting) to be applied
202 *> to each singular block. Not referenced if ICOMPQ = 0.
203 *> \endverbatim
204 *>
205 *> \param[out] GIVPTR
206 *> \verbatim
207 *> GIVPTR is INTEGER
208 *> The number of Givens rotations which took place in this
209 *> subproblem. Not referenced if ICOMPQ = 0.
210 *> \endverbatim
211 *>
212 *> \param[out] GIVCOL
213 *> \verbatim
214 *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
215 *> Each pair of numbers indicates a pair of columns to take place
216 *> in a Givens rotation. Not referenced if ICOMPQ = 0.
217 *> \endverbatim
218 *>
219 *> \param[in] LDGCOL
220 *> \verbatim
221 *> LDGCOL is INTEGER
222 *> The leading dimension of GIVCOL, must be at least N.
223 *> \endverbatim
224 *>
225 *> \param[out] GIVNUM
226 *> \verbatim
227 *> GIVNUM is REAL array, dimension ( LDGNUM, 2 )
228 *> Each number indicates the C or S value to be used in the
229 *> corresponding Givens rotation. Not referenced if ICOMPQ = 0.
230 *> \endverbatim
231 *>
232 *> \param[in] LDGNUM
233 *> \verbatim
234 *> LDGNUM is INTEGER
235 *> The leading dimension of GIVNUM, must be at least N.
236 *> \endverbatim
237 *>
238 *> \param[out] C
239 *> \verbatim
240 *> C is REAL
241 *> C contains garbage if SQRE =0 and the C-value of a Givens
242 *> rotation related to the right null space if SQRE = 1.
243 *> \endverbatim
244 *>
245 *> \param[out] S
246 *> \verbatim
247 *> S is REAL
248 *> S contains garbage if SQRE =0 and the S-value of a Givens
249 *> rotation related to the right null space if SQRE = 1.
250 *> \endverbatim
251 *>
252 *> \param[out] INFO
253 *> \verbatim
254 *> INFO is INTEGER
255 *> = 0: successful exit.
256 *> < 0: if INFO = -i, the i-th argument had an illegal value.
257 *> \endverbatim
258 *
259 * Authors:
260 * ========
261 *
262 *> \author Univ. of Tennessee
263 *> \author Univ. of California Berkeley
264 *> \author Univ. of Colorado Denver
265 *> \author NAG Ltd.
266 *
267 *> \ingroup OTHERauxiliary
268 *
269 *> \par Contributors:
270 * ==================
271 *>
272 *> Ming Gu and Huan Ren, Computer Science Division, University of
273 *> California at Berkeley, USA
274 *>
275 * =====================================================================
276  SUBROUTINE slasd7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
277  $ VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
278  $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
279  $ C, S, INFO )
280 *
281 * -- LAPACK auxiliary routine --
282 * -- LAPACK is a software package provided by Univ. of Tennessee, --
283 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
284 *
285 * .. Scalar Arguments ..
286  INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
287  $ NR, SQRE
288  REAL ALPHA, BETA, C, S
289 * ..
290 * .. Array Arguments ..
291  INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
292  $ IDXQ( * ), PERM( * )
293  REAL D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
294  $ vf( * ), vfw( * ), vl( * ), vlw( * ), z( * ),
295  $ zw( * )
296 * ..
297 *
298 * =====================================================================
299 *
300 * .. Parameters ..
301  REAL ZERO, ONE, TWO, EIGHT
302  PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0, two = 2.0e+0,
303  $ eight = 8.0e+0 )
304 * ..
305 * .. Local Scalars ..
306 *
307  INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
308  $ NLP1, NLP2
309  REAL EPS, HLFTOL, TAU, TOL, Z1
310 * ..
311 * .. External Subroutines ..
312  EXTERNAL scopy, slamrg, srot, xerbla
313 * ..
314 * .. External Functions ..
315  REAL SLAMCH, SLAPY2
316  EXTERNAL SLAMCH, SLAPY2
317 * ..
318 * .. Intrinsic Functions ..
319  INTRINSIC abs, max
320 * ..
321 * .. Executable Statements ..
322 *
323 * Test the input parameters.
324 *
325  info = 0
326  n = nl + nr + 1
327  m = n + sqre
328 *
329  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
330  info = -1
331  ELSE IF( nl.LT.1 ) THEN
332  info = -2
333  ELSE IF( nr.LT.1 ) THEN
334  info = -3
335  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
336  info = -4
337  ELSE IF( ldgcol.LT.n ) THEN
338  info = -22
339  ELSE IF( ldgnum.LT.n ) THEN
340  info = -24
341  END IF
342  IF( info.NE.0 ) THEN
343  CALL xerbla( 'SLASD7', -info )
344  RETURN
345  END IF
346 *
347  nlp1 = nl + 1
348  nlp2 = nl + 2
349  IF( icompq.EQ.1 ) THEN
350  givptr = 0
351  END IF
352 *
353 * Generate the first part of the vector Z and move the singular
354 * values in the first part of D one position backward.
355 *
356  z1 = alpha*vl( nlp1 )
357  vl( nlp1 ) = zero
358  tau = vf( nlp1 )
359  DO 10 i = nl, 1, -1
360  z( i+1 ) = alpha*vl( i )
361  vl( i ) = zero
362  vf( i+1 ) = vf( i )
363  d( i+1 ) = d( i )
364  idxq( i+1 ) = idxq( i ) + 1
365  10 CONTINUE
366  vf( 1 ) = tau
367 *
368 * Generate the second part of the vector Z.
369 *
370  DO 20 i = nlp2, m
371  z( i ) = beta*vf( i )
372  vf( i ) = zero
373  20 CONTINUE
374 *
375 * Sort the singular values into increasing order
376 *
377  DO 30 i = nlp2, n
378  idxq( i ) = idxq( i ) + nlp1
379  30 CONTINUE
380 *
381 * DSIGMA, IDXC, IDXC, and ZW are used as storage space.
382 *
383  DO 40 i = 2, n
384  dsigma( i ) = d( idxq( i ) )
385  zw( i ) = z( idxq( i ) )
386  vfw( i ) = vf( idxq( i ) )
387  vlw( i ) = vl( idxq( i ) )
388  40 CONTINUE
389 *
390  CALL slamrg( nl, nr, dsigma( 2 ), 1, 1, idx( 2 ) )
391 *
392  DO 50 i = 2, n
393  idxi = 1 + idx( i )
394  d( i ) = dsigma( idxi )
395  z( i ) = zw( idxi )
396  vf( i ) = vfw( idxi )
397  vl( i ) = vlw( idxi )
398  50 CONTINUE
399 *
400 * Calculate the allowable deflation tolerance
401 *
402  eps = slamch( 'Epsilon' )
403  tol = max( abs( alpha ), abs( beta ) )
404  tol = eight*eight*eps*max( abs( d( n ) ), tol )
405 *
406 * There are 2 kinds of deflation -- first a value in the z-vector
407 * is small, second two (or more) singular values are very close
408 * together (their difference is small).
409 *
410 * If the value in the z-vector is small, we simply permute the
411 * array so that the corresponding singular value is moved to the
412 * end.
413 *
414 * If two values in the D-vector are close, we perform a two-sided
415 * rotation designed to make one of the corresponding z-vector
416 * entries zero, and then permute the array so that the deflated
417 * singular value is moved to the end.
418 *
419 * If there are multiple singular values then the problem deflates.
420 * Here the number of equal singular values are found. As each equal
421 * singular value is found, an elementary reflector is computed to
422 * rotate the corresponding singular subspace so that the
423 * corresponding components of Z are zero in this new basis.
424 *
425  k = 1
426  k2 = n + 1
427  DO 60 j = 2, n
428  IF( abs( z( j ) ).LE.tol ) THEN
429 *
430 * Deflate due to small z component.
431 *
432  k2 = k2 - 1
433  idxp( k2 ) = j
434  IF( j.EQ.n )
435  $ GO TO 100
436  ELSE
437  jprev = j
438  GO TO 70
439  END IF
440  60 CONTINUE
441  70 CONTINUE
442  j = jprev
443  80 CONTINUE
444  j = j + 1
445  IF( j.GT.n )
446  $ GO TO 90
447  IF( abs( z( j ) ).LE.tol ) THEN
448 *
449 * Deflate due to small z component.
450 *
451  k2 = k2 - 1
452  idxp( k2 ) = j
453  ELSE
454 *
455 * Check if singular values are close enough to allow deflation.
456 *
457  IF( abs( d( j )-d( jprev ) ).LE.tol ) THEN
458 *
459 * Deflation is possible.
460 *
461  s = z( jprev )
462  c = z( j )
463 *
464 * Find sqrt(a**2+b**2) without overflow or
465 * destructive underflow.
466 *
467  tau = slapy2( c, s )
468  z( j ) = tau
469  z( jprev ) = zero
470  c = c / tau
471  s = -s / tau
472 *
473 * Record the appropriate Givens rotation
474 *
475  IF( icompq.EQ.1 ) THEN
476  givptr = givptr + 1
477  idxjp = idxq( idx( jprev )+1 )
478  idxj = idxq( idx( j )+1 )
479  IF( idxjp.LE.nlp1 ) THEN
480  idxjp = idxjp - 1
481  END IF
482  IF( idxj.LE.nlp1 ) THEN
483  idxj = idxj - 1
484  END IF
485  givcol( givptr, 2 ) = idxjp
486  givcol( givptr, 1 ) = idxj
487  givnum( givptr, 2 ) = c
488  givnum( givptr, 1 ) = s
489  END IF
490  CALL srot( 1, vf( jprev ), 1, vf( j ), 1, c, s )
491  CALL srot( 1, vl( jprev ), 1, vl( j ), 1, c, s )
492  k2 = k2 - 1
493  idxp( k2 ) = jprev
494  jprev = j
495  ELSE
496  k = k + 1
497  zw( k ) = z( jprev )
498  dsigma( k ) = d( jprev )
499  idxp( k ) = jprev
500  jprev = j
501  END IF
502  END IF
503  GO TO 80
504  90 CONTINUE
505 *
506 * Record the last singular value.
507 *
508  k = k + 1
509  zw( k ) = z( jprev )
510  dsigma( k ) = d( jprev )
511  idxp( k ) = jprev
512 *
513  100 CONTINUE
514 *
515 * Sort the singular values into DSIGMA. The singular values which
516 * were not deflated go into the first K slots of DSIGMA, except
517 * that DSIGMA(1) is treated separately.
518 *
519  DO 110 j = 2, n
520  jp = idxp( j )
521  dsigma( j ) = d( jp )
522  vfw( j ) = vf( jp )
523  vlw( j ) = vl( jp )
524  110 CONTINUE
525  IF( icompq.EQ.1 ) THEN
526  DO 120 j = 2, n
527  jp = idxp( j )
528  perm( j ) = idxq( idx( jp )+1 )
529  IF( perm( j ).LE.nlp1 ) THEN
530  perm( j ) = perm( j ) - 1
531  END IF
532  120 CONTINUE
533  END IF
534 *
535 * The deflated singular values go back into the last N - K slots of
536 * D.
537 *
538  CALL scopy( n-k, dsigma( k+1 ), 1, d( k+1 ), 1 )
539 *
540 * Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
541 * VL(M).
542 *
543  dsigma( 1 ) = zero
544  hlftol = tol / two
545  IF( abs( dsigma( 2 ) ).LE.hlftol )
546  $ dsigma( 2 ) = hlftol
547  IF( m.GT.n ) THEN
548  z( 1 ) = slapy2( z1, z( m ) )
549  IF( z( 1 ).LE.tol ) THEN
550  c = one
551  s = zero
552  z( 1 ) = tol
553  ELSE
554  c = z1 / z( 1 )
555  s = -z( m ) / z( 1 )
556  END IF
557  CALL srot( 1, vf( m ), 1, vf( 1 ), 1, c, s )
558  CALL srot( 1, vl( m ), 1, vl( 1 ), 1, c, s )
559  ELSE
560  IF( abs( z1 ).LE.tol ) THEN
561  z( 1 ) = tol
562  ELSE
563  z( 1 ) = z1
564  END IF
565  END IF
566 *
567 * Restore Z, VF, and VL.
568 *
569  CALL scopy( k-1, zw( 2 ), 1, z( 2 ), 1 )
570  CALL scopy( n-1, vfw( 2 ), 1, vf( 2 ), 1 )
571  CALL scopy( n-1, vlw( 2 ), 1, vl( 2 ), 1 )
572 *
573  RETURN
574 *
575 * End of SLASD7
576 *
577  END
subroutine slasd7(ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL, VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, INFO)
SLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to def...
Definition: slasd7.f:280
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slamrg(N1, N2, A, STRD1, STRD2, INDEX)
SLAMRG creates a permutation list to merge the entries of two independently sorted sets into a single...
Definition: slamrg.f:99
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82