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