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