Chimie Théorique » OpenPath

      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 )

*

*  -- LAPACK auxiliary routine (version 3.2) --

*  -- LAPACK is a software package provided by Univ. of Tennessee,    --

*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--

*     November 2006

*

*     .. Scalar Arguments ..

      INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,

     $                   NR, SQRE

      DOUBLE PRECISION   ALPHA, BETA, C, S

*     ..

*     .. Array Arguments ..

      INTEGER            GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),

     $                   IDXQ( * ), PERM( * )

      DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),

     $                   VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),

     $                   ZW( * )

*     ..

*

*  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.

*

*  Arguments

*  =========

*

*  ICOMPQ  (input) 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.

*

*  NL     (input) INTEGER

*         The row dimension of the upper block. NL >= 1.

*

*  NR     (input) INTEGER

*         The row dimension of the lower block. NR >= 1.

*

*  SQRE   (input) 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.

*

*  K      (output) INTEGER

*         Contains the dimension of the non-deflated matrix, this is

*         the order of the related secular equation. 1 <= K <=N.

*

*  D      (input/output) 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.

*

*  Z      (output) DOUBLE PRECISION array, dimension ( M )

*         On exit Z contains the updating row vector in the secular

*         equation.

*

*  ZW     (workspace) DOUBLE PRECISION array, dimension ( M )

*         Workspace for Z.

*

*  VF     (input/output) 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.

*

*  VFW    (workspace) DOUBLE PRECISION array, dimension ( M )

*         Workspace for VF.

*

*  VL     (input/output) 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.

*

*  VLW    (workspace) DOUBLE PRECISION array, dimension ( M )

*         Workspace for VL.

*

*  ALPHA  (input) DOUBLE PRECISION

*         Contains the diagonal element associated with the added row.

*

*  BETA   (input) DOUBLE PRECISION

*         Contains the off-diagonal element associated with the added

*         row.

*

*  DSIGMA (output) DOUBLE PRECISION array, dimension ( N )

*         Contains a copy of the diagonal elements (K-1 singular values

*         and one zero) in the secular equation.

*

*  IDX    (workspace) INTEGER array, dimension ( N )

*         This will contain the permutation used to sort the contents of

*         D into ascending order.

*

*  IDXP   (workspace) 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.

*

*  IDXQ   (input) 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.

*

*  PERM   (output) INTEGER array, dimension ( N )

*         The permutations (from deflation and sorting) to be applied

*         to each singular block. Not referenced if ICOMPQ = 0.

*

*  GIVPTR (output) INTEGER

*         The number of Givens rotations which took place in this

*         subproblem. Not referenced if ICOMPQ = 0.

*

*  GIVCOL (output) 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.

*

*  LDGCOL (input) INTEGER

*         The leading dimension of GIVCOL, must be at least N.

*

*  GIVNUM (output) 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.

*

*  LDGNUM (input) INTEGER

*         The leading dimension of GIVNUM, must be at least N.

*

*  C      (output) 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.

*

*  S      (output) 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.

*

*  INFO   (output) INTEGER

*         = 0:  successful exit.

*         < 0:  if INFO = -i, the i-th argument had an illegal value.

*

*  Further Details

*  ===============

*

*  Based on contributions by

*     Ming Gu and Huan Ren, Computer Science Division, University of

*     California at Berkeley, USA

*

*  =====================================================================

*

*     .. Parameters ..

      DOUBLE PRECISION   ZERO, ONE, TWO, EIGHT

      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,

     $                   EIGHT = 8.0D+0 )

*     ..

*     .. Local Scalars ..

*

      INTEGER            I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,

     $                   NLP1, NLP2

      DOUBLE PRECISION   EPS, HLFTOL, TAU, TOL, Z1

*     ..

*     .. External Subroutines ..

      EXTERNAL           DCOPY, DLAMRG, DROT, XERBLA

*     ..

*     .. External Functions ..

      DOUBLE PRECISION   DLAMCH, DLAPY2

      EXTERNAL           DLAMCH, DLAPY2

