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