LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlamrg(N1, N2, A, DTRD1, DTRD2, INDEX)
DLAMRG creates a permutation list to merge the entries of two independently sorted sets into a single...
Definition: dlamrg.f:101
subroutine dlasd7(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)
DLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to def...
Definition: dlasd7.f:282
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53