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