*     ..

*     .. Intrinsic Functions ..

      INTRINSIC          ABS, MAX

*     ..

*     .. Executable Statements ..

*

*     Test the input parameters.

*

      INFO = 0

      N = NL + NR + 1

      M = N + SQRE

*

      IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN

         INFO = -1

      ELSE IF( NL.LT.1 ) THEN

         INFO = -2

      ELSE IF( NR.LT.1 ) THEN

         INFO = -3

      ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN

         INFO = -4

      ELSE IF( LDGCOL.LT.N ) THEN

         INFO = -22

      ELSE IF( LDGNUM.LT.N ) THEN

         INFO = -24

      END IF

      IF( INFO.NE.0 ) THEN

         CALL XERBLA( 'DLASD7', -INFO )

         RETURN

      END IF

*

      NLP1 = NL + 1

      NLP2 = NL + 2

      IF( ICOMPQ.EQ.1 ) THEN

         GIVPTR = 0

      END IF

*

*     Generate the first part of the vector Z and move the singular

*     values in the first part of D one position backward.

*

      Z1 = ALPHA*VL( NLP1 )

      VL( NLP1 ) = ZERO

      TAU = VF( NLP1 )

      DO 10 I = NL, 1, -1

         Z( I+1 ) = ALPHA*VL( I )

         VL( I ) = ZERO

         VF( I+1 ) = VF( I )

         D( I+1 ) = D( I )

         IDXQ( I+1 ) = IDXQ( I ) + 1

   10 CONTINUE

      VF( 1 ) = TAU

*

*     Generate the second part of the vector Z.

*

      DO 20 I = NLP2, M

         Z( I ) = BETA*VF( I )

         VF( I ) = ZERO

   20 CONTINUE

*

*     Sort the singular values into increasing order

*

      DO 30 I = NLP2, N

         IDXQ( I ) = IDXQ( I ) + NLP1

   30 CONTINUE

*

*     DSIGMA, IDXC, IDXC, and ZW are used as storage space.

*

      DO 40 I = 2, N

         DSIGMA( I ) = D( IDXQ( I ) )

         ZW( I ) = Z( IDXQ( I ) )

         VFW( I ) = VF( IDXQ( I ) )

         VLW( I ) = VL( IDXQ( I ) )

   40 CONTINUE

*

      CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )

*

      DO 50 I = 2, N

         IDXI = 1 + IDX( I )

         D( I ) = DSIGMA( IDXI )

         Z( I ) = ZW( IDXI )

         VF( I ) = VFW( IDXI )

         VL( I ) = VLW( IDXI )

   50 CONTINUE

*

*     Calculate the allowable deflation tolerence

*

      EPS = DLAMCH( 'Epsilon' )

      TOL = MAX( ABS( ALPHA ), ABS( BETA ) )

      TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL )

*

*     There are 2 kinds of deflation -- first a value in the z-vector

*     is small, second two (or more) singular values are very close

*     together (their difference is small).

*

*     If the value in the z-vector is small, we simply permute the

*     array so that the corresponding singular value is moved to the

*     end.

*

*     If two values in the D-vector are close, we perform a two-sided

*     rotation designed to make one of the corresponding z-vector

*     entries zero, and then permute the array so that the deflated

*     singular value is moved to the end.

*

*     If there are multiple singular values then the problem deflates.

*     Here the number of equal singular values are found.  As each equal

*     singular value is found, an elementary reflector is computed to

*     rotate the corresponding singular subspace so that the

*     corresponding components of Z are zero in this new basis.

*

      K = 1

      K2 = N + 1

      DO 60 J = 2, N

         IF( ABS( Z( J ) ).LE.TOL ) THEN

*

*           Deflate due to small z component.

