[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.10 and 1.11

version 1.10, 2012/08/22 09:48:40 version 1.11, 2012/12/14 12:30:34
Line 311 Line 311
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver 
 *> \author NAG Ltd.   *> \author NAG Ltd. 
 *  *
 *> \date November 2011  *> \date September 2012
 *  *
 *> \ingroup complex16OTHERcomputational  *> \ingroup complex16OTHERcomputational
 *  *
Line 329 Line 329
      $                   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.4.2) --
 *  -- 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  *     September 2012
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          JOBZ, RANGE        CHARACTER          JOBZ, RANGE
Line 560 Line 560
                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.10  
changed lines
  Added in v.1.11

CVSweb interface <joel.bertrand@systella.fr>