1 SUBROUTINE ZGTCON( NORM, N, DL, D, DU, DU2, IPIV, ANORM, RCOND,
2 $ WORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
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**H)*inv(U**H).
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
2 $ WORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
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**H)*inv(U**H).
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