LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlasd7.f
Go to the documentation of this file.
1 *> \brief \b DLASD7 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 DLASD7 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd7.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd7.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd7.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLASD7( 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 * DOUBLE PRECISION ALPHA, BETA, C, S
30 * ..
31 * .. Array Arguments ..
32 * INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
33 * $ IDXQ( * ), PERM( * )
34 * DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
35 * $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
36 * $ ZW( * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> DLASD7 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 *> DLASD7 is called from DLASD6.
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension ( M )
117 *> Workspace for Z.
118 *> \endverbatim
119 *>
120 *> \param[in,out] VF
121 *> \verbatim
122 *> VF is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension ( M )
133 *> Workspace for VF.
134 *> \endverbatim
135 *>
136 *> \param[in,out] VL
137 *> \verbatim
138 *> VL is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension ( M )
149 *> Workspace for VL.
150 *> \endverbatim
151 *>
152 *> \param[in] ALPHA
153 *> \verbatim
154 *> ALPHA is DOUBLE PRECISION
155 *> Contains the diagonal element associated with the added row.
156 *> \endverbatim
157 *>
158 *> \param[in] BETA
159 *> \verbatim
160 *> BETA is DOUBLE PRECISION
161 *> Contains the off-diagonal element associated with the added
162 *> row.
163 *> \endverbatim
164 *>
165 *> \param[out] DSIGMA
166 *> \verbatim
167 *> DSIGMA is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 *> \date September 2012
268 *
269 *> \ingroup auxOTHERauxiliary
270 *
271 *> \par Contributors:
272 * ==================
273 *>
274 *> Ming Gu and Huan Ren, Computer Science Division, University of
275 *> California at Berkeley, USA
276 *>
277 * =====================================================================
278  SUBROUTINE dlasd7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
279  $ vlw, alpha, beta, dsigma, idx, idxp, idxq,
280  $ perm, givptr, givcol, ldgcol, givnum, ldgnum,
281  $ c, s, info )
282 *
283 * -- LAPACK auxiliary routine (version 3.4.2) --
284 * -- LAPACK is a software package provided by Univ. of Tennessee, --
285 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
286 * September 2012
287 *
288 * .. Scalar Arguments ..
289  INTEGER givptr, icompq, info, k, ldgcol, ldgnum, nl,
290  $ nr, sqre
291  DOUBLE PRECISION alpha, beta, c, s
292 * ..
293 * .. Array Arguments ..
294  INTEGER givcol( ldgcol, * ), idx( * ), idxp( * ),
295  $ idxq( * ), perm( * )
296  DOUBLE PRECISION d( * ), dsigma( * ), givnum( ldgnum, * ),
297  $ vf( * ), vfw( * ), vl( * ), vlw( * ), z( * ),
298  $ zw( * )
299 * ..
300 *
301 * =====================================================================
302 *
303 * .. Parameters ..
304  DOUBLE PRECISION zero, one, two, eight
305  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
306  $ eight = 8.0d+0 )
307 * ..
308 * .. Local Scalars ..
309 *
310  INTEGER i, idxi, idxj, idxjp, j, jp, jprev, k2, m, n,
311  $ nlp1, nlp2
312  DOUBLE PRECISION eps, hlftol, tau, tol, z1
313 * ..
314 * .. External Subroutines ..
315  EXTERNAL dcopy, dlamrg, drot, xerbla
316 * ..
317 * .. External Functions ..
318  DOUBLE PRECISION dlamch, dlapy2
319  EXTERNAL dlamch, dlapy2
320 * ..
321 * .. Intrinsic Functions ..
322  INTRINSIC abs, max
323 * ..
324 * .. Executable Statements ..
325 *
326 * Test the input parameters.
327 *
328  info = 0
329  n = nl + nr + 1
330  m = n + sqre
331 *
332  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
333  info = -1
334  ELSE IF( nl.LT.1 ) THEN
335  info = -2
336  ELSE IF( nr.LT.1 ) THEN
337  info = -3
338  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
339  info = -4
340  ELSE IF( ldgcol.LT.n ) THEN
341  info = -22
342  ELSE IF( ldgnum.LT.n ) THEN
343  info = -24
344  END IF
345  IF( info.NE.0 ) THEN
346  CALL xerbla( 'DLASD7', -info )
347  return
348  END IF
349 *
350  nlp1 = nl + 1
351  nlp2 = nl + 2
352  IF( icompq.EQ.1 ) THEN
353  givptr = 0
354  END IF
355 *
356 * Generate the first part of the vector Z and move the singular
357 * values in the first part of D one position backward.
358 *
359  z1 = alpha*vl( nlp1 )
360  vl( nlp1 ) = zero
361  tau = vf( nlp1 )
362  DO 10 i = nl, 1, -1
363  z( i+1 ) = alpha*vl( i )
364  vl( i ) = zero
365  vf( i+1 ) = vf( i )
366  d( i+1 ) = d( i )
367  idxq( i+1 ) = idxq( i ) + 1
368  10 continue
369  vf( 1 ) = tau
370 *
371 * Generate the second part of the vector Z.
372 *
373  DO 20 i = nlp2, m
374  z( i ) = beta*vf( i )
375  vf( i ) = zero
376  20 continue
377 *
378 * Sort the singular values into increasing order
379 *
380  DO 30 i = nlp2, n
381  idxq( i ) = idxq( i ) + nlp1
382  30 continue
383 *
384 * DSIGMA, IDXC, IDXC, and ZW are used as storage space.
385 *
386  DO 40 i = 2, n
387  dsigma( i ) = d( idxq( i ) )
388  zw( i ) = z( idxq( i ) )
389  vfw( i ) = vf( idxq( i ) )
390  vlw( i ) = vl( idxq( i ) )
391  40 continue
392 *
393  CALL dlamrg( nl, nr, dsigma( 2 ), 1, 1, idx( 2 ) )
394 *
395  DO 50 i = 2, n
396  idxi = 1 + idx( i )
397  d( i ) = dsigma( idxi )
398  z( i ) = zw( idxi )
399  vf( i ) = vfw( idxi )
400  vl( i ) = vlw( idxi )
401  50 continue
402 *
403 * Calculate the allowable deflation tolerence
404 *
405  eps = dlamch( 'Epsilon' )
406  tol = max( abs( alpha ), abs( beta ) )
407  tol = eight*eight*eps*max( abs( d( n ) ), tol )
408 *
409 * There are 2 kinds of deflation -- first a value in the z-vector
410 * is small, second two (or more) singular values are very close
411 * together (their difference is small).
412 *
413 * If the value in the z-vector is small, we simply permute the
414 * array so that the corresponding singular value is moved to the
415 * end.
416 *
417 * If two values in the D-vector are close, we perform a two-sided
418 * rotation designed to make one of the corresponding z-vector
419 * entries zero, and then permute the array so that the deflated
420 * singular value is moved to the end.
421 *
422 * If there are multiple singular values then the problem deflates.
423 * Here the number of equal singular values are found. As each equal
424 * singular value is found, an elementary reflector is computed to
425 * rotate the corresponding singular subspace so that the
426 * corresponding components of Z are zero in this new basis.
427 *
428  k = 1
429  k2 = n + 1
430  DO 60 j = 2, n
431  IF( abs( z( j ) ).LE.tol ) THEN
432 *
433 * Deflate due to small z component.
434 *
435  k2 = k2 - 1
436  idxp( k2 ) = j
437  IF( j.EQ.n )
438  $ go to 100
439  ELSE
440  jprev = j
441  go to 70
442  END IF
443  60 continue
444  70 continue
445  j = jprev
446  80 continue
447  j = j + 1
448  IF( j.GT.n )
449  $ go to 90
450  IF( abs( z( j ) ).LE.tol ) THEN
451 *
452 * Deflate due to small z component.
453 *
454  k2 = k2 - 1
455  idxp( k2 ) = j
456  ELSE
457 *
458 * Check if singular values are close enough to allow deflation.
459 *
460  IF( abs( d( j )-d( jprev ) ).LE.tol ) THEN
461 *
462 * Deflation is possible.
463 *
464  s = z( jprev )
465  c = z( j )
466 *
467 * Find sqrt(a**2+b**2) without overflow or
468 * destructive underflow.
469 *
470  tau = dlapy2( c, s )
471  z( j ) = tau
472  z( jprev ) = zero
473  c = c / tau
474  s = -s / tau
475 *
476 * Record the appropriate Givens rotation
477 *
478  IF( icompq.EQ.1 ) THEN
479  givptr = givptr + 1
480  idxjp = idxq( idx( jprev )+1 )
481  idxj = idxq( idx( j )+1 )
482  IF( idxjp.LE.nlp1 ) THEN
483  idxjp = idxjp - 1
484  END IF
485  IF( idxj.LE.nlp1 ) THEN
486  idxj = idxj - 1
487  END IF
488  givcol( givptr, 2 ) = idxjp
489  givcol( givptr, 1 ) = idxj
490  givnum( givptr, 2 ) = c
491  givnum( givptr, 1 ) = s
492  END IF
493  CALL drot( 1, vf( jprev ), 1, vf( j ), 1, c, s )
494  CALL drot( 1, vl( jprev ), 1, vl( j ), 1, c, s )
495  k2 = k2 - 1
496  idxp( k2 ) = jprev
497  jprev = j
498  ELSE
499  k = k + 1
500  zw( k ) = z( jprev )
501  dsigma( k ) = d( jprev )
502  idxp( k ) = jprev
503  jprev = j
504  END IF
505  END IF
506  go to 80
507  90 continue
508 *
509 * Record the last singular value.
510 *
511  k = k + 1
512  zw( k ) = z( jprev )
513  dsigma( k ) = d( jprev )
514  idxp( k ) = jprev
515 *
516  100 continue
517 *
518 * Sort the singular values into DSIGMA. The singular values which
519 * were not deflated go into the first K slots of DSIGMA, except
520 * that DSIGMA(1) is treated separately.
521 *
522  DO 110 j = 2, n
523  jp = idxp( j )
524  dsigma( j ) = d( jp )
525  vfw( j ) = vf( jp )
526  vlw( j ) = vl( jp )
527  110 continue
528  IF( icompq.EQ.1 ) THEN
529  DO 120 j = 2, n
530  jp = idxp( j )
531  perm( j ) = idxq( idx( jp )+1 )
532  IF( perm( j ).LE.nlp1 ) THEN
533  perm( j ) = perm( j ) - 1
534  END IF
535  120 continue
536  END IF
537 *
538 * The deflated singular values go back into the last N - K slots of
539 * D.
540 *
541  CALL dcopy( n-k, dsigma( k+1 ), 1, d( k+1 ), 1 )
542 *
543 * Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
544 * VL(M).
545 *
546  dsigma( 1 ) = zero
547  hlftol = tol / two
548  IF( abs( dsigma( 2 ) ).LE.hlftol )
549  $ dsigma( 2 ) = hlftol
550  IF( m.GT.n ) THEN
551  z( 1 ) = dlapy2( z1, z( m ) )
552  IF( z( 1 ).LE.tol ) THEN
553  c = one
554  s = zero
555  z( 1 ) = tol
556  ELSE
557  c = z1 / z( 1 )
558  s = -z( m ) / z( 1 )
559  END IF
560  CALL drot( 1, vf( m ), 1, vf( 1 ), 1, c, s )
561  CALL drot( 1, vl( m ), 1, vl( 1 ), 1, c, s )
562  ELSE
563  IF( abs( z1 ).LE.tol ) THEN
564  z( 1 ) = tol
565  ELSE
566  z( 1 ) = z1
567  END IF
568  END IF
569 *
570 * Restore Z, VF, and VL.
571 *
572  CALL dcopy( k-1, zw( 2 ), 1, z( 2 ), 1 )
573  CALL dcopy( n-1, vfw( 2 ), 1, vf( 2 ), 1 )
574  CALL dcopy( n-1, vlw( 2 ), 1, vl( 2 ), 1 )
575 *
576  return
577 *
578 * End of DLASD7
579 *
580  END