LAPACK 3.12.0
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 276 of file dlasd7.f.

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*
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:99
dlapy2
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:63
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: