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