[BACK]Return to dbdsvdx.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dbdsvdx.f between versions 1.3 and 1.4

version 1.3, 2016/08/27 15:34:19 version 1.4, 2017/06/17 10:53:46
Line 2 Line 2
 *  *
 *  =========== DOCUMENTATION ===========  *  =========== DOCUMENTATION ===========
 *  *
 * Online html documentation available at   * Online html documentation available at
 *            http://www.netlib.org/lapack/explore-html/   *            http://www.netlib.org/lapack/explore-html/
 *  *
 *> \htmlonly  *> \htmlonly
 *> Download DBDSVDX + dependencies   *> Download DBDSVDX + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsvdx.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsvdx.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsvdx.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsvdx.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsvdx.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsvdx.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *     SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,   *     SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
 *    $                    NS, S, Z, LDZ, WORK, IWORK, INFO )  *    $                    NS, S, Z, LDZ, WORK, IWORK, INFO )
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
Line 28 Line 28
 *     ..  *     ..
 *     .. Array Arguments ..  *     .. Array Arguments ..
 *      INTEGER            IWORK( * )  *      INTEGER            IWORK( * )
 *      DOUBLE PRECISION   D( * ), E( * ), S( * ), WORK( * ),   *      DOUBLE PRECISION   D( * ), E( * ), S( * ), WORK( * ),
 *                         Z( LDZ, * )  *                         Z( LDZ, * )
 *       ..  *       ..
 *    *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
 *>  *>
 *> \verbatim  *> \verbatim
 *>  *>
 *>  DBDSVDX computes the singular value decomposition (SVD) of a real  *>  DBDSVDX computes the singular value decomposition (SVD) of a real
 *>  N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT,   *>  N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT,
 *>  where S is a diagonal matrix with non-negative diagonal elements   *>  where S is a diagonal matrix with non-negative diagonal elements
 *>  (the singular values of B), and U and VT are orthogonal matrices   *>  (the singular values of B), and U and VT are orthogonal matrices
 *>  of left and right singular vectors, respectively.  *>  of left and right singular vectors, respectively.
 *>  *>
 *>  Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ]   *>  Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ]
 *>  and superdiagonal E = [ e_1 e_2 ... e_N-1 ], DBDSVDX computes the   *>  and superdiagonal E = [ e_1 e_2 ... e_N-1 ], DBDSVDX computes the
 *>  singular value decompositon of B through the eigenvalues and   *>  singular value decompositon of B through the eigenvalues and
 *>  eigenvectors of the N*2-by-N*2 tridiagonal matrix  *>  eigenvectors of the N*2-by-N*2 tridiagonal matrix
 *>               *>
 *>        |  0  d_1                |        *>        |  0  d_1                |
 *>        | d_1  0  e_1            |           *>        | d_1  0  e_1            |
 *>  TGK = |     e_1  0  d_2        |   *>  TGK = |     e_1  0  d_2        |
 *>        |         d_2  .   .     |        *>        |         d_2  .   .     |
 *>        |              .   .   . |  *>        |              .   .   . |
 *>  *>
 *>  If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then   *>  If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then
 *>  (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) /   *>  (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) /
 *>  sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and   *>  sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and
 *>  P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ].   *>  P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ].
 *>  *>
 *>  Given a TGK matrix, one can either a) compute -s,-v and change signs   *>  Given a TGK matrix, one can either a) compute -s,-v and change signs
 *>  so that the singular values (and corresponding vectors) are already in   *>  so that the singular values (and corresponding vectors) are already in
 *>  descending order (as in DGESVD/DGESDD) or b) compute s,v and reorder   *>  descending order (as in DGESVD/DGESDD) or b) compute s,v and reorder
 *>  the values (and corresponding vectors). DBDSVDX implements a) by   *>  the values (and corresponding vectors). DBDSVDX implements a) by
 *>  calling DSTEVX (bisection plus inverse iteration, to be replaced   *>  calling DSTEVX (bisection plus inverse iteration, to be replaced
 *>  with a version of the Multiple Relative Robust Representation   *>  with a version of the Multiple Relative Robust Representation
 *>  algorithm. (See P. Willems and B. Lang, A framework for the MR^3   *>  algorithm. (See P. Willems and B. Lang, A framework for the MR^3
 *>  algorithm: theory and implementation, SIAM J. Sci. Comput.,   *>  algorithm: theory and implementation, SIAM J. Sci. Comput.,
 *>  35:740-766, 2013.)  *>  35:740-766, 2013.)
 *> \endverbatim  *> \endverbatim
 *  *
Line 101 Line 101
 *>          N is INTEGER  *>          N is INTEGER
 *>          The order of the bidiagonal matrix.  N >= 0.  *>          The order of the bidiagonal matrix.  N >= 0.
 *> \endverbatim  *> \endverbatim
 *>   *>
 *> \param[in] D  *> \param[in] D
 *> \verbatim  *> \verbatim
 *>          D is DOUBLE PRECISION array, dimension (N)  *>          D is DOUBLE PRECISION array, dimension (N)
 *>          The n diagonal elements of the bidiagonal matrix B.  *>          The n diagonal elements of the bidiagonal matrix B.
 *> \endverbatim  *> \endverbatim
 *>   *>
 *> \param[in] E  *> \param[in] E
 *> \verbatim  *> \verbatim
 *>          E is DOUBLE PRECISION array, dimension (max(1,N-1))  *>          E is DOUBLE PRECISION array, dimension (max(1,N-1))
Line 167 Line 167
 *> \verbatim  *> \verbatim
 *>          Z is DOUBLE PRECISION array, dimension (2*N,K) )  *>          Z is DOUBLE PRECISION array, dimension (2*N,K) )
 *>          If JOBZ = 'V', then if INFO = 0 the first NS columns of Z  *>          If JOBZ = 'V', then if INFO = 0 the first NS columns of Z
 *>          contain the singular vectors of the matrix B corresponding to   *>          contain the singular vectors of the matrix B corresponding to
 *>          the selected singular values, with U in rows 1 to N and V  *>          the selected singular values, with U in rows 1 to N and V
 *>          in rows N+1 to N*2, i.e.  *>          in rows N+1 to N*2, i.e.
 *>          Z = [ U ]   *>          Z = [ U ]
 *>              [ V ]  *>              [ V ]
 *>          If JOBZ = 'N', then Z is not referenced.      *>          If JOBZ = 'N', then Z is not referenced.
 *>          Note: The user must ensure that at least K = NS+1 columns are   *>          Note: The user must ensure that at least K = NS+1 columns are
 *>          supplied in the array Z; if RANGE = 'V', the exact value of   *>          supplied in the array Z; if RANGE = 'V', the exact value of
 *>          NS is not known in advance and an upper bound must be used.  *>          NS is not known in advance and an upper bound must be used.
 *> \endverbatim  *> \endverbatim
 *>  *>
