1 SUBROUTINE ZPPT01( UPLO, N, A, AFAC, RWORK, RESID )
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 CHARACTER UPLO
9 INTEGER N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION RWORK( * )
14 COMPLEX*16 A( * ), AFAC( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZPPT01 reconstructs a Hermitian positive definite packed matrix A
21 * from its L*L' or U'*U factorization and computes the residual
22 * norm( L*L' - A ) / ( N * norm(A) * EPS ) or
23 * norm( U'*U - A ) / ( N * norm(A) * EPS ),
24 * where EPS is the machine epsilon, L' is the conjugate transpose of
25 * L, and U' is the conjugate transpose of U.
26 *
27 * Arguments
28 * ==========
29 *
30 * UPLO (input) CHARACTER*1
31 * Specifies whether the upper or lower triangular part of the
32 * Hermitian matrix A is stored:
33 * = 'U': Upper triangular
34 * = 'L': Lower triangular
35 *
36 * N (input) INTEGER
37 * The number of rows and columns of the matrix A. N >= 0.
38 *
39 * A (input) COMPLEX*16 array, dimension (N*(N+1)/2)
40 * The original Hermitian matrix A, stored as a packed
41 * triangular matrix.
42 *
43 * AFAC (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
44 * On entry, the factor L or U from the L*L' or U'*U
45 * factorization of A, stored as a packed triangular matrix.
46 * Overwritten with the reconstructed matrix, and then with the
47 * difference L*L' - A (or U'*U - A).
48 *
49 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
50 *
51 * RESID (output) DOUBLE PRECISION
52 * If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
53 * If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
54 *
55 * =====================================================================
56 *
57 * .. Parameters ..
58 DOUBLE PRECISION ZERO, ONE
59 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
60 * ..
61 * .. Local Scalars ..
62 INTEGER I, K, KC
63 DOUBLE PRECISION ANORM, EPS, TR
64 COMPLEX*16 TC
65 * ..
66 * .. External Functions ..
67 LOGICAL LSAME
68 DOUBLE PRECISION DLAMCH, ZLANHP
69 COMPLEX*16 ZDOTC
70 EXTERNAL LSAME, DLAMCH, ZLANHP, ZDOTC
71 * ..
72 * .. External Subroutines ..
73 EXTERNAL ZHPR, ZSCAL, ZTPMV
74 * ..
75 * .. Intrinsic Functions ..
76 INTRINSIC DBLE, DIMAG
77 * ..
78 * .. Executable Statements ..
79 *
80 * Quick exit if N = 0
81 *
82 IF( N.LE.0 ) THEN
83 RESID = ZERO
84 RETURN
85 END IF
86 *
87 * Exit with RESID = 1/EPS if ANORM = 0.
88 *
89 EPS = DLAMCH( 'Epsilon' )
90 ANORM = ZLANHP( '1', UPLO, N, A, RWORK )
91 IF( ANORM.LE.ZERO ) THEN
92 RESID = ONE / EPS
93 RETURN
94 END IF
95 *
96 * Check the imaginary parts of the diagonal elements and return with
97 * an error code if any are nonzero.
98 *
99 KC = 1
100 IF( LSAME( UPLO, 'U' ) ) THEN
101 DO 10 K = 1, N
102 IF( DIMAG( AFAC( KC ) ).NE.ZERO ) THEN
103 RESID = ONE / EPS
104 RETURN
105 END IF
106 KC = KC + K + 1
107 10 CONTINUE
108 ELSE
109 DO 20 K = 1, N
110 IF( DIMAG( AFAC( KC ) ).NE.ZERO ) THEN
111 RESID = ONE / EPS
112 RETURN
113 END IF
114 KC = KC + N - K + 1
115 20 CONTINUE
116 END IF
117 *
118 * Compute the product U'*U, overwriting U.
119 *
120 IF( LSAME( UPLO, 'U' ) ) THEN
121 KC = ( N*( N-1 ) ) / 2 + 1
122 DO 30 K = N, 1, -1
123 *
124 * Compute the (K,K) element of the result.
125 *
126 TR = ZDOTC( K, AFAC( KC ), 1, AFAC( KC ), 1 )
127 AFAC( KC+K-1 ) = TR
128 *
129 * Compute the rest of column K.
130 *
131 IF( K.GT.1 ) THEN
132 CALL ZTPMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
133 $ AFAC( KC ), 1 )
134 KC = KC - ( K-1 )
135 END IF
136 30 CONTINUE
137 *
138 * Compute the difference L*L' - A
139 *
140 KC = 1
141 DO 50 K = 1, N
142 DO 40 I = 1, K - 1
143 AFAC( KC+I-1 ) = AFAC( KC+I-1 ) - A( KC+I-1 )
144 40 CONTINUE
145 AFAC( KC+K-1 ) = AFAC( KC+K-1 ) - DBLE( A( KC+K-1 ) )
146 KC = KC + K
147 50 CONTINUE
148 *
149 * Compute the product L*L', overwriting L.
150 *
151 ELSE
152 KC = ( N*( N+1 ) ) / 2
153 DO 60 K = N, 1, -1
154 *
155 * Add a multiple of column K of the factor L to each of
156 * columns K+1 through N.
157 *
158 IF( K.LT.N )
159 $ CALL ZHPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1,
160 $ AFAC( KC+N-K+1 ) )
161 *
162 * Scale column K by the diagonal element.
163 *
164 TC = AFAC( KC )
165 CALL ZSCAL( N-K+1, TC, AFAC( KC ), 1 )
166 *
167 KC = KC - ( N-K+2 )
168 60 CONTINUE
169 *
170 * Compute the difference U'*U - A
171 *
172 KC = 1
173 DO 80 K = 1, N
174 AFAC( KC ) = AFAC( KC ) - DBLE( A( KC ) )
175 DO 70 I = K + 1, N
176 AFAC( KC+I-K ) = AFAC( KC+I-K ) - A( KC+I-K )
177 70 CONTINUE
178 KC = KC + N - K + 1
179 80 CONTINUE
180 END IF
181 *
182 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
183 *
184 RESID = ZLANHP( '1', UPLO, N, AFAC, RWORK )
185 *
186 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
187 *
188 RETURN
189 *
190 * End of ZPPT01
191 *
192 END
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 CHARACTER UPLO
9 INTEGER N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION RWORK( * )
14 COMPLEX*16 A( * ), AFAC( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZPPT01 reconstructs a Hermitian positive definite packed matrix A
21 * from its L*L' or U'*U factorization and computes the residual
22 * norm( L*L' - A ) / ( N * norm(A) * EPS ) or
23 * norm( U'*U - A ) / ( N * norm(A) * EPS ),
24 * where EPS is the machine epsilon, L' is the conjugate transpose of
25 * L, and U' is the conjugate transpose of U.
26 *
27 * Arguments
28 * ==========
29 *
30 * UPLO (input) CHARACTER*1
31 * Specifies whether the upper or lower triangular part of the
32 * Hermitian matrix A is stored:
33 * = 'U': Upper triangular
34 * = 'L': Lower triangular
35 *
36 * N (input) INTEGER
37 * The number of rows and columns of the matrix A. N >= 0.
38 *
39 * A (input) COMPLEX*16 array, dimension (N*(N+1)/2)
40 * The original Hermitian matrix A, stored as a packed
41 * triangular matrix.
42 *
43 * AFAC (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
44 * On entry, the factor L or U from the L*L' or U'*U
45 * factorization of A, stored as a packed triangular matrix.
46 * Overwritten with the reconstructed matrix, and then with the
47 * difference L*L' - A (or U'*U - A).
48 *
49 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
50 *
51 * RESID (output) DOUBLE PRECISION
52 * If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
53 * If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
54 *
55 * =====================================================================
56 *
57 * .. Parameters ..
58 DOUBLE PRECISION ZERO, ONE
59 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
60 * ..
61 * .. Local Scalars ..
62 INTEGER I, K, KC
63 DOUBLE PRECISION ANORM, EPS, TR
64 COMPLEX*16 TC
65 * ..
66 * .. External Functions ..
67 LOGICAL LSAME
68 DOUBLE PRECISION DLAMCH, ZLANHP
69 COMPLEX*16 ZDOTC
70 EXTERNAL LSAME, DLAMCH, ZLANHP, ZDOTC
71 * ..
72 * .. External Subroutines ..
73 EXTERNAL ZHPR, ZSCAL, ZTPMV
74 * ..
75 * .. Intrinsic Functions ..
76 INTRINSIC DBLE, DIMAG
77 * ..
78 * .. Executable Statements ..
79 *
80 * Quick exit if N = 0
81 *
82 IF( N.LE.0 ) THEN
83 RESID = ZERO
84 RETURN
85 END IF
86 *
87 * Exit with RESID = 1/EPS if ANORM = 0.
88 *
89 EPS = DLAMCH( 'Epsilon' )
90 ANORM = ZLANHP( '1', UPLO, N, A, RWORK )
91 IF( ANORM.LE.ZERO ) THEN
92 RESID = ONE / EPS
93 RETURN
94 END IF
95 *
96 * Check the imaginary parts of the diagonal elements and return with
97 * an error code if any are nonzero.
98 *
99 KC = 1
100 IF( LSAME( UPLO, 'U' ) ) THEN
101 DO 10 K = 1, N
102 IF( DIMAG( AFAC( KC ) ).NE.ZERO ) THEN
103 RESID = ONE / EPS
104 RETURN
105 END IF
106 KC = KC + K + 1
107 10 CONTINUE
108 ELSE
109 DO 20 K = 1, N
110 IF( DIMAG( AFAC( KC ) ).NE.ZERO ) THEN
111 RESID = ONE / EPS
112 RETURN
113 END IF
114 KC = KC + N - K + 1
115 20 CONTINUE
116 END IF
117 *
118 * Compute the product U'*U, overwriting U.
119 *
120 IF( LSAME( UPLO, 'U' ) ) THEN
121 KC = ( N*( N-1 ) ) / 2 + 1
122 DO 30 K = N, 1, -1
123 *
124 * Compute the (K,K) element of the result.
125 *
126 TR = ZDOTC( K, AFAC( KC ), 1, AFAC( KC ), 1 )
127 AFAC( KC+K-1 ) = TR
128 *
129 * Compute the rest of column K.
130 *
131 IF( K.GT.1 ) THEN
132 CALL ZTPMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
133 $ AFAC( KC ), 1 )
134 KC = KC - ( K-1 )
135 END IF
136 30 CONTINUE
137 *
138 * Compute the difference L*L' - A
139 *
140 KC = 1
141 DO 50 K = 1, N
142 DO 40 I = 1, K - 1
143 AFAC( KC+I-1 ) = AFAC( KC+I-1 ) - A( KC+I-1 )
144 40 CONTINUE
145 AFAC( KC+K-1 ) = AFAC( KC+K-1 ) - DBLE( A( KC+K-1 ) )
146 KC = KC + K
147 50 CONTINUE
148 *
149 * Compute the product L*L', overwriting L.
150 *
151 ELSE
152 KC = ( N*( N+1 ) ) / 2
153 DO 60 K = N, 1, -1
154 *
155 * Add a multiple of column K of the factor L to each of
156 * columns K+1 through N.
157 *
158 IF( K.LT.N )
159 $ CALL ZHPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1,
160 $ AFAC( KC+N-K+1 ) )
161 *
162 * Scale column K by the diagonal element.
163 *
164 TC = AFAC( KC )
165 CALL ZSCAL( N-K+1, TC, AFAC( KC ), 1 )
166 *
167 KC = KC - ( N-K+2 )
168 60 CONTINUE
169 *
170 * Compute the difference U'*U - A
171 *
172 KC = 1
173 DO 80 K = 1, N
174 AFAC( KC ) = AFAC( KC ) - DBLE( A( KC ) )
175 DO 70 I = K + 1, N
176 AFAC( KC+I-K ) = AFAC( KC+I-K ) - A( KC+I-K )
177 70 CONTINUE
178 KC = KC + N - K + 1
179 80 CONTINUE
180 END IF
181 *
182 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
183 *
184 RESID = ZLANHP( '1', UPLO, N, AFAC, RWORK )
185 *
186 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
187 *
188 RETURN
189 *
190 * End of ZPPT01
191 *
192 END