1 SUBROUTINE CHET01( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC,
2 $ RWORK, RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER LDA, LDAFAC, LDC, N
11 REAL RESID
12 * ..
13 * .. Array Arguments ..
14 INTEGER IPIV( * )
15 REAL RWORK( * )
16 COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * CHET01 reconstructs a Hermitian indefinite matrix A from its
23 * block L*D*L' or U*D*U' factorization and computes the residual
24 * norm( C - A ) / ( N * norm(A) * EPS ),
25 * where C is the reconstructed matrix, EPS is the machine epsilon,
26 * L' is the conjugate transpose of L, and U' is the conjugate transpose
27 * of U.
28 *
29 * Arguments
30 * ==========
31 *
32 * UPLO (input) CHARACTER*1
33 * Specifies whether the upper or lower triangular part of the
34 * Hermitian matrix A is stored:
35 * = 'U': Upper triangular
36 * = 'L': Lower triangular
37 *
38 * N (input) INTEGER
39 * The number of rows and columns of the matrix A. N >= 0.
40 *
41 * A (input) COMPLEX array, dimension (LDA,N)
42 * The original Hermitian matrix A.
43 *
44 * LDA (input) INTEGER
45 * The leading dimension of the array A. LDA >= max(1,N)
46 *
47 * AFAC (input) COMPLEX array, dimension (LDAFAC,N)
48 * The factored form of the matrix A. AFAC contains the block
49 * diagonal matrix D and the multipliers used to obtain the
50 * factor L or U from the block L*D*L' or U*D*U' factorization
51 * as computed by CHETRF.
52 *
53 * LDAFAC (input) INTEGER
54 * The leading dimension of the array AFAC. LDAFAC >= max(1,N).
55 *
56 * IPIV (input) INTEGER array, dimension (N)
57 * The pivot indices from CHETRF.
58 *
59 * C (workspace) COMPLEX array, dimension (LDC,N)
60 *
61 * LDC (integer) INTEGER
62 * The leading dimension of the array C. LDC >= max(1,N).
63 *
64 * RWORK (workspace) REAL array, dimension (N)
65 *
66 * RESID (output) REAL
67 * If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
68 * If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
69 *
70 * =====================================================================
71 *
72 * .. Parameters ..
73 REAL ZERO, ONE
74 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
75 COMPLEX CZERO, CONE
76 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
77 $ CONE = ( 1.0E+0, 0.0E+0 ) )
78 * ..
79 * .. Local Scalars ..
80 INTEGER I, INFO, J
81 REAL ANORM, EPS
82 * ..
83 * .. External Functions ..
84 LOGICAL LSAME
85 REAL CLANHE, SLAMCH
86 EXTERNAL LSAME, CLANHE, SLAMCH
87 * ..
88 * .. External Subroutines ..
89 EXTERNAL CLAVHE, CLASET
90 * ..
91 * .. Intrinsic Functions ..
92 INTRINSIC AIMAG, REAL
93 * ..
94 * .. Executable Statements ..
95 *
96 * Quick exit if N = 0.
97 *
98 IF( N.LE.0 ) THEN
99 RESID = ZERO
100 RETURN
101 END IF
102 *
103 * Determine EPS and the norm of A.
104 *
105 EPS = SLAMCH( 'Epsilon' )
106 ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
107 *
108 * Check the imaginary parts of the diagonal elements and return with
109 * an error code if any are nonzero.
110 *
111 DO 10 J = 1, N
112 IF( AIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
113 RESID = ONE / EPS
114 RETURN
115 END IF
116 10 CONTINUE
117 *
118 * Initialize C to the identity matrix.
119 *
120 CALL CLASET( 'Full', N, N, CZERO, CONE, C, LDC )
121 *
122 * Call CLAVHE to form the product D * U' (or D * L' ).
123 *
124 CALL CLAVHE( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC, LDAFAC,
125 $ IPIV, C, LDC, INFO )
126 *
127 * Call CLAVHE again to multiply by U (or L ).
128 *
129 CALL CLAVHE( UPLO, 'No transpose', 'Unit', N, N, AFAC, LDAFAC,
130 $ IPIV, C, LDC, INFO )
131 *
132 * Compute the difference C - A .
133 *
134 IF( LSAME( UPLO, 'U' ) ) THEN
135 DO 30 J = 1, N
136 DO 20 I = 1, J - 1
137 C( I, J ) = C( I, J ) - A( I, J )
138 20 CONTINUE
139 C( J, J ) = C( J, J ) - REAL( A( J, J ) )
140 30 CONTINUE
141 ELSE
142 DO 50 J = 1, N
143 C( J, J ) = C( J, J ) - REAL( A( J, J ) )
144 DO 40 I = J + 1, N
145 C( I, J ) = C( I, J ) - A( I, J )
146 40 CONTINUE
147 50 CONTINUE
148 END IF
149 *
150 * Compute norm( C - A ) / ( N * norm(A) * EPS )
151 *
152 RESID = CLANHE( '1', UPLO, N, C, LDC, RWORK )
153 *
154 IF( ANORM.LE.ZERO ) THEN
155 IF( RESID.NE.ZERO )
156 $ RESID = ONE / EPS
157 ELSE
158 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
159 END IF
160 *
161 RETURN
162 *
163 * End of CHET01
164 *
165 END
2 $ RWORK, RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER LDA, LDAFAC, LDC, N
11 REAL RESID
12 * ..
13 * .. Array Arguments ..
14 INTEGER IPIV( * )
15 REAL RWORK( * )
16 COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * CHET01 reconstructs a Hermitian indefinite matrix A from its
23 * block L*D*L' or U*D*U' factorization and computes the residual
24 * norm( C - A ) / ( N * norm(A) * EPS ),
25 * where C is the reconstructed matrix, EPS is the machine epsilon,
26 * L' is the conjugate transpose of L, and U' is the conjugate transpose
27 * of U.
28 *
29 * Arguments
30 * ==========
31 *
32 * UPLO (input) CHARACTER*1
33 * Specifies whether the upper or lower triangular part of the
34 * Hermitian matrix A is stored:
35 * = 'U': Upper triangular
36 * = 'L': Lower triangular
37 *
38 * N (input) INTEGER
39 * The number of rows and columns of the matrix A. N >= 0.
40 *
41 * A (input) COMPLEX array, dimension (LDA,N)
42 * The original Hermitian matrix A.
43 *
44 * LDA (input) INTEGER
45 * The leading dimension of the array A. LDA >= max(1,N)
46 *
47 * AFAC (input) COMPLEX array, dimension (LDAFAC,N)
48 * The factored form of the matrix A. AFAC contains the block
49 * diagonal matrix D and the multipliers used to obtain the
50 * factor L or U from the block L*D*L' or U*D*U' factorization
51 * as computed by CHETRF.
52 *
53 * LDAFAC (input) INTEGER
54 * The leading dimension of the array AFAC. LDAFAC >= max(1,N).
55 *
56 * IPIV (input) INTEGER array, dimension (N)
57 * The pivot indices from CHETRF.
58 *
59 * C (workspace) COMPLEX array, dimension (LDC,N)
60 *
61 * LDC (integer) INTEGER
62 * The leading dimension of the array C. LDC >= max(1,N).
63 *
64 * RWORK (workspace) REAL array, dimension (N)
65 *
66 * RESID (output) REAL
67 * If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
68 * If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
69 *
70 * =====================================================================
71 *
72 * .. Parameters ..
73 REAL ZERO, ONE
74 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
75 COMPLEX CZERO, CONE
76 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
77 $ CONE = ( 1.0E+0, 0.0E+0 ) )
78 * ..
79 * .. Local Scalars ..
80 INTEGER I, INFO, J
81 REAL ANORM, EPS
82 * ..
83 * .. External Functions ..
84 LOGICAL LSAME
85 REAL CLANHE, SLAMCH
86 EXTERNAL LSAME, CLANHE, SLAMCH
87 * ..
88 * .. External Subroutines ..
89 EXTERNAL CLAVHE, CLASET
90 * ..
91 * .. Intrinsic Functions ..
92 INTRINSIC AIMAG, REAL
93 * ..
94 * .. Executable Statements ..
95 *
96 * Quick exit if N = 0.
97 *
98 IF( N.LE.0 ) THEN
99 RESID = ZERO
100 RETURN
101 END IF
102 *
103 * Determine EPS and the norm of A.
104 *
105 EPS = SLAMCH( 'Epsilon' )
106 ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
107 *
108 * Check the imaginary parts of the diagonal elements and return with
109 * an error code if any are nonzero.
110 *
111 DO 10 J = 1, N
112 IF( AIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
113 RESID = ONE / EPS
114 RETURN
115 END IF
116 10 CONTINUE
117 *
118 * Initialize C to the identity matrix.
119 *
120 CALL CLASET( 'Full', N, N, CZERO, CONE, C, LDC )
121 *
122 * Call CLAVHE to form the product D * U' (or D * L' ).
123 *
124 CALL CLAVHE( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC, LDAFAC,
125 $ IPIV, C, LDC, INFO )
126 *
127 * Call CLAVHE again to multiply by U (or L ).
128 *
129 CALL CLAVHE( UPLO, 'No transpose', 'Unit', N, N, AFAC, LDAFAC,
130 $ IPIV, C, LDC, INFO )
131 *
132 * Compute the difference C - A .
133 *
134 IF( LSAME( UPLO, 'U' ) ) THEN
135 DO 30 J = 1, N
136 DO 20 I = 1, J - 1
137 C( I, J ) = C( I, J ) - A( I, J )
138 20 CONTINUE
139 C( J, J ) = C( J, J ) - REAL( A( J, J ) )
140 30 CONTINUE
141 ELSE
142 DO 50 J = 1, N
143 C( J, J ) = C( J, J ) - REAL( A( J, J ) )
144 DO 40 I = J + 1, N
145 C( I, J ) = C( I, J ) - A( I, J )
146 40 CONTINUE
147 50 CONTINUE
148 END IF
149 *
150 * Compute norm( C - A ) / ( N * norm(A) * EPS )
151 *
152 RESID = CLANHE( '1', UPLO, N, C, LDC, RWORK )
153 *
154 IF( ANORM.LE.ZERO ) THEN
155 IF( RESID.NE.ZERO )
156 $ RESID = ONE / EPS
157 ELSE
158 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
159 END IF
160 *
161 RETURN
162 *
163 * End of CHET01
164 *
165 END