![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgbcon.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to zgbcon.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
2: $ WORK, RWORK, 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, 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').
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').
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