Annotation of rpl/lapack/lapack/dgtcon.f, revision 1.6
1.1 bertrand 1: SUBROUTINE DGTCON( NORM, N, DL, D, DU, DU2, IPIV, ANORM, RCOND,
2: $ WORK, IWORK, INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
10: *
11: * .. Scalar Arguments ..
12: CHARACTER NORM
13: INTEGER INFO, N
14: DOUBLE PRECISION ANORM, RCOND
15: * ..
16: * .. Array Arguments ..
17: INTEGER IPIV( * ), IWORK( * )
18: DOUBLE PRECISION D( * ), DL( * ), DU( * ), DU2( * ), WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * DGTCON estimates the reciprocal of the condition number of a real
25: * tridiagonal matrix A using the LU factorization as computed by
26: * DGTTRF.
27: *
28: * An estimate is obtained for norm(inv(A)), and the reciprocal of the
29: * condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
30: *
31: * Arguments
32: * =========
33: *
34: * NORM (input) CHARACTER*1
35: * Specifies whether the 1-norm condition number or the
36: * infinity-norm condition number is required:
37: * = '1' or 'O': 1-norm;
38: * = 'I': Infinity-norm.
39: *
40: * N (input) INTEGER
41: * The order of the matrix A. N >= 0.
42: *
43: * DL (input) DOUBLE PRECISION array, dimension (N-1)
44: * The (n-1) multipliers that define the matrix L from the
45: * LU factorization of A as computed by DGTTRF.
46: *
47: * D (input) DOUBLE PRECISION array, dimension (N)
48: * The n diagonal elements of the upper triangular matrix U from
49: * the LU factorization of A.
50: *
51: * DU (input) DOUBLE PRECISION array, dimension (N-1)
52: * The (n-1) elements of the first superdiagonal of U.
53: *
54: * DU2 (input) DOUBLE PRECISION array, dimension (N-2)
55: * The (n-2) elements of the second superdiagonal of U.
56: *
57: * IPIV (input) INTEGER array, dimension (N)
58: * The pivot indices; for 1 <= i <= n, row i of the matrix was
59: * interchanged with row IPIV(i). IPIV(i) will always be either
60: * i or i+1; IPIV(i) = i indicates a row interchange was not
61: * required.
62: *
63: * ANORM (input) DOUBLE PRECISION
64: * If NORM = '1' or 'O', the 1-norm of the original matrix A.
65: * If NORM = 'I', the infinity-norm of the original matrix A.
66: *
67: * RCOND (output) DOUBLE PRECISION
68: * The reciprocal of the condition number of the matrix A,
69: * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
70: * estimate of the 1-norm of inv(A) computed in this routine.
71: *
72: * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
73: *
74: * IWORK (workspace) INTEGER array, dimension (N)
75: *
76: * INFO (output) INTEGER
77: * = 0: successful exit
78: * < 0: if INFO = -i, the i-th argument had an illegal value
79: *
80: * =====================================================================
81: *
82: * .. Parameters ..
83: DOUBLE PRECISION ONE, ZERO
84: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
85: * ..
86: * .. Local Scalars ..
87: LOGICAL ONENRM
88: INTEGER I, KASE, KASE1
89: DOUBLE PRECISION AINVNM
90: * ..
91: * .. Local Arrays ..
92: INTEGER ISAVE( 3 )
93: * ..
94: * .. External Functions ..
95: LOGICAL LSAME
96: EXTERNAL LSAME
97: * ..
98: * .. External Subroutines ..
99: EXTERNAL DGTTRS, DLACN2, XERBLA
100: * ..
101: * .. Executable Statements ..
102: *
103: * Test the input arguments.
104: *
105: INFO = 0
106: ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
107: IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
108: INFO = -1
109: ELSE IF( N.LT.0 ) THEN
110: INFO = -2
111: ELSE IF( ANORM.LT.ZERO ) THEN
112: INFO = -8
113: END IF
114: IF( INFO.NE.0 ) THEN
115: CALL XERBLA( 'DGTCON', -INFO )
116: RETURN
117: END IF
118: *
119: * Quick return if possible
120: *
121: RCOND = ZERO
122: IF( N.EQ.0 ) THEN
123: RCOND = ONE
124: RETURN
125: ELSE IF( ANORM.EQ.ZERO ) THEN
126: RETURN
127: END IF
128: *
129: * Check that D(1:N) is non-zero.
130: *
131: DO 10 I = 1, N
132: IF( D( I ).EQ.ZERO )
133: $ RETURN
134: 10 CONTINUE
135: *
136: AINVNM = ZERO
137: IF( ONENRM ) THEN
138: KASE1 = 1
139: ELSE
140: KASE1 = 2
141: END IF
142: KASE = 0
143: 20 CONTINUE
144: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
145: IF( KASE.NE.0 ) THEN
146: IF( KASE.EQ.KASE1 ) THEN
147: *
148: * Multiply by inv(U)*inv(L).
149: *
150: CALL DGTTRS( 'No transpose', N, 1, DL, D, DU, DU2, IPIV,
151: $ WORK, N, INFO )
152: ELSE
153: *
154: * Multiply by inv(L')*inv(U').
155: *
156: CALL DGTTRS( 'Transpose', N, 1, DL, D, DU, DU2, IPIV, WORK,
157: $ N, INFO )
158: END IF
159: GO TO 20
160: END IF
161: *
162: * Compute the estimate of the reciprocal condition number.
163: *
164: IF( AINVNM.NE.ZERO )
165: $ RCOND = ( ONE / AINVNM ) / ANORM
166: *
167: RETURN
168: *
169: * End of DGTCON
170: *
171: END
CVSweb interface <joel.bertrand@systella.fr>