Annotation of rpl/lapack/lapack/dpbcon.f, revision 1.2
1.1 bertrand 1: SUBROUTINE DPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
2: $ 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 UPLO
13: INTEGER INFO, KD, LDAB, N
14: DOUBLE PRECISION ANORM, RCOND
15: * ..
16: * .. Array Arguments ..
17: INTEGER IWORK( * )
18: DOUBLE PRECISION AB( LDAB, * ), WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * DPBCON estimates the reciprocal of the condition number (in the
25: * 1-norm) of a real symmetric positive definite band matrix using the
26: * Cholesky factorization A = U**T*U or A = L*L**T computed by DPBTRF.
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: * UPLO (input) CHARACTER*1
35: * = 'U': Upper triangular factor stored in AB;
36: * = 'L': Lower triangular factor stored in AB.
37: *
38: * N (input) INTEGER
39: * The order of the matrix A. N >= 0.
40: *
41: * KD (input) INTEGER
42: * The number of superdiagonals of the matrix A if UPLO = 'U',
43: * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
44: *
45: * AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
46: * The triangular factor U or L from the Cholesky factorization
47: * A = U**T*U or A = L*L**T of the band matrix A, stored in the
48: * first KD+1 rows of the array. The j-th column of U or L is
49: * stored in the j-th column of the array AB as follows:
50: * if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;
51: * if UPLO ='L', AB(1+i-j,j) = L(i,j) for j<=i<=min(n,j+kd).
52: *
53: * LDAB (input) INTEGER
54: * The leading dimension of the array AB. LDAB >= KD+1.
55: *
56: * ANORM (input) DOUBLE PRECISION
57: * The 1-norm (or infinity-norm) of the symmetric band matrix A.
58: *
59: * RCOND (output) DOUBLE PRECISION
60: * The reciprocal of the condition number of the matrix A,
61: * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
62: * estimate of the 1-norm of inv(A) computed in this routine.
63: *
64: * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
65: *
66: * IWORK (workspace) INTEGER array, dimension (N)
67: *
68: * INFO (output) INTEGER
69: * = 0: successful exit
70: * < 0: if INFO = -i, the i-th argument had an illegal value
71: *
72: * =====================================================================
73: *
74: * .. Parameters ..
75: DOUBLE PRECISION ONE, ZERO
76: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
77: * ..
78: * .. Local Scalars ..
79: LOGICAL UPPER
80: CHARACTER NORMIN
81: INTEGER IX, KASE
82: DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
83: * ..
84: * .. Local Arrays ..
85: INTEGER ISAVE( 3 )
86: * ..
87: * .. External Functions ..
88: LOGICAL LSAME
89: INTEGER IDAMAX
90: DOUBLE PRECISION DLAMCH
91: EXTERNAL LSAME, IDAMAX, DLAMCH
92: * ..
93: * .. External Subroutines ..
94: EXTERNAL DLACN2, DLATBS, DRSCL, XERBLA
95: * ..
96: * .. Intrinsic Functions ..
97: INTRINSIC ABS
98: * ..
99: * .. Executable Statements ..
100: *
101: * Test the input parameters.
102: *
103: INFO = 0
104: UPPER = LSAME( UPLO, 'U' )
105: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
106: INFO = -1
107: ELSE IF( N.LT.0 ) THEN
108: INFO = -2
109: ELSE IF( KD.LT.0 ) THEN
110: INFO = -3
111: ELSE IF( LDAB.LT.KD+1 ) THEN
112: INFO = -5
113: ELSE IF( ANORM.LT.ZERO ) THEN
114: INFO = -6
115: END IF
116: IF( INFO.NE.0 ) THEN
117: CALL XERBLA( 'DPBCON', -INFO )
118: RETURN
119: END IF
120: *
121: * Quick return if possible
122: *
123: RCOND = ZERO
124: IF( N.EQ.0 ) THEN
125: RCOND = ONE
126: RETURN
127: ELSE IF( ANORM.EQ.ZERO ) THEN
128: RETURN
129: END IF
130: *
131: SMLNUM = DLAMCH( 'Safe minimum' )
132: *
133: * Estimate the 1-norm of the inverse.
134: *
135: KASE = 0
136: NORMIN = 'N'
137: 10 CONTINUE
138: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
139: IF( KASE.NE.0 ) THEN
140: IF( UPPER ) THEN
141: *
142: * Multiply by inv(U').
143: *
144: CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
145: $ KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ),
146: $ INFO )
147: NORMIN = 'Y'
148: *
149: * Multiply by inv(U).
150: *
151: CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
152: $ KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ),
153: $ INFO )
154: ELSE
155: *
156: * Multiply by inv(L).
157: *
158: CALL DLATBS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
159: $ KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ),
160: $ INFO )
161: NORMIN = 'Y'
162: *
163: * Multiply by inv(L').
164: *
165: CALL DLATBS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N,
166: $ KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ),
167: $ INFO )
168: END IF
169: *
170: * Multiply by 1/SCALE if doing so will not cause overflow.
171: *
172: SCALE = SCALEL*SCALEU
173: IF( SCALE.NE.ONE ) THEN
174: IX = IDAMAX( N, WORK, 1 )
175: IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
176: $ GO TO 20
177: CALL DRSCL( N, SCALE, WORK, 1 )
178: END IF
179: GO TO 10
180: END IF
181: *
182: * Compute the estimate of the reciprocal condition number.
183: *
184: IF( AINVNM.NE.ZERO )
185: $ RCOND = ( ONE / AINVNM ) / ANORM
186: *
187: 20 CONTINUE
188: *
189: RETURN
190: *
191: * End of DPBCON
192: *
193: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dpbcon.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack