LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlasd7()

subroutine dlasd7 ( integer icompq,
integer nl,
integer nr,
integer sqre,
integer k,
double precision, dimension( * ) d,
double precision, dimension( * ) z,
double precision, dimension( * ) zw,
double precision, dimension( * ) vf,
double precision, dimension( * ) vfw,
double precision, dimension( * ) vl,
double precision, dimension( * ) vlw,
double precision alpha,
double precision beta,
double precision, dimension( * ) dsigma,
integer, dimension( * ) idx,
integer, dimension( * ) idxp,
integer, dimension( * ) idxq,
integer, dimension( * ) perm,
integer givptr,
integer, dimension( ldgcol, * ) givcol,
integer ldgcol,
double precision, dimension( ldgnum, * ) givnum,
integer ldgnum,
double precision c,
double precision s,
integer info )

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.

Download DLASD7 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DLASD7 merges the two sets of singular values together into a single
!> sorted set. Then it tries to deflate the size of the problem. There
!> are two ways in which deflation can occur:  when two or more singular
!> values are close together or if there is a tiny entry in the Z
!> vector. For each such occurrence the order of the related
!> secular equation problem is reduced by one.
!>
!> DLASD7 is called from DLASD6.
!>
Parameters
[in]ICOMPQ
!>          ICOMPQ is INTEGER
!>          Specifies whether singular vectors are to be computed
!>          in compact form, as follows:
!>          = 0: Compute singular values only.
!>          = 1: Compute singular vectors of upper
!>               bidiagonal matrix in compact form.
!>
[in]NL
!>          NL is INTEGER
!>         The row dimension of the upper block. NL >= 1.
!>
[in]NR
!>          NR is INTEGER
!>         The row dimension of the lower block. NR >= 1.
!>
[in]SQRE
!>          SQRE is INTEGER
!>         = 0: the lower block is an NR-by-NR square matrix.
!>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
!>
!>         The bidiagonal matrix has
!>         N = NL + NR + 1 rows and
!>         M = N + SQRE >= N columns.
!>
[out]K
!>          K is INTEGER
!>         Contains the dimension of the non-deflated matrix, this is
!>         the order of the related secular equation. 1 <= K <=N.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension ( N )
!>         On entry D contains the singular values of the two submatrices
!>         to be combined. On exit D contains the trailing (N-K) updated
!>         singular values (those which were deflated) sorted into
!>         increasing order.
!>
[out]Z
!>          Z is DOUBLE PRECISION array, dimension ( M )
!>         On exit Z contains the updating row vector in the secular
!>         equation.
!>
[out]ZW
!>          ZW is DOUBLE PRECISION array, dimension ( M )
!>         Workspace for Z.
!>
[in,out]VF
!>          VF is DOUBLE PRECISION array, dimension ( M )
!>         On entry, VF(1:NL+1) contains the first components of all
!>         right singular vectors of the upper block; and VF(NL+2:M)
!>         contains the first components of all right singular vectors
!>         of the lower block. On exit, VF contains the first components
!>         of all right singular vectors of the bidiagonal matrix.
!>
[out]VFW
!>          VFW is DOUBLE PRECISION array, dimension ( M )
!>         Workspace for VF.
!>
[in,out]VL
!>          VL is DOUBLE PRECISION array, dimension ( M )
!>         On entry, VL(1:NL+1) contains the  last components of all
!>         right singular vectors of the upper block; and VL(NL+2:M)
!>         contains the last components of all right singular vectors
!>         of the lower block. On exit, VL contains the last components
!>         of all right singular vectors of the bidiagonal matrix.
!>
[out]VLW
!>          VLW is DOUBLE PRECISION array, dimension ( M )
!>         Workspace for VL.
!>
[in]ALPHA
!>          ALPHA is DOUBLE PRECISION
!>         Contains the diagonal element associated with the added row.
!>
[in]BETA
!>          BETA is DOUBLE PRECISION
!>         Contains the off-diagonal element associated with the added
!>         row.
!>
[out]DSIGMA
!>          DSIGMA is DOUBLE PRECISION array, dimension ( N )
!>         Contains a copy of the diagonal elements (K-1 singular values
!>         and one zero) in the secular equation.
!>
[out]IDX
!>          IDX is INTEGER array, dimension ( N )
!>         This will contain the permutation used to sort the contents of
!>         D into ascending order.
!>
[out]IDXP
!>          IDXP is INTEGER array, dimension ( N )
!>         This will contain the permutation used to place deflated
!>         values of D at the end of the array. On output IDXP(2:K)
!>         points to the nondeflated D-values and IDXP(K+1:N)
!>         points to the deflated singular values.
!>
[in]IDXQ
!>          IDXQ is INTEGER array, dimension ( N )
!>         This contains the permutation which separately sorts the two
!>         sub-problems in D into ascending order.  Note that entries in
!>         the first half of this permutation must first be moved one
!>         position backward; and entries in the second half
!>         must first have NL+1 added to their values.
!>
[out]PERM
!>          PERM is INTEGER array, dimension ( N )
!>         The permutations (from deflation and sorting) to be applied
!>         to each singular block. Not referenced if ICOMPQ = 0.
!>
[out]GIVPTR
!>          GIVPTR is INTEGER
!>         The number of Givens rotations which took place in this
!>         subproblem. Not referenced if ICOMPQ = 0.
!>
[out]GIVCOL
!>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
!>         Each pair of numbers indicates a pair of columns to take place
!>         in a Givens rotation. Not referenced if ICOMPQ = 0.
!>
[in]LDGCOL
!>          LDGCOL is INTEGER
!>         The leading dimension of GIVCOL, must be at least N.
!>
[out]GIVNUM
!>          GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
!>         Each number indicates the C or S value to be used in the
!>         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
!>
[in]LDGNUM
!>          LDGNUM is INTEGER
!>         The leading dimension of GIVNUM, must be at least N.
!>
[out]C
!>          C is DOUBLE PRECISION
!>         C contains garbage if SQRE =0 and the C-value of a Givens
!>         rotation related to the right null space if SQRE = 1.
!>
[out]S
!>          S is DOUBLE PRECISION
!>         S contains garbage if SQRE =0 and the S-value of a Givens
!>         rotation related to the right null space if SQRE = 1.
!>
[out]INFO
!>          INFO is INTEGER
!>         = 0:  successful exit.
!>         < 0:  if INFO = -i, the i-th argument had an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 274 of file dlasd7.f.

