1 SUBROUTINE DPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
2 $ 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 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**T).
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**T).
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
2 $ 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 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**T).
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**T).
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