1 SUBROUTINE ZPOT01( UPLO, N, A, LDA, AFAC, LDAFAC, 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 LDA, LDAFAC, N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION RWORK( * )
14 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZPOT01 reconstructs a Hermitian positive definite matrix A from
21 * 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 L,
25 * 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 (LDA,N)
40 * The original Hermitian matrix A.
41 *
42 * LDA (input) INTEGER
43 * The leading dimension of the array A. LDA >= max(1,N)
44 *
45 * AFAC (input/output) COMPLEX*16 array, dimension (LDAFAC,N)
46 * On entry, the factor L or U from the L*L' or U'*U
47 * factorization of A.
48 * Overwritten with the reconstructed matrix, and then with the
49 * difference L*L' - A (or U'*U - A).
50 *
51 * LDAFAC (input) INTEGER
52 * The leading dimension of the array AFAC. LDAFAC >= max(1,N).
53 *
54 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
55 *
56 * RESID (output) DOUBLE PRECISION
57 * If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
58 * If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ZERO, ONE
64 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
65 * ..
66 * .. Local Scalars ..
67 INTEGER I, J, K
68 DOUBLE PRECISION ANORM, EPS, TR
69 COMPLEX*16 TC
70 * ..
71 * .. External Functions ..
72 LOGICAL LSAME
73 DOUBLE PRECISION DLAMCH, ZLANHE
74 COMPLEX*16 ZDOTC
75 EXTERNAL LSAME, DLAMCH, ZLANHE, ZDOTC
76 * ..
77 * .. External Subroutines ..
78 EXTERNAL ZHER, ZSCAL, ZTRMV
79 * ..
80 * .. Intrinsic Functions ..
81 INTRINSIC DBLE, DIMAG
82 * ..
83 * .. Executable Statements ..
84 *
85 * Quick exit if N = 0.
86 *
87 IF( N.LE.0 ) THEN
88 RESID = ZERO
89 RETURN
90 END IF
91 *
92 * Exit with RESID = 1/EPS if ANORM = 0.
93 *
94 EPS = DLAMCH( 'Epsilon' )
95 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
96 IF( ANORM.LE.ZERO ) THEN
97 RESID = ONE / EPS
98 RETURN
99 END IF
100 *
101 * Check the imaginary parts of the diagonal elements and return with
102 * an error code if any are nonzero.
103 *
104 DO 10 J = 1, N
105 IF( DIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
106 RESID = ONE / EPS
107 RETURN
108 END IF
109 10 CONTINUE
110 *
111 * Compute the product U'*U, overwriting U.
112 *
113 IF( LSAME( UPLO, 'U' ) ) THEN
114 DO 20 K = N, 1, -1
115 *
116 * Compute the (K,K) element of the result.
117 *
118 TR = ZDOTC( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
119 AFAC( K, K ) = TR
120 *
121 * Compute the rest of column K.
122 *
123 CALL ZTRMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
124 $ LDAFAC, AFAC( 1, K ), 1 )
125 *
126 20 CONTINUE
127 *
128 * Compute the product L*L', overwriting L.
129 *
130 ELSE
131 DO 30 K = N, 1, -1
132 *
133 * Add a multiple of column K of the factor L to each of
134 * columns K+1 through N.
135 *
136 IF( K+1.LE.N )
137 $ CALL ZHER( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
138 $ AFAC( K+1, K+1 ), LDAFAC )
139 *
140 * Scale column K by the diagonal element.
141 *
142 TC = AFAC( K, K )
143 CALL ZSCAL( N-K+1, TC, AFAC( K, K ), 1 )
144 *
145 30 CONTINUE
146 END IF
147 *
148 * Compute the difference L*L' - A (or U'*U - A).
149 *
150 IF( LSAME( UPLO, 'U' ) ) THEN
151 DO 50 J = 1, N
152 DO 40 I = 1, J - 1
153 AFAC( I, J ) = AFAC( I, J ) - A( I, J )
154 40 CONTINUE
155 AFAC( J, J ) = AFAC( J, J ) - DBLE( A( J, J ) )
156 50 CONTINUE
157 ELSE
158 DO 70 J = 1, N
159 AFAC( J, J ) = AFAC( J, J ) - DBLE( A( J, J ) )
160 DO 60 I = J + 1, N
161 AFAC( I, J ) = AFAC( I, J ) - A( I, J )
162 60 CONTINUE
163 70 CONTINUE
164 END IF
165 *
166 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
167 *
168 RESID = ZLANHE( '1', UPLO, N, AFAC, LDAFAC, RWORK )
169 *
170 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
171 *
172 RETURN
173 *
174 * End of ZPOT01
175 *
176 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 LDA, LDAFAC, N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION RWORK( * )
14 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZPOT01 reconstructs a Hermitian positive definite matrix A from
21 * 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 L,
25 * 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 (LDA,N)
40 * The original Hermitian matrix A.
41 *
42 * LDA (input) INTEGER
43 * The leading dimension of the array A. LDA >= max(1,N)
44 *
45 * AFAC (input/output) COMPLEX*16 array, dimension (LDAFAC,N)
46 * On entry, the factor L or U from the L*L' or U'*U
47 * factorization of A.
48 * Overwritten with the reconstructed matrix, and then with the
49 * difference L*L' - A (or U'*U - A).
50 *
51 * LDAFAC (input) INTEGER
52 * The leading dimension of the array AFAC. LDAFAC >= max(1,N).
53 *
54 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
55 *
56 * RESID (output) DOUBLE PRECISION
57 * If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
58 * If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ZERO, ONE
64 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
65 * ..
66 * .. Local Scalars ..
67 INTEGER I, J, K
68 DOUBLE PRECISION ANORM, EPS, TR
69 COMPLEX*16 TC
70 * ..
71 * .. External Functions ..
72 LOGICAL LSAME
73 DOUBLE PRECISION DLAMCH, ZLANHE
74 COMPLEX*16 ZDOTC
75 EXTERNAL LSAME, DLAMCH, ZLANHE, ZDOTC
76 * ..
77 * .. External Subroutines ..
78 EXTERNAL ZHER, ZSCAL, ZTRMV
79 * ..
80 * .. Intrinsic Functions ..
81 INTRINSIC DBLE, DIMAG
82 * ..
83 * .. Executable Statements ..
84 *
85 * Quick exit if N = 0.
86 *
87 IF( N.LE.0 ) THEN
88 RESID = ZERO
89 RETURN
90 END IF
91 *
92 * Exit with RESID = 1/EPS if ANORM = 0.
93 *
94 EPS = DLAMCH( 'Epsilon' )
95 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
96 IF( ANORM.LE.ZERO ) THEN
97 RESID = ONE / EPS
98 RETURN
99 END IF
100 *
101 * Check the imaginary parts of the diagonal elements and return with
102 * an error code if any are nonzero.
103 *
104 DO 10 J = 1, N
105 IF( DIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
106 RESID = ONE / EPS
107 RETURN
108 END IF
109 10 CONTINUE
110 *
111 * Compute the product U'*U, overwriting U.
112 *
113 IF( LSAME( UPLO, 'U' ) ) THEN
114 DO 20 K = N, 1, -1
115 *
116 * Compute the (K,K) element of the result.
117 *
118 TR = ZDOTC( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
119 AFAC( K, K ) = TR
120 *
121 * Compute the rest of column K.
122 *
123 CALL ZTRMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
124 $ LDAFAC, AFAC( 1, K ), 1 )
125 *
126 20 CONTINUE
127 *
128 * Compute the product L*L', overwriting L.
129 *
130 ELSE
131 DO 30 K = N, 1, -1
132 *
133 * Add a multiple of column K of the factor L to each of
134 * columns K+1 through N.
135 *
136 IF( K+1.LE.N )
137 $ CALL ZHER( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
138 $ AFAC( K+1, K+1 ), LDAFAC )
139 *
140 * Scale column K by the diagonal element.
141 *
142 TC = AFAC( K, K )
143 CALL ZSCAL( N-K+1, TC, AFAC( K, K ), 1 )
144 *
145 30 CONTINUE
146 END IF
147 *
148 * Compute the difference L*L' - A (or U'*U - A).
149 *
150 IF( LSAME( UPLO, 'U' ) ) THEN
151 DO 50 J = 1, N
152 DO 40 I = 1, J - 1
153 AFAC( I, J ) = AFAC( I, J ) - A( I, J )
154 40 CONTINUE
155 AFAC( J, J ) = AFAC( J, J ) - DBLE( A( J, J ) )
156 50 CONTINUE
157 ELSE
158 DO 70 J = 1, N
159 AFAC( J, J ) = AFAC( J, J ) - DBLE( A( J, J ) )
160 DO 60 I = J + 1, N
161 AFAC( I, J ) = AFAC( I, J ) - A( I, J )
162 60 CONTINUE
163 70 CONTINUE
164 END IF
165 *
166 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
167 *
168 RESID = ZLANHE( '1', UPLO, N, AFAC, LDAFAC, RWORK )
169 *
170 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
171 *
172 RETURN
173 *
174 * End of ZPOT01
175 *
176 END