chbgvx.f.html

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

     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Form a split Cholesky factorization of B.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CPBSTF.264"></a><a href="cpbstf.f.html#CPBSTF.1">CPBSTF</a>( UPLO, N, KB, BB, LDBB, INFO )
      IF( INFO.NE.0 ) THEN
         INFO = N + INFO
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Transform problem to standard eigenvalue problem.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CHBGST.272"></a><a href="chbgst.f.html#CHBGST.1">CHBGST</a>( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Q, LDQ,
     $             WORK, RWORK, IINFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve the standard eigenvalue problem.
</span><span class="comment">*</span><span class="comment">     Reduce Hermitian band matrix to tridiagonal form.
</span><span class="comment">*</span><span class="comment">
</span>      INDD = 1
      INDE = INDD + N
      INDRWK = INDE + N
      INDWRK = 1
      IF( WANTZ ) THEN
         VECT = <span class="string">'U'</span>
      ELSE
         VECT = <span class="string">'N'</span>
      END IF
      CALL <a name="CHBTRD.287"></a><a href="chbtrd.f.html#CHBTRD.1">CHBTRD</a>( VECT, UPLO, N, KA, AB, LDAB, RWORK( INDD ),
     $             RWORK( INDE ), Q, LDQ, WORK( INDWRK ), IINFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If all eigenvalues are desired and ABSTOL is less than or equal
</span><span class="comment">*</span><span class="comment">     to zero, then call <a name="SSTERF.291"></a><a href="ssterf.f.html#SSTERF.1">SSTERF</a> or <a name="CSTEQR.291"></a><a href="csteqr.f.html#CSTEQR.1">CSTEQR</a>.  If this fails for some
</span><span class="comment">*</span><span class="comment">     eigenvalue, then try <a name="SSTEBZ.292"></a><a href="sstebz.f.html#SSTEBZ.1">SSTEBZ</a>.
</span><span class="comment">*</span><span class="comment">
</span>      TEST = .FALSE.
      IF( INDEIG ) THEN
         IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
            TEST = .TRUE.
         END IF
      END IF
      IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN
         CALL SCOPY( N, RWORK( INDD ), 1, W, 1 )
         INDEE = INDRWK + 2*N
         CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 )
         IF( .NOT.WANTZ ) THEN
            CALL <a name="SSTERF.305"></a><a href="ssterf.f.html#SSTERF.1">SSTERF</a>( N, W, RWORK( INDEE ), INFO )
         ELSE
            CALL <a name="CLACPY.307"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'A'</span>, N, N, Q, LDQ, Z, LDZ )
            CALL <a name="CSTEQR.308"></a><a href="csteqr.f.html#CSTEQR.1">CSTEQR</a>( JOBZ, N, W, RWORK( INDEE ), Z, LDZ,
     $                   RWORK( INDRWK ), INFO )
            IF( INFO.EQ.0 ) THEN
               DO 10 I = 1, N
                  IFAIL( I ) = 0
   10          CONTINUE
            END IF
         END IF
         IF( INFO.EQ.0 ) THEN
            M = N
            GO TO 30
         END IF
         INFO = 0
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Otherwise, call <a name="SSTEBZ.323"></a><a href="sstebz.f.html#SSTEBZ.1">SSTEBZ</a> and, if eigenvectors are desired,
</span><span class="comment">*</span><span class="comment">     call <a name="CSTEIN.324"></a><a href="cstein.f.html#CSTEIN.1">CSTEIN</a>.
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         ORDER = <span class="string">'B'</span>
      ELSE
         ORDER = <span class="string">'E'</span>
      END IF
      INDIBL = 1
      INDISP = INDIBL + N
      INDIWK = INDISP + N
      CALL <a name="SSTEBZ.334"></a><a href="sstebz.f.html#SSTEBZ.1">SSTEBZ</a>( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL,
     $             RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W,
     $             IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ),
     $             IWORK( INDIWK ), INFO )
<span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         CALL <a name="CSTEIN.340"></a><a href="cstein.f.html#CSTEIN.1">CSTEIN</a>( N, RWORK( INDD ), RWORK( INDE ), M, W,
     $                IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
     $                RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Apply unitary matrix used in reduction to tridiagonal
</span><span class="comment">*</span><span class="comment">        form to eigenvectors returned by <a name="CSTEIN.345"></a><a href="cstein.f.html#CSTEIN.1">CSTEIN</a>.
</span><span class="comment">*</span><span class="comment">
</span>         DO 20 J = 1, M
            CALL CCOPY( N, Z( 1, J ), 1, WORK( 1 ), 1 )
            CALL CGEMV( <span class="string">'N'</span>, N, N, CONE, Q, LDQ, WORK, 1, CZERO,
     $                  Z( 1, J ), 1 )
   20    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If eigenvalues are not in order, then sort them, along with
</span><span class="comment">*</span><span class="comment">     eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         DO 50 J = 1, M - 1
            I = 0
            TMP1 = W( J )
            DO 40 JJ = J + 1, M
               IF( W( JJ ).LT.TMP1 ) THEN
                  I = JJ
                  TMP1 = W( JJ )
               END IF
   40       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( I.NE.0 ) THEN
               ITMP1 = IWORK( INDIBL+I-1 )
               W( I ) = W( J )
               IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
               W( J ) = TMP1
               IWORK( INDIBL+J-1 ) = ITMP1
               CALL CSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
               IF( INFO.NE.0 ) THEN
                  ITMP1 = IFAIL( I )
                  IFAIL( I ) = IFAIL( J )
                  IFAIL( J ) = ITMP1
               END IF
            END IF
   50    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CHBGVX.388"></a><a href="chbgvx.f.html#CHBGVX.1">CHBGVX</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>