Annotation of rpl/lapack/lapack/zgtcon.f, revision 1.5
1.1 bertrand 1: SUBROUTINE ZGTCON( NORM, N, DL, D, DU, DU2, IPIV, ANORM, RCOND,
2: $ WORK, 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 ZLACN2 in place of ZLACON, 10 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( * )
18: COMPLEX*16 D( * ), DL( * ), DU( * ), DU2( * ), WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * ZGTCON estimates the reciprocal of the condition number of a complex
25: * tridiagonal matrix A using the LU factorization as computed by
26: * ZGTTRF.
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) COMPLEX*16 array, dimension (N-1)
44: * The (n-1) multipliers that define the matrix L from the
45: * LU factorization of A as computed by ZGTTRF.
46: *
47: * D (input) COMPLEX*16 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) COMPLEX*16 array, dimension (N-1)
52: * The (n-1) elements of the first superdiagonal of U.
53: *
54: * DU2 (input) COMPLEX*16 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) COMPLEX*16 array, dimension (2*N)
73: *
74: * INFO (output) INTEGER
75: * = 0: successful exit
76: * < 0: if INFO = -i, the i-th argument had an illegal value
77: *
78: * =====================================================================
79: *
80: * .. Parameters ..
81: DOUBLE PRECISION ONE, ZERO
82: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
83: * ..
84: * .. Local Scalars ..
85: LOGICAL ONENRM
86: INTEGER I, KASE, KASE1
87: DOUBLE PRECISION AINVNM
88: * ..
89: * .. Local Arrays ..
90: INTEGER ISAVE( 3 )
91: * ..
92: * .. External Functions ..
93: LOGICAL LSAME
94: EXTERNAL LSAME
95: * ..
96: * .. External Subroutines ..
97: EXTERNAL XERBLA, ZGTTRS, ZLACN2
98: * ..
99: * .. Intrinsic Functions ..
100: INTRINSIC DCMPLX
101: * ..
102: * .. Executable Statements ..
103: *
104: * Test the input arguments.
105: *
106: INFO = 0
107: ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
108: IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
109: INFO = -1
110: ELSE IF( N.LT.0 ) THEN
111: INFO = -2
112: ELSE IF( ANORM.LT.ZERO ) THEN
113: INFO = -8
114: END IF
115: IF( INFO.NE.0 ) THEN
116: CALL XERBLA( 'ZGTCON', -INFO )
117: RETURN
118: END IF
119: *
120: * Quick return if possible
121: *
122: RCOND = ZERO
123: IF( N.EQ.0 ) THEN
124: RCOND = ONE
125: RETURN
126: ELSE IF( ANORM.EQ.ZERO ) THEN
127: RETURN
128: END IF
129: *
130: * Check that D(1:N) is non-zero.
131: *
132: DO 10 I = 1, N
133: IF( D( I ).EQ.DCMPLX( ZERO ) )
134: $ RETURN
135: 10 CONTINUE
136: *
137: AINVNM = ZERO
138: IF( ONENRM ) THEN
139: KASE1 = 1
140: ELSE
141: KASE1 = 2
142: END IF
143: KASE = 0
144: 20 CONTINUE
145: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
146: IF( KASE.NE.0 ) THEN
147: IF( KASE.EQ.KASE1 ) THEN
148: *
149: * Multiply by inv(U)*inv(L).
150: *
151: CALL ZGTTRS( 'No transpose', N, 1, DL, D, DU, DU2, IPIV,
152: $ WORK, N, INFO )
153: ELSE
154: *
155: * Multiply by inv(L')*inv(U').
156: *
157: CALL ZGTTRS( 'Conjugate transpose', N, 1, DL, D, DU, DU2,
158: $ IPIV, WORK, N, INFO )
159: END IF
160: GO TO 20
161: END IF
162: *
163: * Compute the estimate of the reciprocal condition number.
164: *
165: IF( AINVNM.NE.ZERO )
166: $ RCOND = ( ONE / AINVNM ) / ANORM
167: *
168: RETURN
169: *
170: * End of ZGTCON
171: *
172: END
CVSweb interface <joel.bertrand@systella.fr>