LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> \ingroup lasd7
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 dlasd7( 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 DOUBLE PRECISION ALPHA, BETA, C, S
289* ..
290* .. Array Arguments ..
291 INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
292 $ IDXQ( * ), PERM( * )
293 DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
294 $ vf( * ), vfw( * ), vl( * ), vlw( * ), z( * ),
295 $ zw( * )
296* ..
297*
298* =====================================================================
299*
300* .. Parameters ..
301 DOUBLE PRECISION ZERO, ONE, TWO, EIGHT
302 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
303 $ eight = 8.0d+0 )
304* ..
305* .. Local Scalars ..
306*
307 INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
308 $ NLP1, NLP2
309 DOUBLE PRECISION EPS, HLFTOL, TAU, TOL, Z1
310* ..
311* .. External Subroutines ..
312 EXTERNAL dcopy, dlamrg, drot, xerbla
313* ..
314* .. External Functions ..
315 DOUBLE PRECISION DLAMCH, DLAPY2
316 EXTERNAL DLAMCH, DLAPY2
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( 'DLASD7', -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 dlamrg( 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 = dlamch( '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 = dlapy2( 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 drot( 1, vf( jprev ), 1, vf( j ), 1, c, s )
491 CALL drot( 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 dcopy( 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 ) = dlapy2( 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 drot( 1, vf( m ), 1, vf( 1 ), 1, c, s )
558 CALL drot( 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 dcopy( k-1, zw( 2 ), 1, z( 2 ), 1 )
570 CALL dcopy( n-1, vfw( 2 ), 1, vf( 2 ), 1 )
571 CALL dcopy( n-1, vlw( 2 ), 1, vl( 2 ), 1 )
572*
573 RETURN
574*
575* End of DLASD7
576*
577 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
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:99
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:280
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92