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 + -
显示快捷键?