Line 194 Line 194
 *> \verbatim  *> \verbatim
 *>          IWORK is INTEGER array, dimension (12*N)  *>          IWORK is INTEGER array, dimension (12*N)
 *>          If JOBZ = 'V', then if INFO = 0, the first NS elements of  *>          If JOBZ = 'V', then if INFO = 0, the first NS elements of
 *>          IWORK are zero. If INFO > 0, then IWORK contains the indices   *>          IWORK are zero. If INFO > 0, then IWORK contains the indices
 *>          of the eigenvectors that failed to converge in DSTEVX.  *>          of the eigenvectors that failed to converge in DSTEVX.
 *> \endverbatim  *> \endverbatim
 *>  *>
Line 213 Line 213
 *  Authors:  *  Authors:
 *  ========  *  ========
 *  *
 *> \author Univ. of Tennessee   *> \author Univ. of Tennessee
 *> \author Univ. of California Berkeley   *> \author Univ. of California Berkeley
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.   *> \author NAG Ltd.
 *  *
 *> \date June 2016  *> \date June 2016
 *  *
 *> \ingroup doubleOTHEReigen  *> \ingroup doubleOTHEReigen
 *  *
 *  =====================================================================     *  =====================================================================
       SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,         SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
      $                    NS, S, Z, LDZ, WORK, IWORK, INFO)       $                    NS, S, Z, LDZ, WORK, IWORK, INFO)
 *  *
 *  -- LAPACK driver routine (version 3.6.1) --  *  -- LAPACK driver routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *     November 2016  *     December 2016
 *    *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          JOBZ, RANGE, UPLO        CHARACTER          JOBZ, RANGE, UPLO
       INTEGER            IL, INFO, IU, LDZ, N, NS        INTEGER            IL, INFO, IU, LDZ, N, NS