279*
280* -- LAPACK auxiliary routine --
281* -- LAPACK is a software package provided by Univ. of Tennessee, --
282* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
283*
284* .. Scalar Arguments ..
285 INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
286 $ NR, SQRE
287 DOUBLE PRECISION ALPHA, BETA, C, S
288* ..
289* .. Array Arguments ..
290 INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
291 $ IDXQ( * ), PERM( * )
292 DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
293 $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
294 $ ZW( * )
295* ..
296*
297* =====================================================================
298*
299* .. Parameters ..
300 DOUBLE PRECISION ZERO, ONE, TWO, EIGHT
301 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
302 $ eight = 8.0d+0 )
303* ..
304* .. Local Scalars ..
305*
306 INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
307 $ NLP1, NLP2
308 DOUBLE PRECISION EPS, HLFTOL, TAU, TOL, Z1
309* ..
310* .. External Subroutines ..
311 EXTERNAL dcopy, dlamrg, drot, xerbla
312* ..
313* .. External Functions ..
314 DOUBLE PRECISION DLAMCH, DLAPY2
315 EXTERNAL dlamch, dlapy2
316* ..
317* .. Intrinsic Functions ..
318 INTRINSIC abs, max
319* ..
320* .. Executable Statements ..
321*
322* Test the input parameters.
323*
324 info = 0
325 n = nl + nr + 1
326 m = n + sqre
327*
328 IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
329 info = -1
330 ELSE IF( nl.LT.1 ) THEN
331 info = -2
332 ELSE IF( nr.LT.1 ) THEN
333 info = -3
334 ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
335 info = -4
336 ELSE IF( ldgcol.LT.n ) THEN
337 info = -22
338 ELSE IF( ldgnum.LT.n ) THEN
339 info = -24
340 END IF
341 IF( info.NE.0 ) THEN
342 CALL xerbla( 'DLASD7', -info )
343 RETURN
344 END IF
345*
346 nlp1 = nl + 1
347 nlp2 = nl + 2
348 IF( icompq.EQ.1 ) THEN
349 givptr = 0
350 END IF
351*
352* Generate the first part of the vector Z and move the singular
353* values in the first part of D one position backward.
354*
355 z1 = alpha*vl( nlp1 )
356 vl( nlp1 ) = zero
357 tau = vf( nlp1 )
358 DO 10 i = nl, 1, -1
359 z( i+1 ) = alpha*vl( i )
360 vl( i ) = zero
361 vf( i+1 ) = vf( i )
362 d( i+1 ) = d( i )
363 idxq( i+1 ) = idxq( i ) + 1
364 10 CONTINUE
365 vf( 1 ) = tau
366*
367* Generate the second part of the vector Z.
368*
369 DO 20 i = nlp2, m
370 z( i ) = beta*vf( i )
371 vf( i ) = zero
372 20 CONTINUE
373*
374* Sort the singular values into increasing order
375*
376 DO 30 i = nlp2, n
377 idxq( i ) = idxq( i ) + nlp1
378 30 CONTINUE
379*
380* DSIGMA, IDXC, IDXC, and ZW are used as storage space.
381*
382 DO 40 i = 2, n
383 dsigma( i ) = d( idxq( i ) )
384 zw( i ) = z( idxq( i ) )
385 vfw( i ) = vf( idxq( i ) )
386 vlw( i ) = vl( idxq( i ) )
387 40 CONTINUE
388*
389 CALL dlamrg( nl, nr, dsigma( 2 ), 1, 1, idx( 2 ) )
390*
391 DO 50 i = 2, n
392 idxi = 1 + idx( i )
393 d( i ) = dsigma( idxi )
394 z( i ) = zw( idxi )
395 vf( i ) = vfw( idxi )
396 vl( i ) = vlw( idxi )
397 50 CONTINUE
398*
399* Calculate the allowable deflation tolerance
400*
401 eps = dlamch( 'Epsilon' )
402 tol = max( abs( alpha ), abs( beta ) )
403 tol = eight*eight*eps*max( abs( d( n ) ), tol )
404*
405* There are 2 kinds of deflation -- first a value in the z-vector
406* is small, second two (or more) singular values are very close
407* together (their difference is small).
408*
409* If the value in the z-vector is small, we simply permute the
410* array so that the corresponding singular value is moved to the
411* end.
412*
413* If two values in the D-vector are close, we perform a two-sided
414* rotation designed to make one of the corresponding z-vector
415* entries zero, and then permute the array so that the deflated
416* singular value is moved to the end.
417*
418* If there are multiple singular values then the problem deflates.
419* Here the number of equal singular values are found. As each equal
420* singular value is found, an elementary reflector is computed to
421* rotate the corresponding singular subspace so that the
422* corresponding components of Z are zero in this new basis.
423*
424 k = 1
425 k2 = n + 1
426 DO 60 j = 2, n
427 IF( abs( z( j ) ).LE.tol ) THEN
428*
429* Deflate due to small z component.
430*
431 k2 = k2 - 1
432 idxp( k2 ) = j
433 IF( j.EQ.n )
434 $ GO TO 100
435 ELSE
436 jprev = j
437 GO TO 70
438 END IF
439 60 CONTINUE
440 70 CONTINUE
441 j = jprev
442 80 CONTINUE
443 j = j + 1
444 IF( j.GT.n )
445 $ GO TO 90
446 IF( abs( z( j ) ).LE.tol ) THEN
447*
448* Deflate due to small z component.
449*
450 k2 = k2 - 1
451 idxp( k2 ) = j
452 ELSE
453*
454* Check if singular values are close enough to allow deflation.
455*
456 IF( abs( d( j )-d( jprev ) ).LE.tol ) THEN
457*
458* Deflation is possible.
459*
460 s = z( jprev )
461 c = z( j )
462*
463* Find sqrt(a**2+b**2) without overflow or
464* destructive underflow.
465*
466 tau = dlapy2( c, s )
467 z( j ) = tau
468 z( jprev ) = zero
469 c = c / tau
470 s = -s / tau
471*
472* Record the appropriate Givens rotation
473*
474 IF( icompq.EQ.1 ) THEN
475 givptr = givptr + 1
476 idxjp = idxq( idx( jprev )+1 )
477 idxj = idxq( idx( j )+1 )
478 IF( idxjp.LE.nlp1 ) THEN
479 idxjp = idxjp - 1
480 END IF
481 IF( idxj.LE.nlp1 ) THEN
482 idxj = idxj - 1
483 END IF
484 givcol( givptr, 2 ) = idxjp
485 givcol( givptr, 1 ) = idxj
486 givnum( givptr, 2 ) = c
487 givnum( givptr, 1 ) = s
488 END IF
489 CALL drot( 1, vf( jprev ), 1, vf( j ), 1, c, s )
490 CALL drot( 1, vl( jprev ), 1, vl( j ), 1, c, s )
491 k2 = k2 - 1
492 idxp( k2 ) = jprev
493 jprev = j
494 ELSE
495 k = k + 1
496 zw( k ) = z( jprev )
497 dsigma( k ) = d( jprev )
498 idxp( k ) = jprev
499 jprev = j
500 END IF
501 END IF
502 GO TO 80
503 90 CONTINUE
504*
505* Record the last singular value.
506*
507 k = k + 1
508 zw( k ) = z( jprev )
509 dsigma( k ) = d( jprev )
510 idxp( k ) = jprev
511*
512 100 CONTINUE
513*
514* Sort the singular values into DSIGMA. The singular values which
515* were not deflated go into the first K slots of DSIGMA, except
516* that DSIGMA(1) is treated separately.
517*
518 DO 110 j = 2, n
519 jp = idxp( j )
520 dsigma( j ) = d( jp )
521 vfw( j ) = vf( jp )
522 vlw( j ) = vl( jp )
523 110 CONTINUE
524 IF( icompq.EQ.1 ) THEN
525 DO 120 j = 2, n
526 jp = idxp( j )
527 perm( j ) = idxq( idx( jp )+1 )
528 IF( perm( j ).LE.nlp1 ) THEN
529 perm( j ) = perm( j ) - 1
530 END IF
531 120 CONTINUE
532 END IF
533*
534* The deflated singular values go back into the last N - K slots of
535* D.
536*
537 CALL dcopy( n-k, dsigma( k+1 ), 1, d( k+1 ), 1 )
538*
539* Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
540* VL(M).
541*
542 dsigma( 1 ) = zero
543 hlftol = tol / two
544 IF( abs( dsigma( 2 ) ).LE.hlftol )
545 $ dsigma( 2 ) = hlftol
546 IF( m.GT.n ) THEN
547 z( 1 ) = dlapy2( z1, z( m ) )
548 IF( z( 1 ).LE.tol ) THEN
549 c = one
550 s = zero
551 z( 1 ) = tol
552 ELSE
553 c = z1 / z( 1 )
554 s = -z( m ) / z( 1 )
555 END IF
556 CALL drot( 1, vf( m ), 1, vf( 1 ), 1, c, s )
557 CALL drot( 1, vl( m ), 1, vl( 1 ), 1, c, s )
558 ELSE
559 IF( abs( z1 ).LE.tol ) THEN
560 z( 1 ) = tol
561 ELSE
562 z( 1 ) = z1
563 END IF
564 END IF
565*
566* Restore Z, VF, and VL.
567*
568 CALL dcopy( k-1, zw( 2 ), 1, z( 2 ), 1 )
569 CALL dcopy( n-1, vfw( 2 ), 1, vf( 2 ), 1 )
570 CALL dcopy( n-1, vlw( 2 ), 1, vl( 2 ), 1 )
571*
572 RETURN
573*
574* End of DLASD7
575*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlamrg
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:97
dlapy2
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:61
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
Here is the call graph for this function:
Here is the caller graph for this function: