📄 dlaqtr.f.html
字号:
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>dlaqtr.f</title>
<meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!--
body { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default a { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.string { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.string a { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.comment { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.comment a { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
--></style>
</head>
<body>
<pre>
SUBROUTINE <a name="DLAQTR.1"></a><a href="dlaqtr.f.html#DLAQTR.1">DLAQTR</a>( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK,
$ INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> -- LAPACK auxiliary routine (version 3.1) --
</span><span class="comment">*</span><span class="comment"> Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment"> November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Scalar Arguments ..
</span> LOGICAL LREAL, LTRAN
INTEGER INFO, LDT, N
DOUBLE PRECISION SCALE, W
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Array Arguments ..
</span> DOUBLE PRECISION B( * ), T( LDT, * ), WORK( * ), X( * )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Purpose
</span><span class="comment">*</span><span class="comment"> =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> <a name="DLAQTR.20"></a><a href="dlaqtr.f.html#DLAQTR.1">DLAQTR</a> solves the real quasi-triangular system
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> op(T)*p = scale*c, if LREAL = .TRUE.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> or the complex quasi-triangular systems
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> op(T + iB)*(p+iq) = scale*(c+id), if LREAL = .FALSE.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> in real arithmetic, where T is upper quasi-triangular.
</span><span class="comment">*</span><span class="comment"> If LREAL = .FALSE., then the first diagonal block of T must be
</span><span class="comment">*</span><span class="comment"> 1 by 1, B is the specially structured matrix
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> B = [ b(1) b(2) ... b(n) ]
</span><span class="comment">*</span><span class="comment"> [ w ]
</span><span class="comment">*</span><span class="comment"> [ w ]
</span><span class="comment">*</span><span class="comment"> [ . ]
</span><span class="comment">*</span><span class="comment"> [ w ]
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> op(A) = A or A', A' denotes the conjugate transpose of
</span><span class="comment">*</span><span class="comment"> matrix A.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> On input, X = [ c ]. On output, X = [ p ].
</span><span class="comment">*</span><span class="comment"> [ d ] [ q ]
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> This subroutine is designed for the condition number estimation
</span><span class="comment">*</span><span class="comment"> in routine <a name="DTRSNA.45"></a><a href="dtrsna.f.html#DTRSNA.1">DTRSNA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Arguments
</span><span class="comment">*</span><span class="comment"> =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LTRAN (input) LOGICAL
</span><span class="comment">*</span><span class="comment"> On entry, LTRAN specifies the option of conjugate transpose:
</span><span class="comment">*</span><span class="comment"> = .FALSE., op(T+i*B) = T+i*B,
</span><span class="comment">*</span><span class="comment"> = .TRUE., op(T+i*B) = (T+i*B)'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LREAL (input) LOGICAL
</span><span class="comment">*</span><span class="comment"> On entry, LREAL specifies the input matrix structure:
</span><span class="comment">*</span><span class="comment"> = .FALSE., the input is complex
</span><span class="comment">*</span><span class="comment"> = .TRUE., the input is real
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> N (input) INTEGER
</span><span class="comment">*</span><span class="comment"> On entry, N specifies the order of T+i*B. N >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> T (input) DOUBLE PRECISION array, dimension (LDT,N)
</span><span class="comment">*</span><span class="comment"> On entry, T contains a matrix in Schur canonical form.
</span><span class="comment">*</span><span class="comment"> If LREAL = .FALSE., then the first diagonal block of T mu
</span><span class="comment">*</span><span class="comment"> be 1 by 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDT (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the matrix T. LDT >= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> B (input) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment"> On entry, B contains the elements to form the matrix
</span><span class="comment">*</span><span class="comment"> B as described above.
</span><span class="comment">*</span><span class="comment"> If LREAL = .TRUE., B is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> W (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> On entry, W is the diagonal element of the matrix B.
</span><span class="comment">*</span><span class="comment"> If LREAL = .TRUE., W is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> SCALE (output) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> On exit, SCALE is the scale factor.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> X (input/output) DOUBLE PRECISION array, dimension (2*N)
</span><span class="comment">*</span><span class="comment"> On entry, X contains the right hand side of the system.
</span><span class="comment">*</span><span class="comment"> On exit, X is overwritten by the solution.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> WORK (workspace) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> INFO (output) INTEGER
</span><span class="comment">*</span><span class="comment"> On exit, INFO is set to
</span><span class="comment">*</span><span class="comment"> 0: successful exit.
</span><span class="comment">*</span><span class="comment"> 1: the some diagonal 1 by 1 block has been perturbed by
</span><span class="comment">*</span><span class="comment"> a small number SMIN to keep nonsingularity.
</span><span class="comment">*</span><span class="comment"> 2: the some diagonal 2 by 2 block has been perturbed by
</span><span class="comment">*</span><span class="comment"> a small number in <a name="DLALN2.95"></a><a href="dlaln2.f.html#DLALN2.1">DLALN2</a> to keep nonsingularity.
</span><span class="comment">*</span><span class="comment"> NOTE: In the interests of speed, this routine does not
</span><span class="comment">*</span><span class="comment"> check the inputs for errors.
</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> DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> LOGICAL NOTRAN
INTEGER I, IERR, J, J1, J2, JNEXT, K, N1, N2
DOUBLE PRECISION BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW,
$ SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Arrays ..
</span> DOUBLE PRECISION D( 2, 2 ), V( 2, 2 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> INTEGER IDAMAX
DOUBLE PRECISION DASUM, DDOT, <a name="DLAMCH.116"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANGE.116"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
EXTERNAL IDAMAX, DASUM, DDOT, <a name="DLAMCH.117"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANGE.117"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL DAXPY, <a name="DLADIV.120"></a><a href="dladiv.f.html#DLADIV.1">DLADIV</a>, <a name="DLALN2.120"></a><a href="dlaln2.f.html#DLALN2.1">DLALN2</a>, DSCAL
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC ABS, MAX
<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"> Do not test the input parameters for errors
</span><span class="comment">*</span><span class="comment">
</span> NOTRAN = .NOT.LTRAN
INFO = 0
<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"> Set constants to control overflow
</span><span class="comment">*</span><span class="comment">
</span> EPS = <a name="DLAMCH.139"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
SMLNUM = <a name="DLAMCH.140"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> ) / EPS
BIGNUM = ONE / SMLNUM
<span class="comment">*</span><span class="comment">
</span> XNORM = <a name="DLANGE.143"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>( <span class="string">'M'</span>, N, N, T, LDT, D )
IF( .NOT.LREAL )
$ XNORM = MAX( XNORM, ABS( W ), <a name="DLANGE.145"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>( <span class="string">'M'</span>, N, 1, B, N, D ) )
SMIN = MAX( SMLNUM, EPS*XNORM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute 1-norm of each column of strictly upper triangular
</span><span class="comment">*</span><span class="comment"> part of T to control overflow in triangular solver.
</span><span class="comment">*</span><span class="comment">
</span> WORK( 1 ) = ZERO
DO 10 J = 2, N
WORK( J ) = DASUM( J-1, T( 1, J ), 1 )
10 CONTINUE
<span class="comment">*</span><span class="comment">
</span> IF( .NOT.LREAL ) THEN
DO 20 I = 2, N
WORK( I ) = WORK( I ) + ABS( B( I ) )
20 CONTINUE
END IF
<span class="comment">*</span><span class="comment">
</span> N2 = 2*N
N1 = N
IF( .NOT.LREAL )
$ N1 = N2
K = IDAMAX( N1, X, 1 )
XMAX = ABS( X( K ) )
SCALE = ONE
<span class="comment">*</span><span class="comment">
</span> IF( XMAX.GT.BIGNUM ) THEN
SCALE = BIGNUM / XMAX
CALL DSCAL( N1, SCALE, X, 1 )
XMAX = BIGNUM
END IF
<span class="comment">*</span><span class="comment">
</span> IF( LREAL ) THEN
<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"> Solve T*p = scale*c
</span><span class="comment">*</span><span class="comment">
</span> JNEXT = N
DO 30 J = N, 1, -1
IF( J.GT.JNEXT )
$ GO TO 30
J1 = J
J2 = J
JNEXT = J - 1
IF( J.GT.1 ) THEN
IF( T( J, J-1 ).NE.ZERO ) THEN
J1 = J - 1
JNEXT = J - 2
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> IF( J1.EQ.J2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Meet 1 by 1 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Scale to avoid overflow when computing
</span><span class="comment">*</span><span class="comment"> x(j) = b(j)/T(j,j)
</span><span class="comment">*</span><span class="comment">
</span> XJ = ABS( X( J1 ) )
TJJ = ABS( T( J1, J1 ) )
TMP = T( J1, J1 )
IF( TJJ.LT.SMIN ) THEN
TMP = SMIN
TJJ = SMIN
INFO = 1
END IF
<span class="comment">*</span><span class="comment">
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -