1 SUBROUTINE DGTCON( NORM, N, DL, D, DU, DU2, IPIV, ANORM, RCOND,
2 $ WORK, IWORK, 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 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**T)*inv(U**T).
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
2 $ WORK, IWORK, 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 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**T)*inv(U**T).
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