1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
SUBROUTINE CTPT06( RCOND, RCONDC, UPLO, DIAG, N, AP, RWORK, RAT )
*
* -- LAPACK test routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* February 29, 1992
*
* .. Scalar Arguments ..
CHARACTER DIAG, UPLO
INTEGER N
REAL RAT, RCOND, RCONDC
* ..
* .. Array Arguments ..
REAL RWORK( * )
COMPLEX AP( * )
* ..
*
* Purpose
* =======
*
* CTPT06 computes a test ratio comparing RCOND (the reciprocal
* condition number of the triangular matrix A) and RCONDC, the estimate
* computed by CTPCON. Information about the triangular matrix is used
* if one estimate is zero and the other is non-zero to decide if
* underflow in the estimate is justified.
*
* Arguments
* =========
*
* RCOND (input) REAL
* The estimate of the reciprocal condition number obtained by
* forming the explicit inverse of the matrix A and computing
* RCOND = 1/( norm(A) * norm(inv(A)) ).
*
* RCONDC (input) REAL
* The estimate of the reciprocal condition number computed by
* CTPCON.
*
* UPLO (input) CHARACTER
* Specifies whether the matrix A is upper or lower triangular.
* = 'U': Upper triangular
* = 'L': Lower triangular
*
* DIAG (input) CHARACTER
* Specifies whether or not the matrix A is unit triangular.
* = 'N': Non-unit triangular
* = 'U': Unit triangular
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* AP (input) COMPLEX array, dimension (N*(N+1)/2)
* The upper or lower triangular matrix A, packed columnwise in
* a linear array. The j-th column of A is stored in the array
* AP as follows:
* if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
* if UPLO = 'L',
* AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
*
* RWORK (workspace) REAL array, dimension (N)
*
* RAT (output) REAL
* The test ratio. If both RCOND and RCONDC are nonzero,
* RAT = MAX( RCOND, RCONDC )/MIN( RCOND, RCONDC ) - 1.
* If RAT = 0, the two estimates are exactly the same.
*
* =====================================================================
*
* .. Parameters ..
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
* ..
* .. Local Scalars ..
REAL ANORM, BIGNUM, EPS, RMAX, RMIN
* ..
* .. External Functions ..
REAL CLANTP, SLAMCH
EXTERNAL CLANTP, SLAMCH
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* ..
* .. Executable Statements ..
*
EPS = SLAMCH( 'Epsilon' )
RMAX = MAX( RCOND, RCONDC )
RMIN = MIN( RCOND, RCONDC )
*
* Do the easy cases first.
*
IF( RMIN.LT.ZERO ) THEN
*
* Invalid value for RCOND or RCONDC, return 1/EPS.
*
RAT = ONE / EPS
*
ELSE IF( RMIN.GT.ZERO ) THEN
*
* Both estimates are positive, return RMAX/RMIN - 1.
*
RAT = RMAX / RMIN - ONE
*
ELSE IF( RMAX.EQ.ZERO ) THEN
*
* Both estimates zero.
*
RAT = ZERO
*
ELSE
*
* One estimate is zero, the other is non-zero. If the matrix is
* ill-conditioned, return the nonzero estimate multiplied by
* 1/EPS; if the matrix is badly scaled, return the nonzero
* estimate multiplied by BIGNUM/TMAX, where TMAX is the maximum
* element in absolute value in A.
*
BIGNUM = ONE / SLAMCH( 'Safe minimum' )
ANORM = CLANTP( 'M', UPLO, DIAG, N, AP, RWORK )
*
RAT = RMAX*( MIN( BIGNUM / MAX( ONE, ANORM ), ONE / EPS ) )
END IF
*
RETURN
*
* End of CTPT06
*
END
|
