1 SUBROUTINE ZPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
2 $ PIV, RWORK, RESID, RANK )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Craig Lucas, University of Manchester / NAG Ltd.
6 * October, 2008
7 *
8 * .. Scalar Arguments ..
9 DOUBLE PRECISION RESID
10 INTEGER LDA, LDAFAC, LDPERM, N, RANK
11 CHARACTER UPLO
12 * ..
13 * .. Array Arguments ..
14 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ),
15 $ PERM( LDPERM, * )
16 DOUBLE PRECISION RWORK( * )
17 INTEGER PIV( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * ZPST01 reconstructs an Hermitian positive semidefinite matrix A
24 * from its L or U factors and the permutation matrix P and computes
25 * the residual
26 * norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
27 * norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
28 * where EPS is the machine epsilon, L' is the conjugate transpose of L,
29 * and U' is the conjugate transpose of U.
30 *
31 * Arguments
32 * ==========
33 *
34 * UPLO (input) CHARACTER*1
35 * Specifies whether the upper or lower triangular part of the
36 * Hermitian matrix A is stored:
37 * = 'U': Upper triangular
38 * = 'L': Lower triangular
39 *
40 * N (input) INTEGER
41 * The number of rows and columns of the matrix A. N >= 0.
42 *
43 * A (input) COMPLEX*16 array, dimension (LDA,N)
44 * The original Hermitian matrix A.
45 *
46 * LDA (input) INTEGER
47 * The leading dimension of the array A. LDA >= max(1,N)
48 *
49 * AFAC (input) COMPLEX*16 array, dimension (LDAFAC,N)
50 * The factor L or U from the L*L' or U'*U
51 * factorization of A.
52 *
53 * LDAFAC (input) INTEGER
54 * The leading dimension of the array AFAC. LDAFAC >= max(1,N).
55 *
56 * PERM (output) COMPLEX*16 array, dimension (LDPERM,N)
57 * Overwritten with the reconstructed matrix, and then with the
58 * difference P*L*L'*P' - A (or P*U'*U*P' - A)
59 *
60 * LDPERM (input) INTEGER
61 * The leading dimension of the array PERM.
62 * LDAPERM >= max(1,N).
63 *
64 * PIV (input) INTEGER array, dimension (N)
65 * PIV is such that the nonzero entries are
66 * P( PIV( K ), K ) = 1.
67 *
68 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
69 *
70 * RESID (output) DOUBLE PRECISION
71 * If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
72 * If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
73 *
74 * =====================================================================
75 *
76 * .. Parameters ..
77 DOUBLE PRECISION ZERO, ONE
78 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
79 COMPLEX*16 CZERO
80 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
81 * ..
82 * .. Local Scalars ..
83 COMPLEX*16 TC
84 DOUBLE PRECISION ANORM, EPS, TR
85 INTEGER I, J, K
86 * ..
87 * .. External Functions ..
88 COMPLEX*16 ZDOTC
89 DOUBLE PRECISION DLAMCH, ZLANHE
90 LOGICAL LSAME
91 EXTERNAL ZDOTC, DLAMCH, ZLANHE, LSAME
92 * ..
93 * .. External Subroutines ..
94 EXTERNAL ZHER, ZSCAL, ZTRMV
95 * ..
96 * .. Intrinsic Functions ..
97 INTRINSIC DBLE, DCONJG, DIMAG
98 * ..
99 * .. Executable Statements ..
100 *
101 * Quick exit if N = 0.
102 *
103 IF( N.LE.0 ) THEN
104 RESID = ZERO
105 RETURN
106 END IF
107 *
108 * Exit with RESID = 1/EPS if ANORM = 0.
109 *
110 EPS = DLAMCH( 'Epsilon' )
111 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
112 IF( ANORM.LE.ZERO ) THEN
113 RESID = ONE / EPS
114 RETURN
115 END IF
116 *
117 * Check the imaginary parts of the diagonal elements and return with
118 * an error code if any are nonzero.
119 *
120 DO 100 J = 1, N
121 IF( DIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
122 RESID = ONE / EPS
123 RETURN
124 END IF
125 100 CONTINUE
126 *
127 * Compute the product U'*U, overwriting U.
128 *
129 IF( LSAME( UPLO, 'U' ) ) THEN
130 *
131 IF( RANK.LT.N ) THEN
132 DO 120 J = RANK + 1, N
133 DO 110 I = RANK + 1, J
134 AFAC( I, J ) = CZERO
135 110 CONTINUE
136 120 CONTINUE
137 END IF
138 *
139 DO 130 K = N, 1, -1
140 *
141 * Compute the (K,K) element of the result.
142 *
143 TR = ZDOTC( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
144 AFAC( K, K ) = TR
145 *
146 * Compute the rest of column K.
147 *
148 CALL ZTRMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
149 $ LDAFAC, AFAC( 1, K ), 1 )
150 *
151 130 CONTINUE
152 *
153 * Compute the product L*L', overwriting L.
154 *
155 ELSE
156 *
157 IF( RANK.LT.N ) THEN
158 DO 150 J = RANK + 1, N
159 DO 140 I = J, N
160 AFAC( I, J ) = CZERO
161 140 CONTINUE
162 150 CONTINUE
163 END IF
164 *
165 DO 160 K = N, 1, -1
166 * Add a multiple of column K of the factor L to each of
167 * columns K+1 through N.
168 *
169 IF( K+1.LE.N )
170 $ CALL ZHER( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
171 $ AFAC( K+1, K+1 ), LDAFAC )
172 *
173 * Scale column K by the diagonal element.
174 *
175 TC = AFAC( K, K )
176 CALL ZSCAL( N-K+1, TC, AFAC( K, K ), 1 )
177 160 CONTINUE
178 *
179 END IF
180 *
181 * Form P*L*L'*P' or P*U'*U*P'
182 *
183 IF( LSAME( UPLO, 'U' ) ) THEN
184 *
185 DO 180 J = 1, N
186 DO 170 I = 1, N
187 IF( PIV( I ).LE.PIV( J ) ) THEN
188 IF( I.LE.J ) THEN
189 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
190 ELSE
191 PERM( PIV( I ), PIV( J ) ) = DCONJG( AFAC( J, I ) )
192 END IF
193 END IF
194 170 CONTINUE
195 180 CONTINUE
196 *
197 *
198 ELSE
199 *
200 DO 200 J = 1, N
201 DO 190 I = 1, N
202 IF( PIV( I ).GE.PIV( J ) ) THEN
203 IF( I.GE.J ) THEN
204 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
205 ELSE
206 PERM( PIV( I ), PIV( J ) ) = DCONJG( AFAC( J, I ) )
207 END IF
208 END IF
209 190 CONTINUE
210 200 CONTINUE
211 *
212 END IF
213 *
214 * Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A).
215 *
216 IF( LSAME( UPLO, 'U' ) ) THEN
217 DO 220 J = 1, N
218 DO 210 I = 1, J - 1
219 PERM( I, J ) = PERM( I, J ) - A( I, J )
220 210 CONTINUE
221 PERM( J, J ) = PERM( J, J ) - DBLE( A( J, J ) )
222 220 CONTINUE
223 ELSE
224 DO 240 J = 1, N
225 PERM( J, J ) = PERM( J, J ) - DBLE( A( J, J ) )
226 DO 230 I = J + 1, N
227 PERM( I, J ) = PERM( I, J ) - A( I, J )
228 230 CONTINUE
229 240 CONTINUE
230 END IF
231 *
232 * Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
233 * ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
234 *
235 RESID = ZLANHE( '1', UPLO, N, PERM, LDAFAC, RWORK )
236 *
237 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
238 *
239 RETURN
240 *
241 * End of ZPST01
242 *
243 END
2 $ PIV, RWORK, RESID, RANK )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Craig Lucas, University of Manchester / NAG Ltd.
6 * October, 2008
7 *
8 * .. Scalar Arguments ..
9 DOUBLE PRECISION RESID
10 INTEGER LDA, LDAFAC, LDPERM, N, RANK
11 CHARACTER UPLO
12 * ..
13 * .. Array Arguments ..
14 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ),
15 $ PERM( LDPERM, * )
16 DOUBLE PRECISION RWORK( * )
17 INTEGER PIV( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * ZPST01 reconstructs an Hermitian positive semidefinite matrix A
24 * from its L or U factors and the permutation matrix P and computes
25 * the residual
26 * norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
27 * norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
28 * where EPS is the machine epsilon, L' is the conjugate transpose of L,
29 * and U' is the conjugate transpose of U.
30 *
31 * Arguments
32 * ==========
33 *
34 * UPLO (input) CHARACTER*1
35 * Specifies whether the upper or lower triangular part of the
36 * Hermitian matrix A is stored:
37 * = 'U': Upper triangular
38 * = 'L': Lower triangular
39 *
40 * N (input) INTEGER
41 * The number of rows and columns of the matrix A. N >= 0.
42 *
43 * A (input) COMPLEX*16 array, dimension (LDA,N)
44 * The original Hermitian matrix A.
45 *
46 * LDA (input) INTEGER
47 * The leading dimension of the array A. LDA >= max(1,N)
48 *
49 * AFAC (input) COMPLEX*16 array, dimension (LDAFAC,N)
50 * The factor L or U from the L*L' or U'*U
51 * factorization of A.
52 *
53 * LDAFAC (input) INTEGER
54 * The leading dimension of the array AFAC. LDAFAC >= max(1,N).
55 *
56 * PERM (output) COMPLEX*16 array, dimension (LDPERM,N)
57 * Overwritten with the reconstructed matrix, and then with the
58 * difference P*L*L'*P' - A (or P*U'*U*P' - A)
59 *
60 * LDPERM (input) INTEGER
61 * The leading dimension of the array PERM.
62 * LDAPERM >= max(1,N).
63 *
64 * PIV (input) INTEGER array, dimension (N)
65 * PIV is such that the nonzero entries are
66 * P( PIV( K ), K ) = 1.
67 *
68 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
69 *
70 * RESID (output) DOUBLE PRECISION
71 * If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
72 * If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
73 *
74 * =====================================================================
75 *
76 * .. Parameters ..
77 DOUBLE PRECISION ZERO, ONE
78 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
79 COMPLEX*16 CZERO
80 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
81 * ..
82 * .. Local Scalars ..
83 COMPLEX*16 TC
84 DOUBLE PRECISION ANORM, EPS, TR
85 INTEGER I, J, K
86 * ..
87 * .. External Functions ..
88 COMPLEX*16 ZDOTC
89 DOUBLE PRECISION DLAMCH, ZLANHE
90 LOGICAL LSAME
91 EXTERNAL ZDOTC, DLAMCH, ZLANHE, LSAME
92 * ..
93 * .. External Subroutines ..
94 EXTERNAL ZHER, ZSCAL, ZTRMV
95 * ..
96 * .. Intrinsic Functions ..
97 INTRINSIC DBLE, DCONJG, DIMAG
98 * ..
99 * .. Executable Statements ..
100 *
101 * Quick exit if N = 0.
102 *
103 IF( N.LE.0 ) THEN
104 RESID = ZERO
105 RETURN
106 END IF
107 *
108 * Exit with RESID = 1/EPS if ANORM = 0.
109 *
110 EPS = DLAMCH( 'Epsilon' )
111 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
112 IF( ANORM.LE.ZERO ) THEN
113 RESID = ONE / EPS
114 RETURN
115 END IF
116 *
117 * Check the imaginary parts of the diagonal elements and return with
118 * an error code if any are nonzero.
119 *
120 DO 100 J = 1, N
121 IF( DIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
122 RESID = ONE / EPS
123 RETURN
124 END IF
125 100 CONTINUE
126 *
127 * Compute the product U'*U, overwriting U.
128 *
129 IF( LSAME( UPLO, 'U' ) ) THEN
130 *
131 IF( RANK.LT.N ) THEN
132 DO 120 J = RANK + 1, N
133 DO 110 I = RANK + 1, J
134 AFAC( I, J ) = CZERO
135 110 CONTINUE
136 120 CONTINUE
137 END IF
138 *
139 DO 130 K = N, 1, -1
140 *
141 * Compute the (K,K) element of the result.
142 *
143 TR = ZDOTC( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
144 AFAC( K, K ) = TR
145 *
146 * Compute the rest of column K.
147 *
148 CALL ZTRMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
149 $ LDAFAC, AFAC( 1, K ), 1 )
150 *
151 130 CONTINUE
152 *
153 * Compute the product L*L', overwriting L.
154 *
155 ELSE
156 *
157 IF( RANK.LT.N ) THEN
158 DO 150 J = RANK + 1, N
159 DO 140 I = J, N
160 AFAC( I, J ) = CZERO
161 140 CONTINUE
162 150 CONTINUE
163 END IF
164 *
165 DO 160 K = N, 1, -1
166 * Add a multiple of column K of the factor L to each of
167 * columns K+1 through N.
168 *
169 IF( K+1.LE.N )
170 $ CALL ZHER( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
171 $ AFAC( K+1, K+1 ), LDAFAC )
172 *
173 * Scale column K by the diagonal element.
174 *
175 TC = AFAC( K, K )
176 CALL ZSCAL( N-K+1, TC, AFAC( K, K ), 1 )
177 160 CONTINUE
178 *
179 END IF
180 *
181 * Form P*L*L'*P' or P*U'*U*P'
182 *
183 IF( LSAME( UPLO, 'U' ) ) THEN
184 *
185 DO 180 J = 1, N
186 DO 170 I = 1, N
187 IF( PIV( I ).LE.PIV( J ) ) THEN
188 IF( I.LE.J ) THEN
189 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
190 ELSE
191 PERM( PIV( I ), PIV( J ) ) = DCONJG( AFAC( J, I ) )
192 END IF
193 END IF
194 170 CONTINUE
195 180 CONTINUE
196 *
197 *
198 ELSE
199 *
200 DO 200 J = 1, N
201 DO 190 I = 1, N
202 IF( PIV( I ).GE.PIV( J ) ) THEN
203 IF( I.GE.J ) THEN
204 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
205 ELSE
206 PERM( PIV( I ), PIV( J ) ) = DCONJG( AFAC( J, I ) )
207 END IF
208 END IF
209 190 CONTINUE
210 200 CONTINUE
211 *
212 END IF
213 *
214 * Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A).
215 *
216 IF( LSAME( UPLO, 'U' ) ) THEN
217 DO 220 J = 1, N
218 DO 210 I = 1, J - 1
219 PERM( I, J ) = PERM( I, J ) - A( I, J )
220 210 CONTINUE
221 PERM( J, J ) = PERM( J, J ) - DBLE( A( J, J ) )
222 220 CONTINUE
223 ELSE
224 DO 240 J = 1, N
225 PERM( J, J ) = PERM( J, J ) - DBLE( A( J, J ) )
226 DO 230 I = J + 1, N
227 PERM( I, J ) = PERM( I, J ) - A( I, J )
228 230 CONTINUE
229 240 CONTINUE
230 END IF
231 *
232 * Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
233 * ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
234 *
235 RESID = ZLANHE( '1', UPLO, N, PERM, LDAFAC, RWORK )
236 *
237 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
238 *
239 RETURN
240 *
241 * End of ZPST01
242 *
243 END