1 SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, 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, KL, KU, LDAB, N
14 DOUBLE PRECISION ANORM, RCOND
15 * ..
16 * .. Array Arguments ..
17 INTEGER IPIV( * ), IWORK( * )
18 DOUBLE PRECISION AB( LDAB, * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DGBCON estimates the reciprocal of the condition number of a real
25 * general band matrix A, in either the 1-norm or the infinity-norm,
26 * using the LU factorization computed by DGBTRF.
27 *
28 * An estimate is obtained for norm(inv(A)), and the reciprocal of the
29 * condition number is computed as
30 * RCOND = 1 / ( norm(A) * norm(inv(A)) ).
31 *
32 * Arguments
33 * =========
34 *
35 * NORM (input) CHARACTER*1
36 * Specifies whether the 1-norm condition number or the
37 * infinity-norm condition number is required:
38 * = '1' or 'O': 1-norm;
39 * = 'I': Infinity-norm.
40 *
41 * N (input) INTEGER
42 * The order of the matrix A. N >= 0.
43 *
44 * KL (input) INTEGER
45 * The number of subdiagonals within the band of A. KL >= 0.
46 *
47 * KU (input) INTEGER
48 * The number of superdiagonals within the band of A. KU >= 0.
49 *
50 * AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
51 * Details of the LU factorization of the band matrix A, as
52 * computed by DGBTRF. U is stored as an upper triangular band
53 * matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
54 * the multipliers used during the factorization are stored in
55 * rows KL+KU+2 to 2*KL+KU+1.
56 *
57 * LDAB (input) INTEGER
58 * The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
59 *
60 * IPIV (input) INTEGER array, dimension (N)
61 * The pivot indices; for 1 <= i <= N, row i of the matrix was
62 * interchanged with row IPIV(i).
63 *
64 * ANORM (input) DOUBLE PRECISION
65 * If NORM = '1' or 'O', the 1-norm of the original matrix A.
66 * If NORM = 'I', the infinity-norm of the original matrix A.
67 *
68 * RCOND (output) DOUBLE PRECISION
69 * The reciprocal of the condition number of the matrix A,
70 * computed as RCOND = 1/(norm(A) * norm(inv(A))).
71 *
72 * WORK (workspace) DOUBLE PRECISION array, dimension (3*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 LNOTI, ONENRM
88 CHARACTER NORMIN
89 INTEGER IX, J, JP, KASE, KASE1, KD, LM
90 DOUBLE PRECISION AINVNM, SCALE, SMLNUM, T
91 * ..
92 * .. Local Arrays ..
93 INTEGER ISAVE( 3 )
94 * ..
95 * .. External Functions ..
96 LOGICAL LSAME
97 INTEGER IDAMAX
98 DOUBLE PRECISION DDOT, DLAMCH
99 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
100 * ..
101 * .. External Subroutines ..
102 EXTERNAL DAXPY, DLACN2, DLATBS, DRSCL, XERBLA
103 * ..
104 * .. Intrinsic Functions ..
105 INTRINSIC ABS, MIN
106 * ..
107 * .. Executable Statements ..
108 *
109 * Test the input parameters.
110 *
111 INFO = 0
112 ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
113 IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
114 INFO = -1
115 ELSE IF( N.LT.0 ) THEN
116 INFO = -2
117 ELSE IF( KL.LT.0 ) THEN
118 INFO = -3
119 ELSE IF( KU.LT.0 ) THEN
120 INFO = -4
121 ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
122 INFO = -6
123 ELSE IF( ANORM.LT.ZERO ) THEN
124 INFO = -8
125 END IF
126 IF( INFO.NE.0 ) THEN
127 CALL XERBLA( 'DGBCON', -INFO )
128 RETURN
129 END IF
130 *
131 * Quick return if possible
132 *
133 RCOND = ZERO
134 IF( N.EQ.0 ) THEN
135 RCOND = ONE
136 RETURN
137 ELSE IF( ANORM.EQ.ZERO ) THEN
138 RETURN
139 END IF
140 *
141 SMLNUM = DLAMCH( 'Safe minimum' )
142 *
143 * Estimate the norm of inv(A).
144 *
145 AINVNM = ZERO
146 NORMIN = 'N'
147 IF( ONENRM ) THEN
148 KASE1 = 1
149 ELSE
150 KASE1 = 2
151 END IF
152 KD = KL + KU + 1
153 LNOTI = KL.GT.0
154 KASE = 0
155 10 CONTINUE
156 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
157 IF( KASE.NE.0 ) THEN
158 IF( KASE.EQ.KASE1 ) THEN
159 *
160 * Multiply by inv(L).
161 *
162 IF( LNOTI ) THEN
163 DO 20 J = 1, N - 1
164 LM = MIN( KL, N-J )
165 JP = IPIV( J )
166 T = WORK( JP )
167 IF( JP.NE.J ) THEN
168 WORK( JP ) = WORK( J )
169 WORK( J ) = T
170 END IF
171 CALL DAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
172 20 CONTINUE
173 END IF
174 *
175 * Multiply by inv(U).
176 *
177 CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
178 $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
179 $ INFO )
180 ELSE
181 *
182 * Multiply by inv(U**T).
183 *
184 CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
185 $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
186 $ INFO )
187 *
188 * Multiply by inv(L**T).
189 *
190 IF( LNOTI ) THEN
191 DO 30 J = N - 1, 1, -1
192 LM = MIN( KL, N-J )
193 WORK( J ) = WORK( J ) - DDOT( LM, AB( KD+1, J ), 1,
194 $ WORK( J+1 ), 1 )
195 JP = IPIV( J )
196 IF( JP.NE.J ) THEN
197 T = WORK( JP )
198 WORK( JP ) = WORK( J )
199 WORK( J ) = T
200 END IF
201 30 CONTINUE
202 END IF
203 END IF
204 *
205 * Divide X by 1/SCALE if doing so will not cause overflow.
206 *
207 NORMIN = 'Y'
208 IF( SCALE.NE.ONE ) THEN
209 IX = IDAMAX( N, WORK, 1 )
210 IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
211 $ GO TO 40
212 CALL DRSCL( N, SCALE, WORK, 1 )
213 END IF
214 GO TO 10
215 END IF
216 *
217 * Compute the estimate of the reciprocal condition number.
218 *
219 IF( AINVNM.NE.ZERO )
220 $ RCOND = ( ONE / AINVNM ) / ANORM
221 *
222 40 CONTINUE
223 RETURN
224 *
225 * End of DGBCON
226 *
227 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, KL, KU, LDAB, N
14 DOUBLE PRECISION ANORM, RCOND
15 * ..
16 * .. Array Arguments ..
17 INTEGER IPIV( * ), IWORK( * )
18 DOUBLE PRECISION AB( LDAB, * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DGBCON estimates the reciprocal of the condition number of a real
25 * general band matrix A, in either the 1-norm or the infinity-norm,
26 * using the LU factorization computed by DGBTRF.
27 *
28 * An estimate is obtained for norm(inv(A)), and the reciprocal of the
29 * condition number is computed as
30 * RCOND = 1 / ( norm(A) * norm(inv(A)) ).
31 *
32 * Arguments
33 * =========
34 *
35 * NORM (input) CHARACTER*1
36 * Specifies whether the 1-norm condition number or the
37 * infinity-norm condition number is required:
38 * = '1' or 'O': 1-norm;
39 * = 'I': Infinity-norm.
40 *
41 * N (input) INTEGER
42 * The order of the matrix A. N >= 0.
43 *
44 * KL (input) INTEGER
45 * The number of subdiagonals within the band of A. KL >= 0.
46 *
47 * KU (input) INTEGER
48 * The number of superdiagonals within the band of A. KU >= 0.
49 *
50 * AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
51 * Details of the LU factorization of the band matrix A, as
52 * computed by DGBTRF. U is stored as an upper triangular band
53 * matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
54 * the multipliers used during the factorization are stored in
55 * rows KL+KU+2 to 2*KL+KU+1.
56 *
57 * LDAB (input) INTEGER
58 * The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
59 *
60 * IPIV (input) INTEGER array, dimension (N)
61 * The pivot indices; for 1 <= i <= N, row i of the matrix was
62 * interchanged with row IPIV(i).
63 *
64 * ANORM (input) DOUBLE PRECISION
65 * If NORM = '1' or 'O', the 1-norm of the original matrix A.
66 * If NORM = 'I', the infinity-norm of the original matrix A.
67 *
68 * RCOND (output) DOUBLE PRECISION
69 * The reciprocal of the condition number of the matrix A,
70 * computed as RCOND = 1/(norm(A) * norm(inv(A))).
71 *
72 * WORK (workspace) DOUBLE PRECISION array, dimension (3*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 LNOTI, ONENRM
88 CHARACTER NORMIN
89 INTEGER IX, J, JP, KASE, KASE1, KD, LM
90 DOUBLE PRECISION AINVNM, SCALE, SMLNUM, T
91 * ..
92 * .. Local Arrays ..
93 INTEGER ISAVE( 3 )
94 * ..
95 * .. External Functions ..
96 LOGICAL LSAME
97 INTEGER IDAMAX
98 DOUBLE PRECISION DDOT, DLAMCH
99 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
100 * ..
101 * .. External Subroutines ..
102 EXTERNAL DAXPY, DLACN2, DLATBS, DRSCL, XERBLA
103 * ..
104 * .. Intrinsic Functions ..
105 INTRINSIC ABS, MIN
106 * ..
107 * .. Executable Statements ..
108 *
109 * Test the input parameters.
110 *
111 INFO = 0
112 ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
113 IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
114 INFO = -1
115 ELSE IF( N.LT.0 ) THEN
116 INFO = -2
117 ELSE IF( KL.LT.0 ) THEN
118 INFO = -3
119 ELSE IF( KU.LT.0 ) THEN
120 INFO = -4
121 ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
122 INFO = -6
123 ELSE IF( ANORM.LT.ZERO ) THEN
124 INFO = -8
125 END IF
126 IF( INFO.NE.0 ) THEN
127 CALL XERBLA( 'DGBCON', -INFO )
128 RETURN
129 END IF
130 *
131 * Quick return if possible
132 *
133 RCOND = ZERO
134 IF( N.EQ.0 ) THEN
135 RCOND = ONE
136 RETURN
137 ELSE IF( ANORM.EQ.ZERO ) THEN
138 RETURN
139 END IF
140 *
141 SMLNUM = DLAMCH( 'Safe minimum' )
142 *
143 * Estimate the norm of inv(A).
144 *
145 AINVNM = ZERO
146 NORMIN = 'N'
147 IF( ONENRM ) THEN
148 KASE1 = 1
149 ELSE
150 KASE1 = 2
151 END IF
152 KD = KL + KU + 1
153 LNOTI = KL.GT.0
154 KASE = 0
155 10 CONTINUE
156 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
157 IF( KASE.NE.0 ) THEN
158 IF( KASE.EQ.KASE1 ) THEN
159 *
160 * Multiply by inv(L).
161 *
162 IF( LNOTI ) THEN
163 DO 20 J = 1, N - 1
164 LM = MIN( KL, N-J )
165 JP = IPIV( J )
166 T = WORK( JP )
167 IF( JP.NE.J ) THEN
168 WORK( JP ) = WORK( J )
169 WORK( J ) = T
170 END IF
171 CALL DAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
172 20 CONTINUE
173 END IF
174 *
175 * Multiply by inv(U).
176 *
177 CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
178 $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
179 $ INFO )
180 ELSE
181 *
182 * Multiply by inv(U**T).
183 *
184 CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
185 $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
186 $ INFO )
187 *
188 * Multiply by inv(L**T).
189 *
190 IF( LNOTI ) THEN
191 DO 30 J = N - 1, 1, -1
192 LM = MIN( KL, N-J )
193 WORK( J ) = WORK( J ) - DDOT( LM, AB( KD+1, J ), 1,
194 $ WORK( J+1 ), 1 )
195 JP = IPIV( J )
196 IF( JP.NE.J ) THEN
197 T = WORK( JP )
198 WORK( JP ) = WORK( J )
199 WORK( J ) = T
200 END IF
201 30 CONTINUE
202 END IF
203 END IF
204 *
205 * Divide X by 1/SCALE if doing so will not cause overflow.
206 *
207 NORMIN = 'Y'
208 IF( SCALE.NE.ONE ) THEN
209 IX = IDAMAX( N, WORK, 1 )
210 IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
211 $ GO TO 40
212 CALL DRSCL( N, SCALE, WORK, 1 )
213 END IF
214 GO TO 10
215 END IF
216 *
217 * Compute the estimate of the reciprocal condition number.
218 *
219 IF( AINVNM.NE.ZERO )
220 $ RCOND = ( ONE / AINVNM ) / ANORM
221 *
222 40 CONTINUE
223 RETURN
224 *
225 * End of DGBCON
226 *
227 END