spbstf.f.html

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

</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      UPPER = <a name="LSAME.125"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> )
      IF( .NOT.UPPER .AND. .NOT.<a name="LSAME.126"></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.136"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SPBSTF.136"></a><a href="spbstf.f.html#SPBSTF.1">SPBSTF</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**T*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 = AB( KD+1, J )
            IF( AJJ.LE.ZERO )
     $         GO TO 50
            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 SSCAL( KM, ONE / AJJ, AB( KD+1-KM, J ), 1 )
            CALL SSYR( <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**T*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 = AB( KD+1, J )
            IF( AJJ.LE.ZERO )
     $         GO TO 50
            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 SSCAL( KM, ONE / AJJ, AB( KD, J+1 ), KLD )
               CALL SSYR( <span class="string">'Upper'</span>, KM, -ONE, AB( KD, J+1 ), KLD,
     $                    AB( KD+1, 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**T*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 = AB( 1, J )
            IF( AJJ.LE.ZERO )
     $         GO TO 50
            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 SSCAL( KM, ONE / AJJ, AB( KM+1, J-KM ), KLD )
            CALL SSYR( <span class="string">'Lower'</span>, KM, -ONE, AB( KM+1, J-KM ), KLD,
     $                 AB( 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**T*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 = AB( 1, J )
            IF( AJJ.LE.ZERO )
     $         GO TO 50
            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 SSCAL( KM, ONE / AJJ, AB( 2, J ), 1 )
               CALL SSYR( <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="SPBSTF.248"></a><a href="spbstf.f.html#SPBSTF.1">SPBSTF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>