zlarrv.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 881 行 · 第 1/5 页

</span><span class="comment">*</span><span class="comment">        The eigenvalue approximations will be refined when necessary as
</span><span class="comment">*</span><span class="comment">        high relative accuracy is required for the computation of the
</span><span class="comment">*</span><span class="comment">        corresponding eigenvectors.
</span>         CALL DCOPY( IM, W( WBEGIN ), 1,
     &amp;                   WORK( WBEGIN ), 1 )

<span class="comment">*</span><span class="comment">        We store in W the eigenvalue approximations w.r.t. the original
</span><span class="comment">*</span><span class="comment">        matrix T.
</span>         DO 30 I=1,IM
            W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA
 30      CONTINUE


<span class="comment">*</span><span class="comment">        NDEPTH is the current depth of the representation tree
</span>         NDEPTH = 0
<span class="comment">*</span><span class="comment">        PARITY is either 1 or 0
</span>         PARITY = 1
<span class="comment">*</span><span class="comment">        NCLUS is the number of clusters for the next level of the
</span><span class="comment">*</span><span class="comment">        representation tree, we start with NCLUS = 1 for the root
</span>         NCLUS = 1
         IWORK( IINDC1+1 ) = 1
         IWORK( IINDC1+2 ) = IM

<span class="comment">*</span><span class="comment">        IDONE is the number of eigenvectors already computed in the current
</span><span class="comment">*</span><span class="comment">        block
</span>         IDONE = 0
<span class="comment">*</span><span class="comment">        loop while( IDONE.LT.IM )
</span><span class="comment">*</span><span class="comment">        generate the representation tree for the current block and
</span><span class="comment">*</span><span class="comment">        compute the eigenvectors
</span>   40    CONTINUE
         IF( IDONE.LT.IM ) THEN
<span class="comment">*</span><span class="comment">           This is a crude protection against infinitely deep trees
</span>            IF( NDEPTH.GT.M ) THEN
               INFO = -2
               RETURN
            ENDIF
<span class="comment">*</span><span class="comment">           breadth first processing of the current level of the representation
</span><span class="comment">*</span><span class="comment">           tree: OLDNCL = number of clusters on current level
</span>            OLDNCL = NCLUS
<span class="comment">*</span><span class="comment">           reset NCLUS to count the number of child clusters
</span>            NCLUS = 0
<span class="comment">*</span><span class="comment">
</span>            PARITY = 1 - PARITY
            IF( PARITY.EQ.0 ) THEN
               OLDCLS = IINDC1
               NEWCLS = IINDC2
            ELSE
               OLDCLS = IINDC2
               NEWCLS = IINDC1
            END IF
<span class="comment">*</span><span class="comment">           Process the clusters on the current level
</span>            DO 150 I = 1, OLDNCL
               J = OLDCLS + 2*I
<span class="comment">*</span><span class="comment">              OLDFST, OLDLST = first, last index of current cluster.
</span><span class="comment">*</span><span class="comment">                               cluster indices start with 1 and are relative
</span><span class="comment">*</span><span class="comment">                               to WBEGIN when accessing W, WGAP, WERR, Z
</span>               OLDFST = IWORK( J-1 )
               OLDLST = IWORK( J )
               IF( NDEPTH.GT.0 ) THEN
<span class="comment">*</span><span class="comment">                 Retrieve relatively robust representation (RRR) of cluster
</span><span class="comment">*</span><span class="comment">                 that has been computed at the previous level
</span><span class="comment">*</span><span class="comment">                 The RRR is stored in Z and overwritten once the eigenvectors
</span><span class="comment">*</span><span class="comment">                 have been computed or when the cluster is refined
</span>
                  IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
<span class="comment">*</span><span class="comment">                    Get representation from location of the leftmost evalue
</span><span class="comment">*</span><span class="comment">                    of the cluster
</span>                     J = WBEGIN + OLDFST - 1
                  ELSE
                     IF(WBEGIN+OLDFST-1.LT.DOL) THEN
<span class="comment">*</span><span class="comment">                       Get representation from the left end of Z array
</span>                        J = DOL - 1
                     ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN
<span class="comment">*</span><span class="comment">                       Get representation from the right end of Z array
</span>                        J = DOU
                     ELSE
                        J = WBEGIN + OLDFST - 1
                     ENDIF
                  ENDIF
                  DO 45 K = 1, IN - 1
                     D( IBEGIN+K-1 ) = DBLE( Z( IBEGIN+K-1,
     $                                 J ) )
                     L( IBEGIN+K-1 ) = DBLE( Z( IBEGIN+K-1,
     $                                 J+1 ) )
   45             CONTINUE
                  D( IEND ) = DBLE( Z( IEND, J ) )
                  SIGMA = DBLE( Z( IEND, J+1 ) )

<span class="comment">*</span><span class="comment">                 Set the corresponding entries in Z to zero
</span>                  CALL <a name="ZLASET.425"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'Full'</span>, IN, 2, CZERO, CZERO,
     $                         Z( IBEGIN, J), LDZ )
               END IF