Line 238 Line 238
 *     ..  *     ..
 *     .. Array Arguments ..  *     .. Array Arguments ..
       INTEGER            IWORK( * )        INTEGER            IWORK( * )
       DOUBLE PRECISION   D( * ), E( * ), S( * ), WORK( * ),         DOUBLE PRECISION   D( * ), E( * ), S( * ), WORK( * ),
      $                   Z( LDZ, * )       $                   Z( LDZ, * )
 *     ..  *     ..
 *  *
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
       DOUBLE PRECISION   ZERO, ONE, TEN, HNDRD, MEIGTH         DOUBLE PRECISION   ZERO, ONE, TEN, HNDRD, MEIGTH
       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0,         PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0,
      $                     HNDRD = 100.0D0, MEIGTH = -0.1250D0 )       $                     HNDRD = 100.0D0, MEIGTH = -0.1250D0 )
       DOUBLE PRECISION   FUDGE        DOUBLE PRECISION   FUDGE
       PARAMETER          ( FUDGE = 2.0D0 )        PARAMETER          ( FUDGE = 2.0D0 )
 *     ..  *     ..
 *     .. Local Scalars ..   *     .. Local Scalars ..
       CHARACTER          RNGVX        CHARACTER          RNGVX
       LOGICAL            ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ         LOGICAL            ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ
       INTEGER            I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR,         INTEGER            I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR,
      $                   IETGK, IIFAIL, IIWORK, ILTGK, IROWU, IROWV,        $                   IETGK, IIFAIL, IIWORK, ILTGK, IROWU, IROWV,
      $                   IROWZ, ISBEG, ISPLT, ITEMP, IUTGK, J, K,        $                   IROWZ, ISBEG, ISPLT, ITEMP, IUTGK, J, K,
      $                   NTGK, NRU, NRV, NSL       $                   NTGK, NRU, NRV, NSL
       DOUBLE PRECISION   ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX,        DOUBLE PRECISION   ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX,
      $                   SMIN, SQRT2, THRESH, TOL, ULP,        $                   SMIN, SQRT2, THRESH, TOL, ULP,
      $                   VLTGK, VUTGK, ZJTJI       $                   VLTGK, VUTGK, ZJTJI
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
Line 274 Line 274
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC          ABS, DBLE, SIGN, SQRT        INTRINSIC          ABS, DBLE, SIGN, SQRT
 *     ..  *     ..
 *     .. Executable Statements ..        *     .. Executable Statements ..
 *  *
 *     Test the input parameters.  *     Test the input parameters.
 *  *
Line 321 Line 321
 *  *
       NS = 0        NS = 0
       IF( N.EQ.0 ) RETURN        IF( N.EQ.0 ) RETURN
 *       *
       IF( N.EQ.1 ) THEN        IF( N.EQ.1 ) THEN
          IF( ALLSV .OR. INDSV ) THEN           IF( ALLSV .OR. INDSV ) THEN
             NS = 1              NS = 1
Line 339 Line 339
          RETURN           RETURN
       END IF        END IF
 *  *
       ABSTOL = 2*DLAMCH( 'Safe Minimum' )         ABSTOL = 2*DLAMCH( 'Safe Minimum' )
       ULP = DLAMCH( 'Precision' )        ULP = DLAMCH( 'Precision' )
       EPS = DLAMCH( 'Epsilon' )        EPS = DLAMCH( 'Epsilon' )
       SQRT2 = SQRT( 2.0D0 )        SQRT2 = SQRT( 2.0D0 )
       ORTOL = SQRT( ULP )        ORTOL = SQRT( ULP )
 *     *
 *     Criterion for splitting is taken from DBDSQR when singular  *     Criterion for splitting is taken from DBDSQR when singular
 *     values are computed to relative accuracy TOL. (See J. Demmel and   *     values are computed to relative accuracy TOL. (See J. Demmel and
 *     W. Kahan, Accurate singular values of bidiagonal matrices, SIAM   *     W. Kahan, Accurate singular values of bidiagonal matrices, SIAM
 *     J. Sci. and Stat. Comput., 11:873–912, 1990.)  *     J. Sci. and Stat. Comput., 11:873–912, 1990.)
 *   *
       TOL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )*EPS        TOL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )*EPS
 *  *
 *     Compute approximate maximum, minimum singular values.  *     Compute approximate maximum, minimum singular values.
