clatbs.f.html

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

HTML
933
字号
</span>      IF( .NOT.UPPER .AND. .NOT.<a name="LSAME.216"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) ) THEN
         INFO = -1
      ELSE IF( .NOT.NOTRAN .AND. .NOT.<a name="LSAME.218"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'T'</span> ) .AND. .NOT.
     $         <a name="LSAME.219"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'C'</span> ) ) THEN
         INFO = -2
      ELSE IF( .NOT.NOUNIT .AND. .NOT.<a name="LSAME.221"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( DIAG, <span class="string">'U'</span> ) ) THEN
         INFO = -3
      ELSE IF( .NOT.<a name="LSAME.223"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( NORMIN, <span class="string">'Y'</span> ) .AND. .NOT.
     $         <a name="LSAME.224"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( NORMIN, <span class="string">'N'</span> ) ) THEN
         INFO = -4
      ELSE IF( N.LT.0 ) THEN
         INFO = -5
      ELSE IF( KD.LT.0 ) THEN
         INFO = -6
      ELSE IF( LDAB.LT.KD+1 ) THEN
         INFO = -8
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.234"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CLATBS.234"></a><a href="clatbs.f.html#CLATBS.1">CLATBS</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><span class="comment">*</span><span class="comment">     Determine machine dependent parameters to control overflow.
</span><span class="comment">*</span><span class="comment">
</span>      SMLNUM = <a name="SLAMCH.245"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Safe minimum'</span> )
      BIGNUM = ONE / SMLNUM
      CALL <a name="SLABAD.247"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SMLNUM, BIGNUM )
      SMLNUM = SMLNUM / <a name="SLAMCH.248"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Precision'</span> )
      BIGNUM = ONE / SMLNUM
      SCALE = ONE
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.252"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( NORMIN, <span class="string">'N'</span> ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute the 1-norm of each column, not including the diagonal.
</span><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">           A is upper triangular.
</span><span class="comment">*</span><span class="comment">
</span>            DO 10 J = 1, N
               JLEN = MIN( KD, J-1 )
               CNORM( J ) = SCASUM( JLEN, AB( KD+1-JLEN, J ), 1 )
   10       CONTINUE
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           A is lower triangular.
</span><span class="comment">*</span><span class="comment">
</span>            DO 20 J = 1, N
               JLEN = MIN( KD, N-J )
               IF( JLEN.GT.0 ) THEN
                  CNORM( J ) = SCASUM( JLEN, AB( 2, J ), 1 )
               ELSE
                  CNORM( J ) = ZERO
               END IF
   20       CONTINUE
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale the column norms by TSCAL if the maximum element in CNORM is
</span><span class="comment">*</span><span class="comment">     greater than BIGNUM/2.
</span><span class="comment">*</span><span class="comment">
</span>      IMAX = ISAMAX( N, CNORM, 1 )
      TMAX = CNORM( IMAX )
      IF( TMAX.LE.BIGNUM*HALF ) THEN
         TSCAL = ONE
      ELSE
         TSCAL = HALF / ( SMLNUM*TMAX )
         CALL SSCAL( N, TSCAL, CNORM, 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute a bound on the computed solution vector to see if the
</span><span class="comment">*</span><span class="comment">     Level 2 BLAS routine CTBSV can be used.
</span><span class="comment">*</span><span class="comment">
</span>      XMAX = ZERO
      DO 30 J = 1, N
         XMAX = MAX( XMAX, CABS2( X( J ) ) )
   30 CONTINUE
      XBND = XMAX
      IF( NOTRAN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute the growth in A * x = b.
</span><span class="comment">*</span><span class="comment">
</span>         IF( UPPER ) THEN
            JFIRST = N
            JLAST = 1
            JINC = -1
            MAIND = KD + 1
         ELSE
            JFIRST = 1
            JLAST = N
            JINC = 1
            MAIND = 1
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( TSCAL.NE.ONE ) THEN
            GROW = ZERO
            GO TO 60
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( NOUNIT ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           A is non-unit triangular.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute GROW = 1/G(j) and XBND = 1/M(j).
</span><span class="comment">*</span><span class="comment">           Initially, G(0) = max{x(i), i=1,...,n}.
</span><span class="comment">*</span><span class="comment">
</span>            GROW = HALF / MAX( XBND, SMLNUM )
            XBND = GROW
            DO 40 J = JFIRST, JLAST, JINC
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Exit the loop if the growth factor is too small.
</span><span class="comment">*</span><span class="comment">
</span>               IF( GROW.LE.SMLNUM )
     $            GO TO 60
<span class="comment">*</span><span class="comment">
</span>               TJJS = AB( MAIND, J )
               TJJ = CABS1( TJJS )
<span class="comment">*</span><span class="comment">
</span>               IF( TJJ.GE.SMLNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 M(j) = G(j-1) / abs(A(j,j))
</span><span class="comment">*</span><span class="comment">
</span>                  XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 M(j) could overflow, set XBND to 0.
</span><span class="comment">*</span><span class="comment">
</span>                  XBND = ZERO
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
</span><span class="comment">*</span><span class="comment">
</span>                  GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 G(j) could overflow, set GROW to 0.
</span><span class="comment">*</span><span class="comment">
</span>                  GROW = ZERO
               END IF
   40       CONTINUE
            GROW = XBND
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           A is unit triangular.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
</span><span class="comment">*</span><span class="comment">
</span>            GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
            DO 50 J = JFIRST, JLAST, JINC
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Exit the loop if the growth factor is too small.
</span><span class="comment">*</span><span class="comment">
</span>               IF( GROW.LE.SMLNUM )
     $            GO TO 60
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              G(j) = G(j-1)*( 1 + CNORM(j) )
</span><span class="comment">*</span><span class="comment">
</span>               GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
   50       CONTINUE
         END IF
   60    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute the growth in A**T * x = b  or  A**H * x = b.
</span><span class="comment">*</span><span class="comment">
</span>         IF( UPPER ) THEN
            JFIRST = 1
            JLAST = N
            JINC = 1
            MAIND = KD + 1
         ELSE
            JFIRST = N
            JLAST = 1
            JINC = -1
            MAIND = 1
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( TSCAL.NE.ONE ) THEN
            GROW = ZERO
            GO TO 90
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( NOUNIT ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           A is non-unit triangular.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute GROW = 1/G(j) and XBND = 1/M(j).
</span><span class="comment">*</span><span class="comment">           Initially, M(0) = max{x(i), i=1,...,n}.
</span><span class="comment">*</span><span class="comment">
</span>            GROW = HALF / MAX( XBND, SMLNUM )
            XBND = GROW
            DO 70 J = JFIRST, JLAST, JINC
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Exit the loop if the growth factor is too small.
</span><span class="comment">*</span><span class="comment">
</span>               IF( GROW.LE.SMLNUM )
     $            GO TO 90
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
</span><span class="comment">*</span><span class="comment">
</span>               XJ = ONE + CNORM( J )
               GROW = MIN( GROW, XBND / XJ )
<span class="comment">*</span><span class="comment">
</span>               TJJS = AB( MAIND, J )
               TJJ = CABS1( TJJS )
<span class="comment">*</span><span class="comment">
</span>               IF( TJJ.GE.SMLNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
</span><span class="comment">*</span><span class="comment">
</span>                  IF( XJ.GT.TJJ )
     $               XBND = XBND*( TJJ / XJ )
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 M(j) could overflow, set XBND to 0.
</span><span class="comment">*</span><span class="comment">
</span>                  XBND = ZERO
               END IF
   70       CONTINUE
            GROW = MIN( GROW, XBND )
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           A is unit triangular.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
</span><span class="comment">*</span><span class="comment">

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?