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

Diff for /rpl/lapack/lapack/zstemr.f between versions 1.8 and 1.18

version 1.8, 2011/11/21 20:43:21 version 1.18, 2017/06/17 10:54:28
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 ZSTEMR + dependencies   *> Download ZSTEMR + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstemr.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstemr.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstemr.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstemr.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstemr.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstemr.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
Line 21 Line 21
 *       SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,  *       SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
 *                          M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,  *                          M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
 *                          IWORK, LIWORK, INFO )  *                          IWORK, LIWORK, INFO )
 *   *
 *       .. Scalar Arguments ..  *       .. Scalar Arguments ..
 *       CHARACTER          JOBZ, RANGE  *       CHARACTER          JOBZ, RANGE
 *       LOGICAL            TRYRAC  *       LOGICAL            TRYRAC
Line 33 Line 33
 *       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )  *       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
 *       COMPLEX*16         Z( LDZ, * )  *       COMPLEX*16         Z( LDZ, * )
 *       ..  *       ..
 *    *
 *  *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
Line 153 Line 153
 *> \param[in] VL  *> \param[in] VL
 *> \verbatim  *> \verbatim
 *>          VL is DOUBLE PRECISION  *>          VL is DOUBLE PRECISION
   *>
   *>          If RANGE='V', the lower bound of the interval to
   *>          be searched for eigenvalues. VL < VU.
   *>          Not referenced if RANGE = 'A' or 'I'.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] VU  *> \param[in] VU
 *> \verbatim  *> \verbatim
 *>          VU is DOUBLE PRECISION  *>          VU is DOUBLE PRECISION
 *>  *>
 *>          If RANGE='V', the lower and upper bounds of the interval to  *>          If RANGE='V', the upper bound of the interval to
 *>          be searched for eigenvalues. VL < VU.  *>          be searched for eigenvalues. VL < VU.
 *>          Not referenced if RANGE = 'A' or 'I'.  *>          Not referenced if RANGE = 'A' or 'I'.
 *> \endverbatim  *> \endverbatim
Line 167 Line 171
 *> \param[in] IL  *> \param[in] IL
 *> \verbatim  *> \verbatim
 *>          IL is INTEGER  *>          IL is INTEGER
   *>
   *>          If RANGE='I', the index of the
   *>          smallest eigenvalue to be returned.
   *>          1 <= IL <= IU <= N, if N > 0.
   *>          Not referenced if RANGE = 'A' or 'V'.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] IU  *> \param[in] IU
 *> \verbatim  *> \verbatim
 *>          IU is INTEGER  *>          IU is INTEGER
 *>  *>
 *>          If RANGE='I', the indices (in ascending order) of the  *>          If RANGE='I', the index of the
 *>          smallest and largest eigenvalues to be returned.  *>          largest eigenvalue to be returned.
 *>          1 <= IL <= IU <= N, if N > 0.  *>          1 <= IL <= IU <= N, if N > 0.
 *>          Not referenced if RANGE = 'A' or 'V'.  *>          Not referenced if RANGE = 'A' or 'V'.
 *> \endverbatim  *> \endverbatim
Line 306 Line 315
 *  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 November 2011  *> \date June 2016
 *  *
 *> \ingroup complex16OTHERcomputational  *> \ingroup complex16OTHERcomputational
 *  *
Line 329 Line 338
      $                   M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,       $                   M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
      $                   IWORK, LIWORK, INFO )       $                   IWORK, LIWORK, INFO )
 *  *
 *  -- LAPACK computational routine (version 3.4.0) --  *  -- LAPACK computational 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 2011  *     June 2016
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          JOBZ, RANGE        CHARACTER          JOBZ, RANGE
Line 408 Line 417
       WU = ZERO        WU = ZERO
       IIL = 0        IIL = 0
       IIU = 0        IIU = 0
         NSPLIT = 0
   
       IF( VALEIG ) THEN        IF( VALEIG ) THEN
 *        We do not reference VL, VU in the cases RANGE = 'I','A'  *        We do not reference VL, VU in the cases RANGE = 'I','A'
Line 525 Line 535
                IF (SN.NE.ZERO) THEN                 IF (SN.NE.ZERO) THEN
                   IF (CS.NE.ZERO) THEN                    IF (CS.NE.ZERO) THEN
                      ISUPPZ(2*M-1) = 1                       ISUPPZ(2*M-1) = 1
                      ISUPPZ(2*M-1) = 2                       ISUPPZ(2*M) = 2
                   ELSE                    ELSE
                      ISUPPZ(2*M-1) = 1                       ISUPPZ(2*M-1) = 1
                      ISUPPZ(2*M-1) = 1                       ISUPPZ(2*M) = 1
                   END IF                    END IF
                ELSE                 ELSE
                   ISUPPZ(2*M-1) = 2                    ISUPPZ(2*M-1) = 2