Line 390 Line 390
       IIWORK = IIFAIL + N*2        IIWORK = IIFAIL + N*2
 *  *
 *     Set RNGVX, which corresponds to RANGE for DSTEVX in TGK mode.  *     Set RNGVX, which corresponds to RANGE for DSTEVX in TGK mode.
 *     VL,VU or IL,IU are redefined to conform to implementation a)   *     VL,VU or IL,IU are redefined to conform to implementation a)
 *     described in the leading comments.  *     described in the leading comments.
 *  *
       ILTGK = 0        ILTGK = 0
       IUTGK = 0         IUTGK = 0
       VLTGK = ZERO        VLTGK = ZERO
       VUTGK = ZERO        VUTGK = ZERO
 *  *
       IF( ALLSV ) THEN        IF( ALLSV ) THEN
 *  *
 *        All singular values will be found. We aim at -s (see   *        All singular values will be found. We aim at -s (see
 *        leading comments) with RNGVX = 'I'. IL and IU are set  *        leading comments) with RNGVX = 'I'. IL and IU are set
 *        later (as ILTGK and IUTGK) according to the dimension   *        later (as ILTGK and IUTGK) according to the dimension
 *        of the active submatrix.  *        of the active submatrix.
 *  *
          RNGVX = 'I'           RNGVX = 'I'
Line 419 Line 419
          WORK( IDTGK:IDTGK+2*N-1 ) = ZERO           WORK( IDTGK:IDTGK+2*N-1 ) = ZERO
          CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )           CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )
          CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )           CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
          CALL DSTEVX( 'N', 'V', N*2, WORK( IDTGK ), WORK( IETGK ),            CALL DSTEVX( 'N', 'V', N*2, WORK( IDTGK ), WORK( IETGK ),
      $                VLTGK, VUTGK, ILTGK, ILTGK, ABSTOL, NS, S,       $                VLTGK, VUTGK, ILTGK, ILTGK, ABSTOL, NS, S,
      $                Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),       $                Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),
      $                IWORK( IIFAIL ), INFO )       $                IWORK( IIFAIL ), INFO )
Line 430 Line 430
          END IF           END IF
       ELSE IF( INDSV ) THEN        ELSE IF( INDSV ) THEN
 *  *
 *        Find the IL-th through the IU-th singular values. We aim   *        Find the IL-th through the IU-th singular values. We aim
 *        at -s (see leading comments) and indices are mapped into   *        at -s (see leading comments) and indices are mapped into
 *        values, therefore mimicking DSTEBZ, where  *        values, therefore mimicking DSTEBZ, where
 *  *
 *        GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN  *        GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
 *        GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN  *        GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
 *  *
          ILTGK = IL           ILTGK = IL
          IUTGK = IU            IUTGK = IU
          RNGVX = 'V'           RNGVX = 'V'
          WORK( IDTGK:IDTGK+2*N-1 ) = ZERO           WORK( IDTGK:IDTGK+2*N-1 ) = ZERO
          CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )           CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )
          CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )           CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
          CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ),            CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ),
      $                VLTGK, VLTGK, ILTGK, ILTGK, ABSTOL, NS, S,       $                VLTGK, VLTGK, ILTGK, ILTGK, ABSTOL, NS, S,
      $                Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),       $                Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),
      $                IWORK( IIFAIL ), INFO )       $                IWORK( IIFAIL ), INFO )
Line 451 Line 451
          WORK( IDTGK:IDTGK+2*N-1 ) = ZERO           WORK( IDTGK:IDTGK+2*N-1 ) = ZERO
          CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )           CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )
          CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )           CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
          CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ),            CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ),
      $                VUTGK, VUTGK, IUTGK, IUTGK, ABSTOL, NS, S,       $                VUTGK, VUTGK, IUTGK, IUTGK, ABSTOL, NS, S,
      $                Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),       $                Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),
      $                IWORK( IIFAIL ), INFO )       $                IWORK( IIFAIL ), INFO )
Line 459 Line 459
          VUTGK = MIN( VUTGK, ZERO )           VUTGK = MIN( VUTGK, ZERO )
 *  *
 *        If VLTGK=VUTGK, DSTEVX returns an error message,  *        If VLTGK=VUTGK, DSTEVX returns an error message,
 *        so if needed we change VUTGK slightly.      *        so if needed we change VUTGK slightly.
 *  *
          IF( VLTGK.EQ.VUTGK ) VLTGK = VLTGK - TOL           IF( VLTGK.EQ.VUTGK ) VLTGK = VLTGK - TOL
 *  *
          IF( WANTZ ) CALL DLASET( 'F', N*2, IU-IL+1, ZERO, ZERO, Z, LDZ)           IF( WANTZ ) CALL DLASET( 'F', N*2, IU-IL+1, ZERO, ZERO, Z, LDZ)
       END IF                     END IF
 *  *
 *     Initialize variables and pointers for S, Z, and WORK.  *     Initialize variables and pointers for S, Z, and WORK.
 *  *