<span class="comment">*</span><span class="comment">              Compute DL and DLL of current RRR
</span>               DO 50 J = IBEGIN, IEND-1
                  TMP = D( J )*L( J )
                  WORK( INDLD-1+J ) = TMP
                  WORK( INDLLD-1+J ) = TMP*L( J )
   50          CONTINUE

               IF( NDEPTH.GT.0 ) THEN
<span class="comment">*</span><span class="comment">                 P and Q are index of the first and last eigenvalue to compute
</span><span class="comment">*</span><span class="comment">                 within the current block
</span>                  P = INDEXW( WBEGIN-1+OLDFST )
                  Q = INDEXW( WBEGIN-1+OLDLST )
<span class="comment">*</span><span class="comment">                 Offset for the arrays WORK, WGAP and WERR, i.e., th P-OFFSET
</span><span class="comment">*</span><span class="comment">                 thru' Q-OFFSET elements of these arrays are to be used.
</span>C                  OFFSET = P-OLDFST
                  OFFSET = INDEXW( WBEGIN ) - 1
<span class="comment">*</span><span class="comment">                 perform limited bisection (if necessary) to get approximate
</span><span class="comment">*</span><span class="comment">                 eigenvalues to the precision needed.
</span>                  CALL <a name="DLARRB.447"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a>( IN, D( IBEGIN ),
     $                         WORK(INDLLD+IBEGIN-1),
     $                         P, Q, RTOL1, RTOL2, OFFSET,
     $                         WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN),
     $                         WORK( INDWRK ), IWORK( IINDWK ),
     $                         PIVMIN, SPDIAM, IN, IINFO )
                  IF( IINFO.NE.0 ) THEN
                     INFO = -1
                     RETURN
                  ENDIF
<span class="comment">*</span><span class="comment">                 We also recompute the extremal gaps. W holds all eigenvalues
</span><span class="comment">*</span><span class="comment">                 of the unshifted matrix and must be used for computation
</span><span class="comment">*</span><span class="comment">                 of WGAP, the entries of WORK might stem from RRRs with
</span><span class="comment">*</span><span class="comment">                 different shifts. The gaps from WBEGIN-1+OLDFST to
</span><span class="comment">*</span><span class="comment">                 WBEGIN-1+OLDLST are correctly computed in <a name="DLARRB.461"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a>.
</span><span class="comment">*</span><span class="comment">                 However, we only allow the gaps to become greater since
</span><span class="comment">*</span><span class="comment">                 this is what should happen when we decrease WERR
</span>                  IF( OLDFST.GT.1) THEN
                     WGAP( WBEGIN+OLDFST-2 ) =
     $             MAX(WGAP(WBEGIN+OLDFST-2),
     $                 W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1)
     $                 - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) )
                  ENDIF
                  IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN
                     WGAP( WBEGIN+OLDLST-1 ) =
     $               MAX(WGAP(WBEGIN+OLDLST-1),
     $                   W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST)
     $                   - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) )
                  ENDIF
<span class="comment">*</span><span class="comment">                 Each time the eigenvalues in WORK get refined, we store
</span><span class="comment">*</span><span class="comment">                 the newly found approximation with all shifts applied in W
</span>                  DO 53 J=OLDFST,OLDLST
                     W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA
 53               CONTINUE
               END IF

<span class="comment">*</span><span class="comment">              Process the current node.
</span>               NEWFST = OLDFST
               DO 140 J = OLDFST, OLDLST
                  IF( J.EQ.OLDLST ) THEN
<span class="comment">*</span><span class="comment">                    we are at the right end of the cluster, this is also the
</span><span class="comment">*</span><span class="comment">                    boundary of the child cluster
</span>                     NEWLST = J
                  ELSE IF ( WGAP( WBEGIN + J -1).GE.
     $                    MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN
<span class="comment">*</span><span class="comment">                    the right relative gap is big enough, the child cluster
</span><span class="comment">*</span><span class="comment">                    (NEWFST,..,NEWLST) is well separated from the following
</span>                     NEWLST = J
                   ELSE
<span class="comment">*</span><span class="comment">                    inside a child cluster, the relative gap is not
</span><span class="comment">*</span><span class="comment">                    big enough.
</span>                     GOTO 140
                  END IF

<span class="comment">*</span><span class="comment">                 Compute size of child cluster found
</span>                  NEWSIZ = NEWLST - NEWFST + 1

<span class="comment">*</span><span class="comment">                 NEWFTT is the place in Z where the new RRR or the computed
</span><span class="comment">*</span><span class="comment">                 eigenvector is to be stored
</span>                  IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
<span class="comment">*</span><span class="comment">                    Store representation at location of the leftmost evalue
</span><span class="comment">*</span><span class="comment">                    of the cluster
</span>                     NEWFTT = WBEGIN + NEWFST - 1
                  ELSE
                     IF(WBEGIN+NEWFST-1.LT.DOL) THEN
<span class="comment">*</span><span class="comment">                       Store representation at the left end of Z array