*

            K2 = K2 - 1

            IDXP( K2 ) = J

            IF( J.EQ.N )

     $         GO TO 100

         ELSE

            JPREV = J

            GO TO 70

         END IF

   60 CONTINUE

   70 CONTINUE

      J = JPREV

   80 CONTINUE

      J = J + 1

      IF( J.GT.N )

     $   GO TO 90

      IF( ABS( Z( J ) ).LE.TOL ) THEN

*

*        Deflate due to small z component.

*

         K2 = K2 - 1

         IDXP( K2 ) = J

      ELSE

*

*        Check if singular values are close enough to allow deflation.

*

         IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN

*

*           Deflation is possible.

*

            S = Z( JPREV )

            C = Z( J )

*

*           Find sqrt(a**2+b**2) without overflow or

*           destructive underflow.

*

            TAU = DLAPY2( C, S )

            Z( J ) = TAU

            Z( JPREV ) = ZERO

            C = C / TAU

            S = -S / TAU

*

*           Record the appropriate Givens rotation

*

            IF( ICOMPQ.EQ.1 ) THEN

               GIVPTR = GIVPTR + 1

               IDXJP = IDXQ( IDX( JPREV )+1 )

               IDXJ = IDXQ( IDX( J )+1 )

               IF( IDXJP.LE.NLP1 ) THEN

                  IDXJP = IDXJP - 1

               END IF

               IF( IDXJ.LE.NLP1 ) THEN

                  IDXJ = IDXJ - 1

               END IF

               GIVCOL( GIVPTR, 2 ) = IDXJP

               GIVCOL( GIVPTR, 1 ) = IDXJ

               GIVNUM( GIVPTR, 2 ) = C

               GIVNUM( GIVPTR, 1 ) = S

            END IF

            CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )

            CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )

            K2 = K2 - 1

            IDXP( K2 ) = JPREV

            JPREV = J

         ELSE

            K = K + 1

            ZW( K ) = Z( JPREV )

            DSIGMA( K ) = D( JPREV )

            IDXP( K ) = JPREV

            JPREV = J

         END IF

      END IF

      GO TO 80

   90 CONTINUE

*

*     Record the last singular value.

*

      K = K + 1

      ZW( K ) = Z( JPREV )

      DSIGMA( K ) = D( JPREV )

      IDXP( K ) = JPREV

*

  100 CONTINUE

*

*     Sort the singular values into DSIGMA. The singular values which

*     were not deflated go into the first K slots of DSIGMA, except

*     that DSIGMA(1) is treated separately.

*

      DO 110 J = 2, N

         JP = IDXP( J )

         DSIGMA( J ) = D( JP )

         VFW( J ) = VF( JP )

         VLW( J ) = VL( JP )

  110 CONTINUE

      IF( ICOMPQ.EQ.1 ) THEN

         DO 120 J = 2, N

            JP = IDXP( J )

            PERM( J ) = IDXQ( IDX( JP )+1 )

            IF( PERM( J ).LE.NLP1 ) THEN

               PERM( J ) = PERM( J ) - 1

            END IF

  120    CONTINUE

      END IF

*

*     The deflated singular values go back into the last N - K slots of

*     D.

*

      CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )

*

*     Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and

*     VL(M).

*

      DSIGMA( 1 ) = ZERO

      HLFTOL = TOL / TWO

      IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )

     $   DSIGMA( 2 ) = HLFTOL

      IF( M.GT.N ) THEN

         Z( 1 ) = DLAPY2( Z1, Z( M ) )

         IF( Z( 1 ).LE.TOL ) THEN

            C = ONE

            S = ZERO

            Z( 1 ) = TOL

         ELSE

            C = Z1 / Z( 1 )

            S = -Z( M ) / Z( 1 )

         END IF

         CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )

         CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )

      ELSE

         IF( ABS( Z1 ).LE.TOL ) THEN

            Z( 1 ) = TOL

         ELSE

            Z( 1 ) = Z1

         END IF

      END IF

*

*     Restore Z, VF, and VL.

*

      CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )

      CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )

      CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )

*

      RETURN

*

*     End of DLASD7

*

      END