Line 483 Line 483
       IROWU = 2        IROWU = 2
       IROWV = 1        IROWV = 1
       SPLIT = .FALSE.        SPLIT = .FALSE.
       SVEQ0 = .FALSE.            SVEQ0 = .FALSE.
 *  *
 *     Form the tridiagonal TGK matrix.  *     Form the tridiagonal TGK matrix.
 *  *
Line 494 Line 494
       CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )        CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
 *  *
 *  *
 *     Check for splits in two levels, outer level   *     Check for splits in two levels, outer level
 *     in E and inner level in D.  *     in E and inner level in D.
 *  *
       DO IEPTR = 2, N*2, 2         DO IEPTR = 2, N*2, 2
          IF( WORK( IETGK+IEPTR-1 ).EQ.ZERO ) THEN                  IF( WORK( IETGK+IEPTR-1 ).EQ.ZERO ) THEN
 *  *
 *           Split in E (this piece of B is square) or bottom  *           Split in E (this piece of B is square) or bottom
 *           of the (input bidiagonal) matrix.  *           of the (input bidiagonal) matrix.
 *            *
             ISPLT = IDBEG              ISPLT = IDBEG
             IDEND = IEPTR - 1              IDEND = IEPTR - 1
             DO IDPTR = IDBEG, IDEND, 2              DO IDPTR = IDBEG, IDEND, 2
Line 519 Line 519
                      IF( IDBEG.EQ.IDEND) THEN                       IF( IDBEG.EQ.IDEND) THEN
                         NRU = 1                          NRU = 1
                         NRV = 1                          NRV = 1
                      END IF                           END IF
                   ELSE IF( IDPTR.EQ.IDEND ) THEN                    ELSE IF( IDPTR.EQ.IDEND ) THEN
 *  *
 *                    D=0 at the bottom.  *                    D=0 at the bottom.
 *  *
                      SVEQ0 = .TRUE.                       SVEQ0 = .TRUE.
                      NRU = (IDEND-ISPLT)/2 + 1                        NRU = (IDEND-ISPLT)/2 + 1
                      NRV = NRU                         NRV = NRU
                      IF( ISPLT.NE.IDBEG ) THEN                       IF( ISPLT.NE.IDBEG ) THEN
                         NRU = NRU + 1                          NRU = NRU + 1
                      END IF                           END IF
                   ELSE                    ELSE
                      IF( ISPLT.EQ.IDBEG ) THEN                       IF( ISPLT.EQ.IDBEG ) THEN
 *  *
 *                       Split: top rectangular submatrix.  *                       Split: top rectangular submatrix.
 *   *
                         NRU = (IDPTR-IDBEG)/2                          NRU = (IDPTR-IDBEG)/2
                         NRV = NRU + 1                          NRV = NRU + 1
                      ELSE                       ELSE
Line 542 Line 542
 *                       Split: middle square submatrix.  *                       Split: middle square submatrix.
 *  *
                         NRU = (IDPTR-ISPLT)/2 + 1                          NRU = (IDPTR-ISPLT)/2 + 1
                         NRV = NRU                                             NRV = NRU
                      END IF                       END IF
                   END IF                    END IF
                ELSE IF( IDPTR.EQ.IDEND ) THEN                 ELSE IF( IDPTR.EQ.IDEND ) THEN
Line 560 Line 560
 *                    Split: bottom rectangular submatrix.  *                    Split: bottom rectangular submatrix.
 *  *
                      NRV = (IDEND-ISPLT)/2 + 1                       NRV = (IDEND-ISPLT)/2 + 1
                      NRU = NRV + 1                                         NRU = NRV + 1
                   END IF                    END IF
                END IF                 END IF
 *  *
Line 568 Line 568
 *  *
                IF( NTGK.GT.0 ) THEN                 IF( NTGK.GT.0 ) THEN
 *  *
 *                 Compute eigenvalues/vectors of the active   *                 Compute eigenvalues/vectors of the active
 *                 submatrix according to RANGE:   *                 submatrix according to RANGE:
 *                 if RANGE='A' (ALLSV) then RNGVX = 'I'  *                 if RANGE='A' (ALLSV) then RNGVX = 'I'
 *                 if RANGE='V' (VALSV) then RNGVX = 'V'  *                 if RANGE='V' (VALSV) then RNGVX = 'V'
 *                 if RANGE='I' (INDSV) then RNGVX = 'V'  *                 if RANGE='I' (INDSV) then RNGVX = 'V'
 *  *
                   ILTGK = 1                    ILTGK = 1
                   IUTGK = NTGK / 2                     IUTGK = NTGK / 2
                   IF( ALLSV .OR. VUTGK.EQ.ZERO ) THEN                    IF( ALLSV .OR. VUTGK.EQ.ZERO ) THEN
                      IF( SVEQ0 .OR.                        IF( SVEQ0 .OR.
      $                   SMIN.LT.EPS .OR.        $                   SMIN.LT.EPS .OR.
      $                   MOD(NTGK,2).GT.0 ) THEN       $                   MOD(NTGK,2).GT.0 ) THEN
 *                        Special case: eigenvalue equal to zero or very  *                        Special case: eigenvalue equal to zero or very
 *                        small, additional eigenvector is needed.  *                        small, additional eigenvector is needed.
                          IUTGK = IUTGK + 1                           IUTGK = IUTGK + 1
                      END IF                           END IF
                   END IF                    END IF
 *  *
 *                 Workspace needed by DSTEVX:     *                 Workspace needed by DSTEVX:
 *                 WORK( ITEMP: ): 2*5*NTGK   *                 WORK( ITEMP: ): 2*5*NTGK
 *                 IWORK( 1: ): 2*6*NTGK  *                 IWORK( 1: ): 2*6*NTGK
 *  *
                   CALL DSTEVX( JOBZ, RNGVX, NTGK, WORK( IDTGK+ISPLT-1 ),                     CALL DSTEVX( JOBZ, RNGVX, NTGK, WORK( IDTGK+ISPLT-1 ),
      $                         WORK( IETGK+ISPLT-1 ), VLTGK, VUTGK,        $                         WORK( IETGK+ISPLT-1 ), VLTGK, VUTGK,
      $                         ILTGK, IUTGK, ABSTOL, NSL, S( ISBEG ),        $                         ILTGK, IUTGK, ABSTOL, NSL, S( ISBEG ),
      $                         Z( IROWZ,ICOLZ ), LDZ, WORK( ITEMP ),        $                         Z( IROWZ,ICOLZ ), LDZ, WORK( ITEMP ),
      $                         IWORK( IIWORK ), IWORK( IIFAIL ),       $                         IWORK( IIWORK ), IWORK( IIFAIL ),
      $                         INFO )       $                         INFO )
                   IF( INFO.NE.0 ) THEN                    IF( INFO.NE.0 ) THEN
Line 601 Line 601
                      RETURN                       RETURN
                   END IF                    END IF
                   EMIN = ABS( MAXVAL( S( ISBEG:ISBEG+NSL-1 ) ) )                    EMIN = ABS( MAXVAL( S( ISBEG:ISBEG+NSL-1 ) ) )
 *     *
                   IF( NSL.GT.0 .AND. WANTZ ) THEN                    IF( NSL.GT.0 .AND. WANTZ ) THEN
 *  *
 *                    Normalize u=Z([2,4,...],:) and v=Z([1,3,...],:),  *                    Normalize u=Z([2,4,...],:) and v=Z([1,3,...],:),
Line 615 Line 615
                      IF( NSL.GT.1 .AND.                       IF( NSL.GT.1 .AND.
      $                   VUTGK.EQ.ZERO .AND.       $                   VUTGK.EQ.ZERO .AND.
      $                   MOD(NTGK,2).EQ.0 .AND.       $                   MOD(NTGK,2).EQ.0 .AND.
      $                   EMIN.EQ.0 .AND. .NOT.SPLIT ) THEN        $                   EMIN.EQ.0 .AND. .NOT.SPLIT ) THEN
 *  *
 *                       D=0 at the top or bottom of the active submatrix:  *                       D=0 at the top or bottom of the active submatrix:
 *                       one eigenvalue is equal to zero; concatenate the   *                       one eigenvalue is equal to zero; concatenate the
 *                       eigenvectors corresponding to the two smallest   *                       eigenvectors corresponding to the two smallest
 *                       eigenvalues.  *                       eigenvalues.
 *  *
                         Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) =                          Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) =
      $                  Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) +       $                  Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) +
      $                  Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 )       $                  Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 )
                         Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) =                           Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) =
      $                  ZERO        $                  ZERO
 *                       IF( IUTGK*2.GT.NTGK ) THEN  *                       IF( IUTGK*2.GT.NTGK ) THEN
 *                          Eigenvalue equal to zero or very small.  *                          Eigenvalue equal to zero or very small.
 *                          NSL = NSL - 1  *                          NSL = NSL - 1
 *                       END IF    *                       END IF
                      END IF                       END IF
 *  *
                      DO I = 0, MIN( NSL-1, NRU-1 )                       DO I = 0, MIN( NSL-1, NRU-1 )
