zpbstf.f.html

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

      UPPER = <a name="LSAME.126"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> )
      IF( .NOT.UPPER .AND. .NOT.<a name="LSAME.127"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( KD.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDAB.LT.KD+1 ) THEN
         INFO = -5
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.137"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZPBSTF.137"></a><a href="zpbstf.f.html#ZPBSTF.1">ZPBSTF</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      KLD = MAX( 1, LDAB-1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set the splitting point m.
</span><span class="comment">*</span><span class="comment">
</span>      M = ( N+KD ) / 2
<span class="comment">*</span><span class="comment">
</span>      IF( UPPER ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Factorize A(m+1:n,m+1:n) as L**H*L, and update A(1:m,1:m).
</span><span class="comment">*</span><span class="comment">
</span>         DO 10 J = N, M + 1, -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute s(j,j) and test for non-positive-definiteness.
</span><span class="comment">*</span><span class="comment">
</span>            AJJ = DBLE( AB( KD+1, J ) )
            IF( AJJ.LE.ZERO ) THEN
               AB( KD+1, J ) = AJJ
               GO TO 50
            END IF
            AJJ = SQRT( AJJ )
            AB( KD+1, J ) = AJJ
            KM = MIN( J-1, KD )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute elements j-km:j-1 of the j-th column and update the
</span><span class="comment">*</span><span class="comment">           the leading submatrix within the band.
</span><span class="comment">*</span><span class="comment">
</span>            CALL ZDSCAL( KM, ONE / AJJ, AB( KD+1-KM, J ), 1 )
            CALL ZHER( <span class="string">'Upper'</span>, KM, -ONE, AB( KD+1-KM, J ), 1,
     $                 AB( KD+1, J-KM ), KLD )
   10    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Factorize the updated submatrix A(1:m,1:m) as U**H*U.
</span><span class="comment">*</span><span class="comment">
</span>         DO 20 J = 1, M
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute s(j,j) and test for non-positive-definiteness.
</span><span class="comment">*</span><span class="comment">
</span>            AJJ = DBLE( AB( KD+1, J ) )
            IF( AJJ.LE.ZERO ) THEN
               AB( KD+1, J ) = AJJ
               GO TO 50
            END IF
            AJJ = SQRT( AJJ )
            AB( KD+1, J ) = AJJ
            KM = MIN( KD, M-J )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute elements j+1:j+km of the j-th row and update the
</span><span class="comment">*</span><span class="comment">           trailing submatrix within the band.
</span><span class="comment">*</span><span class="comment">
</span>            IF( KM.GT.0 ) THEN
               CALL ZDSCAL( KM, ONE / AJJ, AB( KD, J+1 ), KLD )
               CALL <a name="ZLACGV.197"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( KM, AB( KD, J+1 ), KLD )
               CALL ZHER( <span class="string">'Upper'</span>, KM, -ONE, AB( KD, J+1 ), KLD,
     $                    AB( KD+1, J+1 ), KLD )
               CALL <a name="ZLACGV.200"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( KM, AB( KD, J+1 ), KLD )
            END IF
   20    CONTINUE
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Factorize A(m+1:n,m+1:n) as L**H*L, and update A(1:m,1:m).
</span><span class="comment">*</span><span class="comment">
</span>         DO 30 J = N, M + 1, -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute s(j,j) and test for non-positive-definiteness.
</span><span class="comment">*</span><span class="comment">
</span>            AJJ = DBLE( AB( 1, J ) )
            IF( AJJ.LE.ZERO ) THEN
               AB( 1, J ) = AJJ
               GO TO 50
            END IF
            AJJ = SQRT( AJJ )
            AB( 1, J ) = AJJ
            KM = MIN( J-1, KD )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute elements j-km:j-1 of the j-th row and update the
</span><span class="comment">*</span><span class="comment">           trailing submatrix within the band.
</span><span class="comment">*</span><span class="comment">
</span>            CALL ZDSCAL( KM, ONE / AJJ, AB( KM+1, J-KM ), KLD )
            CALL <a name="ZLACGV.224"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( KM, AB( KM+1, J-KM ), KLD )
            CALL ZHER( <span class="string">'Lower'</span>, KM, -ONE, AB( KM+1, J-KM ), KLD,
     $                 AB( 1, J-KM ), KLD )
            CALL <a name="ZLACGV.227"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( KM, AB( KM+1, J-KM ), KLD )
   30    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Factorize the updated submatrix A(1:m,1:m) as U**H*U.
</span><span class="comment">*</span><span class="comment">
</span>         DO 40 J = 1, M
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute s(j,j) and test for non-positive-definiteness.
</span><span class="comment">*</span><span class="comment">
</span>            AJJ = DBLE( AB( 1, J ) )
            IF( AJJ.LE.ZERO ) THEN
               AB( 1, J ) = AJJ
               GO TO 50
            END IF
            AJJ = SQRT( AJJ )
            AB( 1, J ) = AJJ
            KM = MIN( KD, M-J )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute elements j+1:j+km of the j-th column and update the
</span><span class="comment">*</span><span class="comment">           trailing submatrix within the band.
</span><span class="comment">*</span><span class="comment">
</span>            IF( KM.GT.0 ) THEN
               CALL ZDSCAL( KM, ONE / AJJ, AB( 2, J ), 1 )
               CALL ZHER( <span class="string">'Lower'</span>, KM, -ONE, AB( 2, J ), 1,
     $                    AB( 1, J+1 ), KLD )
            END IF
   40    CONTINUE
      END IF
      RETURN
<span class="comment">*</span><span class="comment">
</span>   50 CONTINUE
      INFO = J
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZPBSTF.261"></a><a href="zpbstf.f.html#ZPBSTF.1">ZPBSTF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>