clatrs.f.html

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

HTML
904
字号
</span>      INFO = 0
      UPPER = <a name="LSAME.209"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> )
      NOTRAN = <a name="LSAME.210"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'N'</span> )
      NOUNIT = <a name="LSAME.211"></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.215"></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.217"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'T'</span> ) .AND. .NOT.
     $         <a name="LSAME.218"></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.220"></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.222"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( NORMIN, <span class="string">'Y'</span> ) .AND. .NOT.
     $         <a name="LSAME.223"></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.231"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CLATRS.231"></a><a href="clatrs.f.html#CLATRS.1">CLATRS</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.242"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Safe minimum'</span> )
      BIGNUM = ONE / SMLNUM
      CALL <a name="SLABAD.244"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SMLNUM, BIGNUM )
      SMLNUM = SMLNUM / <a name="SLAMCH.245"></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.249"></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 ) = SCASUM( 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 ) = SCASUM( 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/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 CTRSV 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
<span class="comment">*</span><span class="comment">
</span>      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 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 = A( J, 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
         ELSE
            JFIRST = N
            JLAST = 1
            JINC = -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 = A( J, 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

⌨️ 快捷键说明

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