Line 639 Line 639
                            INFO = N*2 + 1                             INFO = N*2 + 1
                            RETURN                             RETURN
                         END IF                          END IF
                         CALL DSCAL( NRU, ONE/NRMU,                           CALL DSCAL( NRU, ONE/NRMU,
      $                              Z( IROWU,ICOLZ+I ), 2 )       $                              Z( IROWU,ICOLZ+I ), 2 )
                         IF( NRMU.NE.ONE .AND.                          IF( NRMU.NE.ONE .AND.
      $                      ABS( NRMU-ORTOL )*SQRT2.GT.ONE )       $                      ABS( NRMU-ORTOL )*SQRT2.GT.ONE )
      $                      THEN       $                      THEN
                            DO J = 0, I-1                             DO J = 0, I-1
                               ZJTJI = -DDOT( NRU, Z( IROWU, ICOLZ+J ),                                 ZJTJI = -DDOT( NRU, Z( IROWU, ICOLZ+J ),
      $                                       2, Z( IROWU, ICOLZ+I ), 2 )       $                                       2, Z( IROWU, ICOLZ+I ), 2 )
                               CALL DAXPY( NRU, ZJTJI,                                 CALL DAXPY( NRU, ZJTJI,
      $                                    Z( IROWU, ICOLZ+J ), 2,       $                                    Z( IROWU, ICOLZ+J ), 2,
      $                                    Z( IROWU, ICOLZ+I ), 2 )       $                                    Z( IROWU, ICOLZ+I ), 2 )
                            END DO                             END DO
                            NRMU = DNRM2( NRU, Z( IROWU, ICOLZ+I ), 2 )                             NRMU = DNRM2( NRU, Z( IROWU, ICOLZ+I ), 2 )
                            CALL DSCAL( NRU, ONE/NRMU,                              CALL DSCAL( NRU, ONE/NRMU,
      $                                 Z( IROWU,ICOLZ+I ), 2 )       $                                 Z( IROWU,ICOLZ+I ), 2 )
                         END IF                          END IF
                      END DO                       END DO
Line 662 Line 662
                            INFO = N*2 + 1                             INFO = N*2 + 1
                            RETURN                             RETURN
                         END IF                          END IF
                         CALL DSCAL( NRV, -ONE/NRMV,                           CALL DSCAL( NRV, -ONE/NRMV,
      $                              Z( IROWV,ICOLZ+I ), 2 )       $                              Z( IROWV,ICOLZ+I ), 2 )
                         IF( NRMV.NE.ONE .AND.                          IF( NRMV.NE.ONE .AND.
      $                      ABS( NRMV-ORTOL )*SQRT2.GT.ONE )       $                      ABS( NRMV-ORTOL )*SQRT2.GT.ONE )
Line 670 Line 670
                            DO J = 0, I-1                             DO J = 0, I-1
                               ZJTJI = -DDOT( NRV, Z( IROWV, ICOLZ+J ),                                ZJTJI = -DDOT( NRV, Z( IROWV, ICOLZ+J ),
      $                                       2, Z( IROWV, ICOLZ+I ), 2 )       $                                       2, Z( IROWV, ICOLZ+I ), 2 )
                               CALL DAXPY( NRU, ZJTJI,                                 CALL DAXPY( NRU, ZJTJI,
      $                                    Z( IROWV, ICOLZ+J ), 2,       $                                    Z( IROWV, ICOLZ+J ), 2,
      $                                    Z( IROWV, ICOLZ+I ), 2 )       $                                    Z( IROWV, ICOLZ+I ), 2 )
                            END DO                             END DO
                            NRMV = DNRM2( NRV, Z( IROWV, ICOLZ+I ), 2 )                             NRMV = DNRM2( NRV, Z( IROWV, ICOLZ+I ), 2 )
                            CALL DSCAL( NRV, ONE/NRMV,                              CALL DSCAL( NRV, ONE/NRMV,
      $                                 Z( IROWV,ICOLZ+I ), 2 )       $                                 Z( IROWV,ICOLZ+I ), 2 )
                         END IF                          END IF
                      END DO                       END DO
