ctrsna.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 381 行 · 第 1/2 页
HTML
381 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>ctrsna.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="CTRSNA.1"></a><a href="ctrsna.f.html#CTRSNA.1">CTRSNA</a>( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
$ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
$ INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> -- LAPACK 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"> Modified to call <a name="CLACN2.9"></a><a href="clacn2.f.html#CLACN2.1">CLACN2</a> in place of <a name="CLACON.9"></a><a href="clacon.f.html#CLACON.1">CLACON</a>, 10 Feb 03, SJH.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Scalar Arguments ..
</span> CHARACTER HOWMNY, JOB
INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Array Arguments ..
</span> LOGICAL SELECT( * )
REAL RWORK( * ), S( * ), SEP( * )
COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
$ WORK( LDWORK, * )
<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="CTRSNA.25"></a><a href="ctrsna.f.html#CTRSNA.1">CTRSNA</a> estimates reciprocal condition numbers for specified
</span><span class="comment">*</span><span class="comment"> eigenvalues and/or right eigenvectors of a complex upper triangular
</span><span class="comment">*</span><span class="comment"> matrix T (or of any matrix Q*T*Q**H with Q unitary).
</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"> JOB (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment"> Specifies whether condition numbers are required for
</span><span class="comment">*</span><span class="comment"> eigenvalues (S) or eigenvectors (SEP):
</span><span class="comment">*</span><span class="comment"> = 'E': for eigenvalues only (S);
</span><span class="comment">*</span><span class="comment"> = 'V': for eigenvectors only (SEP);
</span><span class="comment">*</span><span class="comment"> = 'B': for both eigenvalues and eigenvectors (S and SEP).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> HOWMNY (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment"> = 'A': compute condition numbers for all eigenpairs;
</span><span class="comment">*</span><span class="comment"> = 'S': compute condition numbers for selected eigenpairs
</span><span class="comment">*</span><span class="comment"> specified by the array SELECT.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> SELECT (input) LOGICAL array, dimension (N)
</span><span class="comment">*</span><span class="comment"> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
</span><span class="comment">*</span><span class="comment"> condition numbers are required. To select condition numbers
</span><span class="comment">*</span><span class="comment"> for the j-th eigenpair, SELECT(j) must be set to .TRUE..
</span><span class="comment">*</span><span class="comment"> If HOWMNY = 'A', SELECT is not referenced.
</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"> The order of the matrix T. N >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> T (input) COMPLEX array, dimension (LDT,N)
</span><span class="comment">*</span><span class="comment"> The upper triangular matrix T.
</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 array T. LDT >= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> VL (input) COMPLEX array, dimension (LDVL,M)
</span><span class="comment">*</span><span class="comment"> If JOB = 'E' or 'B', VL must contain left eigenvectors of T
</span><span class="comment">*</span><span class="comment"> (or of any Q*T*Q**H with Q unitary), corresponding to the
</span><span class="comment">*</span><span class="comment"> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
</span><span class="comment">*</span><span class="comment"> must be stored in consecutive columns of VL, as returned by
</span><span class="comment">*</span><span class="comment"> <a name="CHSEIN.64"></a><a href="chsein.f.html#CHSEIN.1">CHSEIN</a> or <a name="CTREVC.64"></a><a href="ctrevc.f.html#CTREVC.1">CTREVC</a>.
</span><span class="comment">*</span><span class="comment"> If JOB = 'V', VL is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDVL (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array VL.
</span><span class="comment">*</span><span class="comment"> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> VR (input) COMPLEX array, dimension (LDVR,M)
</span><span class="comment">*</span><span class="comment"> If JOB = 'E' or 'B', VR must contain right eigenvectors of T
</span><span class="comment">*</span><span class="comment"> (or of any Q*T*Q**H with Q unitary), corresponding to the
</span><span class="comment">*</span><span class="comment"> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
</span><span class="comment">*</span><span class="comment"> must be stored in consecutive columns of VR, as returned by
</span><span class="comment">*</span><span class="comment"> <a name="CHSEIN.76"></a><a href="chsein.f.html#CHSEIN.1">CHSEIN</a> or <a name="CTREVC.76"></a><a href="ctrevc.f.html#CTREVC.1">CTREVC</a>.
</span><span class="comment">*</span><span class="comment"> If JOB = 'V', VR is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDVR (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array VR.
</span><span class="comment">*</span><span class="comment"> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> S (output) REAL array, dimension (MM)
</span><span class="comment">*</span><span class="comment"> If JOB = 'E' or 'B', the reciprocal condition numbers of the
</span><span class="comment">*</span><span class="comment"> selected eigenvalues, stored in consecutive elements of the
</span><span class="comment">*</span><span class="comment"> array. Thus S(j), SEP(j), and the j-th columns of VL and VR
</span><span class="comment">*</span><span class="comment"> all correspond to the same eigenpair (but not in general the
</span><span class="comment">*</span><span class="comment"> j-th eigenpair, unless all eigenpairs are selected).
</span><span class="comment">*</span><span class="comment"> If JOB = 'V', S is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> SEP (output) REAL array, dimension (MM)
</span><span class="comment">*</span><span class="comment"> If JOB = 'V' or 'B', the estimated reciprocal condition
</span><span class="comment">*</span><span class="comment"> numbers of the selected eigenvectors, stored in consecutive
</span><span class="comment">*</span><span class="comment"> elements of the array.
</span><span class="comment">*</span><span class="comment"> If JOB = 'E', SEP is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> MM (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of elements in the arrays S (if JOB = 'E' or 'B')
</span><span class="comment">*</span><span class="comment"> and/or SEP (if JOB = 'V' or 'B'). MM >= M.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> M (output) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of elements of the arrays S and/or SEP actually
</span><span class="comment">*</span><span class="comment"> used to store the estimated condition numbers.
</span><span class="comment">*</span><span class="comment"> If HOWMNY = 'A', M is set to N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> WORK (workspace) COMPLEX array, dimension (LDWORK,N+6)
</span><span class="comment">*</span><span class="comment"> If JOB = 'E', WORK is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDWORK (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array WORK.
</span><span class="comment">*</span><span class="comment"> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> RWORK (workspace) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment"> If JOB = 'E', RWORK is not referenced.
</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"> = 0: successful exit
</span><span class="comment">*</span><span class="comment"> < 0: if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Further Details
</span><span class="comment">*</span><span class="comment"> ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The reciprocal of the condition number of an eigenvalue lambda is
</span><span class="comment">*</span><span class="comment"> defined as
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> S(lambda) = |v'*u| / (norm(u)*norm(v))
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> where u and v are the right and left eigenvectors of T corresponding
</span><span class="comment">*</span><span class="comment"> to lambda; v' denotes the conjugate transpose of v, and norm(u)
</span><span class="comment">*</span><span class="comment"> denotes the Euclidean norm. These reciprocal condition numbers always
</span><span class="comment">*</span><span class="comment"> lie between zero (very badly conditioned) and one (very well
</span><span class="comment">*</span><span class="comment"> conditioned). If n = 1, S(lambda) is defined to be 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> An approximate error bound for a computed eigenvalue W(i) is given by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> EPS * norm(T) / S(i)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> where EPS is the machine precision.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The reciprocal of the condition number of the right eigenvector u
</span><span class="comment">*</span><span class="comment"> corresponding to lambda is defined as follows. Suppose
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> T = ( lambda c )
</span><span class="comment">*</span><span class="comment"> ( 0 T22 )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Then the reciprocal condition number is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> where sigma-min denotes the smallest singular value. We approximate
</span><span class="comment">*</span><span class="comment"> the smallest singular value by the reciprocal of an estimate of the
</span><span class="comment">*</span><span class="comment"> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
</span><span class="comment">*</span><span class="comment"> defined to be abs(T(1,1)).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> An approximate error bound for a computed right eigenvector VR(i)
</span><span class="comment">*</span><span class="comment"> is given by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> EPS * norm(T) / SEP(i)
</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, ONE
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0+0 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> LOGICAL SOMCON, WANTBH, WANTS, WANTSP
CHARACTER NORMIN
INTEGER I, IERR, IX, J, K, KASE, KS
REAL BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
$ XNORM
COMPLEX CDUM, PROD
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?