slatrs.f.html

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

HTML
726 行
字号
1234
</span><span class="comment">*</span><span class="comment">                      1&lt;=i&lt;=j
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  and we can safely call STRSV if 1/M(n) and 1/G(n) are both greater
</span><span class="comment">*</span><span class="comment">  than max(underflow, 1/overflow).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, HALF, ONE
      PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            NOTRAN, NOUNIT, UPPER
      INTEGER            I, IMAX, J, JFIRST, JINC, JLAST
      REAL               BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
     $                   TMAX, TSCAL, USCAL, XBND, XJ, XMAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.182"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            ISAMAX
      REAL               SASUM, SDOT, <a name="SLAMCH.184"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="LSAME.185"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, ISAMAX, SASUM, SDOT, <a name="SLAMCH.185"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           SAXPY, SSCAL, STRSV, <a name="XERBLA.188"></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          ABS, 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>      INFO = 0
      UPPER = <a name="LSAME.196"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> )
      NOTRAN = <a name="LSAME.197"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'N'</span> )
      NOUNIT = <a name="LSAME.198"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( DIAG, <span class="string">'N'</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>      IF( .NOT.UPPER .AND. .NOT.<a name="LSAME.202"></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.204"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'T'</span> ) .AND. .NOT.
     $         <a name="LSAME.205"></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.207"></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.209"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( NORMIN, <span class="string">'Y'</span> ) .AND. .NOT.
     $         <a name="LSAME.210"></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( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -7
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.218"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SLATRS.218"></a><a href="slatrs.f.html#SLATRS.1">SLATRS</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.229"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Safe minimum'</span> ) / <a name="SLAMCH.229"></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.233"></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
               CNORM( J ) = SASUM( J-1, A( 1, 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 - 1
               CNORM( J ) = SASUM( N-J, A( J+1, J ), 1 )
   20       CONTINUE
            CNORM( N ) = ZERO
         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.
</span><span class="comment">*</span><span class="comment">
</span>      IMAX = ISAMAX( N, CNORM, 1 )
      TMAX = CNORM( IMAX )
      IF( TMAX.LE.BIGNUM ) THEN
         TSCAL = ONE
      ELSE
         TSCAL = ONE / ( 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 STRSV can be used.
</span><span class="comment">*</span><span class="comment">
</span>      J = ISAMAX( N, X, 1 )
      XMAX = ABS( X( J ) )
      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
         ELSE
            JFIRST = 1
            JLAST = N
            JINC = 1
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( TSCAL.NE.ONE ) THEN
            GROW = ZERO
            GO TO 50
         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 = ONE / MAX( XBND, SMLNUM )
            XBND = GROW
            DO 30 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 50
<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>               TJJ = ABS( A( J, J ) )
               XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
               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
   30       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, ONE / MAX( XBND, SMLNUM ) )
            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 50
<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 ) ) )
   40       CONTINUE
         END IF
   50    CONTINUE
<span class="comment">*</span><span class="comment">
1234

slatrs.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 slatrs.f.html 源码文件，采用 HTML 编程语言编写，共 726 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。

虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。

⌨️ 快捷键说明

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