Line 684 Line 684
      $                   MOD(NTGK,2).GT.0 ) THEN       $                   MOD(NTGK,2).GT.0 ) THEN
 *  *
 *                       D=0 in the middle of the active submatrix (one  *                       D=0 in the middle of the active submatrix (one
 *                       eigenvalue is equal to zero): save the corresponding   *                       eigenvalue is equal to zero): save the corresponding
 *                       eigenvector for later use (when bottom of the  *                       eigenvector for later use (when bottom of the
 *                       active submatrix is reached).  *                       active submatrix is reached).
 *  *
                         SPLIT = .TRUE.                          SPLIT = .TRUE.
                         Z( IROWZ:IROWZ+NTGK-1,N+1 ) =                           Z( IROWZ:IROWZ+NTGK-1,N+1 ) =
      $                     Z( IROWZ:IROWZ+NTGK-1,NS+NSL )       $                     Z( IROWZ:IROWZ+NTGK-1,NS+NSL )
                         Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) =                           Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) =
      $                     ZERO        $                     ZERO
                      END IF                                            END IF
                   END IF !** WANTZ **!                    END IF !** WANTZ **!
 *                *
                   NSL = MIN( NSL, NRU )                    NSL = MIN( NSL, NRU )
                   SVEQ0 = .FALSE.                    SVEQ0 = .FALSE.
 *  *
Line 706 Line 706
                   END DO                    END DO
 *  *
 *                 Update pointers for TGK, S and Z.  *                 Update pointers for TGK, S and Z.
 *                 *
                   ISBEG = ISBEG + NSL                    ISBEG = ISBEG + NSL
                   IROWZ = IROWZ + NTGK                    IROWZ = IROWZ + NTGK
                   ICOLZ = ICOLZ + NSL                    ICOLZ = ICOLZ + NSL
                   IROWU = IROWZ                    IROWU = IROWZ
                   IROWV = IROWZ + 1                                    IROWV = IROWZ + 1
                   ISPLT = IDPTR + 1                    ISPLT = IDPTR + 1
                   NS = NS + NSL                    NS = NS + NSL
                   NRU = 0                    NRU = 0
                   NRV = 0                           NRV = 0
                END IF !** NTGK.GT.0 **!                  END IF !** NTGK.GT.0 **!
                IF( IROWZ.LT.N*2 .AND. WANTZ ) THEN                 IF( IROWZ.LT.N*2 .AND. WANTZ ) THEN
                   Z( 1:IROWZ-1, ICOLZ ) = ZERO                    Z( 1:IROWZ-1, ICOLZ ) = ZERO
                END IF                 END IF
Line 726 Line 726
 *              Bring back eigenvector corresponding  *              Bring back eigenvector corresponding
 *              to eigenvalue equal to zero.  *              to eigenvalue equal to zero.
 *  *
                Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) =                  Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) =
      $         Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) +       $         Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) +
      $         Z( IDBEG:IDEND-NTGK+1,N+1 )       $         Z( IDBEG:IDEND-NTGK+1,N+1 )
                Z( IDBEG:IDEND-NTGK+1,N+1 ) = 0                 Z( IDBEG:IDEND-NTGK+1,N+1 ) = 0
Line 735 Line 735
             IROWU = IROWU + 1              IROWU = IROWU + 1
             IDBEG = IEPTR + 1              IDBEG = IEPTR + 1
             SVEQ0 = .FALSE.              SVEQ0 = .FALSE.
             SPLIT = .FALSE.                      SPLIT = .FALSE.
          END IF !** Check for split in E **!           END IF !** Check for split in E **!
       END DO !** IEPTR loop **!        END DO !** IEPTR loop **!
 *  *
Line 757 Line 757
             IF( WANTZ ) CALL DSWAP( N*2, Z( 1,K ), 1, Z( 1,NS+1-I ), 1 )              IF( WANTZ ) CALL DSWAP( N*2, Z( 1,K ), 1, Z( 1,NS+1-I ), 1 )
          END IF           END IF
       END DO        END DO
 *     *
 *     If RANGE=I, check for singular values/vectors to be discarded.  *     If RANGE=I, check for singular values/vectors to be discarded.
 *  *
       IF( INDSV ) THEN        IF( INDSV ) THEN
Line 767 Line 767
             IF( WANTZ ) Z( 1:N*2,K+1:NS ) = ZERO              IF( WANTZ ) Z( 1:N*2,K+1:NS ) = ZERO
             NS = K              NS = K
          END IF           END IF
       END IF         END IF
 *  *
 *     Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ).  *     Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ).
 *     If B is a lower diagonal, swap U and V.  *     If B is a lower diagonal, swap U and V.
Removed from v.1.3  
changed lines
  Added in v.1.4

CVSweb interface <joel.bertrand@systella.fr>