Line 549 Line 559
                IF (SN.NE.ZERO) THEN                 IF (SN.NE.ZERO) THEN
                   IF (CS.NE.ZERO) THEN                    IF (CS.NE.ZERO) THEN
                      ISUPPZ(2*M-1) = 1                       ISUPPZ(2*M-1) = 1
                      ISUPPZ(2*M-1) = 2                       ISUPPZ(2*M) = 2
                   ELSE                    ELSE
                      ISUPPZ(2*M-1) = 1                       ISUPPZ(2*M-1) = 1
                      ISUPPZ(2*M-1) = 1                       ISUPPZ(2*M) = 1
                   END IF                    END IF
                ELSE                 ELSE
                   ISUPPZ(2*M-1) = 2                    ISUPPZ(2*M-1) = 2
Line 560 Line 570
                END IF                 END IF
             ENDIF              ENDIF
          ENDIF           ENDIF
          RETURN        ELSE
       END IF  
   
 *     Continue with general N  *        Continue with general N
   
       INDGRS = 1           INDGRS = 1
       INDERR = 2*N + 1           INDERR = 2*N + 1
       INDGP = 3*N + 1           INDGP = 3*N + 1
       INDD = 4*N + 1           INDD = 4*N + 1
       INDE2 = 5*N + 1           INDE2 = 5*N + 1
       INDWRK = 6*N + 1           INDWRK = 6*N + 1
 *  *
       IINSPL = 1           IINSPL = 1
       IINDBL = N + 1           IINDBL = N + 1
       IINDW = 2*N + 1           IINDW = 2*N + 1
       IINDWK = 3*N + 1           IINDWK = 3*N + 1
 *  *
 *     Scale matrix to allowable range, if necessary.  *        Scale matrix to allowable range, if necessary.
 *     The allowable range is related to the PIVMIN parameter; see the  *        The allowable range is related to the PIVMIN parameter; see the
 *     comments in DLARRD.  The preference for scaling small values  *        comments in DLARRD.  The preference for scaling small values
 *     up is heuristic; we expect users' matrices not to be close to the  *        up is heuristic; we expect users' matrices not to be close to the
 *     RMAX threshold.  *        RMAX threshold.
 *  *
       SCALE = ONE           SCALE = ONE
       TNRM = DLANST( 'M', N, D, E )           TNRM = DLANST( 'M', N, D, E )
       IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN           IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
          SCALE = RMIN / TNRM              SCALE = RMIN / TNRM
       ELSE IF( TNRM.GT.RMAX ) THEN           ELSE IF( TNRM.GT.RMAX ) THEN
          SCALE = RMAX / TNRM              SCALE = RMAX / TNRM
       END IF           END IF
       IF( SCALE.NE.ONE ) THEN           IF( SCALE.NE.ONE ) THEN
          CALL DSCAL( N, SCALE, D, 1 )              CALL DSCAL( N, SCALE, D, 1 )
          CALL DSCAL( N-1, SCALE, E, 1 )              CALL DSCAL( N-1, SCALE, E, 1 )
          TNRM = TNRM*SCALE              TNRM = TNRM*SCALE
          IF( VALEIG ) THEN              IF( VALEIG ) THEN
 *           If eigenvalues in interval have to be found,  *              If eigenvalues in interval have to be found,
 *           scale (WL, WU] accordingly  *              scale (WL, WU] accordingly
             WL = WL*SCALE                 WL = WL*SCALE
             WU = WU*SCALE                 WU = WU*SCALE
          ENDIF              ENDIF
       END IF           END IF
 *  *
 *     Compute the desired eigenvalues of the tridiagonal after splitting  *        Compute the desired eigenvalues of the tridiagonal after splitting
 *     into smaller subblocks if the corresponding off-diagonal elements  *        into smaller subblocks if the corresponding off-diagonal elements
 *     are small  *        are small
 *     THRESH is the splitting parameter for DLARRE  *        THRESH is the splitting parameter for DLARRE
 *     A negative THRESH forces the old splitting criterion based on the  *        A negative THRESH forces the old splitting criterion based on the
 *     size of the off-diagonal. A positive THRESH switches to splitting  *        size of the off-diagonal. A positive THRESH switches to splitting
 *     which preserves relative accuracy.  *        which preserves relative accuracy.
 *  *
       IF( TRYRAC ) THEN           IF( TRYRAC ) THEN
 *        Test whether the matrix warrants the more expensive relative approach.  *           Test whether the matrix warrants the more expensive relative approach.
          CALL DLARRR( N, D, E, IINFO )              CALL DLARRR( N, D, E, IINFO )
       ELSE           ELSE
 *        The user does not care about relative accurately eigenvalues  *           The user does not care about relative accurately eigenvalues
          IINFO = -1              IINFO = -1
       ENDIF           ENDIF
 *     Set the splitting criterion  *        Set the splitting criterion
       IF (IINFO.EQ.0) THEN           IF (IINFO.EQ.0) THEN
          THRESH = EPS              THRESH = EPS
       ELSE           ELSE
          THRESH = -EPS              THRESH = -EPS
 *        relative accuracy is desired but T does not guarantee it  *           relative accuracy is desired but T does not guarantee it
          TRYRAC = .FALSE.              TRYRAC = .FALSE.
       ENDIF           ENDIF
 *  *
       IF( TRYRAC ) THEN           IF( TRYRAC ) THEN
 *        Copy original diagonal, needed to guarantee relative accuracy  *           Copy original diagonal, needed to guarantee relative accuracy
          CALL DCOPY(N,D,1,WORK(INDD),1)              CALL DCOPY(N,D,1,WORK(INDD),1)
       ENDIF           ENDIF
 *     Store the squares of the offdiagonal values of T  *        Store the squares of the offdiagonal values of T
       DO 5 J = 1, N-1           DO 5 J = 1, N-1
          WORK( INDE2+J-1 ) = E(J)**2              WORK( INDE2+J-1 ) = E(J)**2
  5    CONTINUE   5    CONTINUE
   
 *     Set the tolerance parameters for bisection  *        Set the tolerance parameters for bisection
       IF( .NOT.WANTZ ) THEN           IF( .NOT.WANTZ ) THEN
 *        DLARRE computes the eigenvalues to full precision.  *           DLARRE computes the eigenvalues to full precision.
          RTOL1 = FOUR * EPS              RTOL1 = FOUR * EPS
          RTOL2 = FOUR * EPS              RTOL2 = FOUR * EPS
       ELSE           ELSE
 *        DLARRE computes the eigenvalues to less than full precision.  *           DLARRE computes the eigenvalues to less than full precision.
 *        ZLARRV will refine the eigenvalue approximations, and we only  *           ZLARRV will refine the eigenvalue approximations, and we only
 *        need less accurate initial bisection in DLARRE.  *           need less accurate initial bisection in DLARRE.
 *        Note: these settings do only affect the subset case and DLARRE  *           Note: these settings do only affect the subset case and DLARRE
          RTOL1 = SQRT(EPS)              RTOL1 = SQRT(EPS)
          RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS )              RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS )
       ENDIF           ENDIF
       CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,           CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,
      $             WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT,       $             WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT,
      $             IWORK( IINSPL ), M, W, WORK( INDERR ),       $             IWORK( IINSPL ), M, W, WORK( INDERR ),
      $             WORK( INDGP ), IWORK( IINDBL ),       $             WORK( INDGP ), IWORK( IINDBL ),
      $             IWORK( IINDW ), WORK( INDGRS ), PIVMIN,       $             IWORK( IINDW ), WORK( INDGRS ), PIVMIN,
      $             WORK( INDWRK ), IWORK( IINDWK ), IINFO )       $             WORK( INDWRK ), IWORK( IINDWK ), IINFO )
       IF( IINFO.NE.0 ) THEN           IF( IINFO.NE.0 ) THEN
          INFO = 10 + ABS( IINFO )              INFO = 10 + ABS( IINFO )
          RETURN              RETURN
       END IF           END IF
 *     Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired  *        Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
 *     part of the spectrum. All desired eigenvalues are contained in  *        part of the spectrum. All desired eigenvalues are contained in
 *     (WL,WU]  *        (WL,WU]
   
   
       IF( WANTZ ) THEN           IF( WANTZ ) THEN
 *  *
 *        Compute the desired eigenvectors corresponding to the computed  *           Compute the desired eigenvectors corresponding to the computed
 *        eigenvalues  *           eigenvalues
 *  *
          CALL ZLARRV( N, WL, WU, D, E,              CALL ZLARRV( N, WL, WU, D, E,
      $                PIVMIN, IWORK( IINSPL ), M,       $                PIVMIN, IWORK( IINSPL ), M,
      $                1, M, MINRGP, RTOL1, RTOL2,       $                1, M, MINRGP, RTOL1, RTOL2,
      $                W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ),       $                W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ),
      $                IWORK( IINDW ), WORK( INDGRS ), Z, LDZ,       $                IWORK( IINDW ), WORK( INDGRS ), Z, LDZ,
      $                ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO )       $                ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO )
          IF( IINFO.NE.0 ) THEN              IF( IINFO.NE.0 ) THEN
             INFO = 20 + ABS( IINFO )                 INFO = 20 + ABS( IINFO )
             RETURN                 RETURN
          END IF              END IF
       ELSE           ELSE
 *        DLARRE computes eigenvalues of the (shifted) root representation  *           DLARRE computes eigenvalues of the (shifted) root representation
 *        ZLARRV returns the eigenvalues of the unshifted matrix.  *           ZLARRV returns the eigenvalues of the unshifted matrix.
 *        However, if the eigenvectors are not desired by the user, we need  *           However, if the eigenvectors are not desired by the user, we need
 *        to apply the corresponding shifts from DLARRE to obtain the  *           to apply the corresponding shifts from DLARRE to obtain the
 *        eigenvalues of the original matrix.  *           eigenvalues of the original matrix.
          DO 20 J = 1, M              DO 20 J = 1, M
             ITMP = IWORK( IINDBL+J-1 )                 ITMP = IWORK( IINDBL+J-1 )
             W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )                 W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
  20      CONTINUE   20      CONTINUE
       END IF           END IF
 *  *
   
       IF ( TRYRAC ) THEN           IF ( TRYRAC ) THEN
 *        Refine computed eigenvalues so that they are relatively accurate  *           Refine computed eigenvalues so that they are relatively accurate
 *        with respect to the original matrix T.  *           with respect to the original matrix T.
          IBEGIN = 1              IBEGIN = 1
          WBEGIN = 1              WBEGIN = 1
          DO 39  JBLK = 1, IWORK( IINDBL+M-1 )              DO 39  JBLK = 1, IWORK( IINDBL+M-1 )
             IEND = IWORK( IINSPL+JBLK-1 )                 IEND = IWORK( IINSPL+JBLK-1 )
             IN = IEND - IBEGIN + 1                 IN = IEND - IBEGIN + 1
             WEND = WBEGIN - 1                 WEND = WBEGIN - 1
 *           check if any eigenvalues have to be refined in this block  *              check if any eigenvalues have to be refined in this block
  36         CONTINUE   36         CONTINUE
             IF( WEND.LT.M ) THEN                 IF( WEND.LT.M ) THEN
                IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN                    IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
                   WEND = WEND + 1                       WEND = WEND + 1
                   GO TO 36                       GO TO 36
                     END IF
                  END IF
                  IF( WEND.LT.WBEGIN ) THEN
                     IBEGIN = IEND + 1
                     GO TO 39
                END IF                 END IF
             END IF  
             IF( WEND.LT.WBEGIN ) THEN  
                IBEGIN = IEND + 1  
                GO TO 39  
             END IF  
   
             OFFSET = IWORK(IINDW+WBEGIN-1)-1                 OFFSET = IWORK(IINDW+WBEGIN-1)-1
             IFIRST = IWORK(IINDW+WBEGIN-1)                 IFIRST = IWORK(IINDW+WBEGIN-1)
             ILAST = IWORK(IINDW+WEND-1)                 ILAST = IWORK(IINDW+WEND-1)
             RTOL2 = FOUR * EPS                 RTOL2 = FOUR * EPS
             CALL DLARRJ( IN,                 CALL DLARRJ( IN,
      $                   WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1),       $                   WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1),
      $                   IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN),       $                   IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN),
      $                   WORK( INDERR+WBEGIN-1 ),       $                   WORK( INDERR+WBEGIN-1 ),
      $                   WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,       $                   WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,
      $                   TNRM, IINFO )       $                   TNRM, IINFO )
             IBEGIN = IEND + 1                 IBEGIN = IEND + 1
             WBEGIN = WEND + 1                 WBEGIN = WEND + 1
  39      CONTINUE   39      CONTINUE
       ENDIF           ENDIF
 *  *
 *     If matrix was scaled, then rescale eigenvalues appropriately.  *        If matrix was scaled, then rescale eigenvalues appropriately.
 *  *
       IF( SCALE.NE.ONE ) THEN           IF( SCALE.NE.ONE ) THEN
          CALL DSCAL( M, ONE / SCALE, W, 1 )              CALL DSCAL( M, ONE / SCALE, W, 1 )
            END IF
       END IF        END IF
 *  *
 *     If eigenvalues are not in increasing order, then sort them,  *     If eigenvalues are not in increasing order, then sort them,
 *     possibly along with eigenvectors.  *     possibly along with eigenvectors.
 *  *
       IF( NSPLIT.GT.1 ) THEN        IF( NSPLIT.GT.1 .OR. N.EQ.2 ) THEN
          IF( .NOT. WANTZ ) THEN           IF( .NOT. WANTZ ) THEN
             CALL DLASRT( 'I', M, W, IINFO )              CALL DLASRT( 'I', M, W, IINFO )
             IF( IINFO.NE.0 ) THEN              IF( IINFO.NE.0 ) THEN
Removed from v.1.8  
changed lines
  Added in v.1.18

CVSweb interface <joel.bertrand@systella.fr>