cgbtf2.f.html

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

</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER            ICAMAX
      EXTERNAL           ICAMAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CGERU, CSCAL, CSWAP, <a name="XERBLA.101"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          MAX, MIN
<span class="comment">*</span><span class="comment">     ..
</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">     KV is the number of superdiagonals in the factor U, allowing for
</span><span class="comment">*</span><span class="comment">     fill-in.
</span><span class="comment">*</span><span class="comment">
</span>      KV = KU + KL
<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
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( KL.LT.0 ) THEN
         INFO = -3
      ELSE IF( KU.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDAB.LT.KL+KV+1 ) THEN
         INFO = -6
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.128"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CGBTF2.128"></a><a href="cgbtf2.f.html#CGBTF2.1">CGBTF2</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( M.EQ.0 .OR. N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Gaussian elimination with partial pivoting
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set fill-in elements in columns KU+2 to KV to zero.
</span><span class="comment">*</span><span class="comment">
</span>      DO 20 J = KU + 2, MIN( KV, N )
         DO 10 I = KV - J + 2, KL
            AB( I, J ) = ZERO
   10    CONTINUE
   20 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     JU is the index of the last column affected by the current stage
</span><span class="comment">*</span><span class="comment">     of the factorization.
</span><span class="comment">*</span><span class="comment">
</span>      JU = 1
<span class="comment">*</span><span class="comment">
</span>      DO 40 J = 1, MIN( M, N )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set fill-in elements in column J+KV to zero.
</span><span class="comment">*</span><span class="comment">
</span>         IF( J+KV.LE.N ) THEN
            DO 30 I = 1, KL
               AB( I, J+KV ) = ZERO
   30       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Find pivot and test for singularity. KM is the number of
</span><span class="comment">*</span><span class="comment">        subdiagonal elements in the current column.
</span><span class="comment">*</span><span class="comment">
</span>         KM = MIN( KL, M-J )
         JP = ICAMAX( KM+1, AB( KV+1, J ), 1 )
         IPIV( J ) = JP + J - 1
         IF( AB( KV+JP, J ).NE.ZERO ) THEN
            JU = MAX( JU, MIN( J+KU+JP-1, N ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Apply interchange to columns J to JU.
</span><span class="comment">*</span><span class="comment">
</span>            IF( JP.NE.1 )
     $         CALL CSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1,
     $                     AB( KV+1, J ), LDAB-1 )
            IF( KM.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute multipliers.
</span><span class="comment">*</span><span class="comment">
</span>               CALL CSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Update trailing submatrix within the band.
</span><span class="comment">*</span><span class="comment">
</span>               IF( JU.GT.J )
     $            CALL CGERU( KM, JU-J, -ONE, AB( KV+2, J ), 1,
     $                        AB( KV, J+1 ), LDAB-1, AB( KV+1, J+1 ),
     $                        LDAB-1 )
            END IF
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If pivot is zero, set INFO to the index of the pivot
</span><span class="comment">*</span><span class="comment">           unless a zero pivot has already been found.
</span><span class="comment">*</span><span class="comment">
</span>            IF( INFO.EQ.0 )
     $         INFO = J
         END IF
   40 CONTINUE
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CGBTF2.200"></a><a href="cgbtf2.f.html#CGBTF2.1">CGBTF2</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>