1 SUBROUTINE ZSYT01( 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 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 INTEGER IPIV( * )
15 DOUBLE PRECISION RWORK( * )
16 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZSYT01 reconstructs a complex symmetric 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 transpose of L, and U' is the transpose of U.
27 *
28 * Arguments
29 * ==========
30 *
31 * UPLO (input) CHARACTER*1
32 * Specifies whether the upper or lower triangular part of the
33 * complex symmetric matrix A is stored:
34 * = 'U': Upper triangular
35 * = 'L': Lower triangular
36 *
37 * N (input) INTEGER
38 * The number of rows and columns of the matrix A. N >= 0.
39 *
40 * A (input) COMPLEX*16 array, dimension (LDA,N)
41 * The original complex symmetric matrix A.
42 *
43 * LDA (input) INTEGER
44 * The leading dimension of the array A. LDA >= max(1,N)
45 *
46 * AFAC (input) COMPLEX*16 array, dimension (LDAFAC,N)
47 * The factored form of the matrix A. AFAC contains the block
48 * diagonal matrix D and the multipliers used to obtain the
49 * factor L or U from the block L*D*L' or U*D*U' factorization
50 * as computed by ZSYTRF.
51 *
52 * LDAFAC (input) INTEGER
53 * The leading dimension of the array AFAC. LDAFAC >= max(1,N).
54 *
55 * IPIV (input) INTEGER array, dimension (N)
56 * The pivot indices from ZSYTRF.
57 *
58 * C (workspace) COMPLEX*16 array, dimension (LDC,N)
59 *
60 * LDC (integer) INTEGER
61 * The leading dimension of the array C. LDC >= max(1,N).
62 *
63 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
64 *
65 * RESID (output) DOUBLE PRECISION
66 * If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
67 * If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
68 *
69 * =====================================================================
70 *
71 * .. Parameters ..
72 DOUBLE PRECISION ZERO, ONE
73 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
74 COMPLEX*16 CZERO, CONE
75 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
76 $ CONE = ( 1.0D+0, 0.0D+0 ) )
77 * ..
78 * .. Local Scalars ..
79 INTEGER I, INFO, J
80 DOUBLE PRECISION ANORM, EPS
81 * ..
82 * .. External Functions ..
83 LOGICAL LSAME
84 DOUBLE PRECISION DLAMCH, ZLANSY
85 EXTERNAL LSAME, DLAMCH, ZLANSY
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL ZLASET, ZLAVSY
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC DBLE
92 * ..
93 * .. Executable Statements ..
94 *
95 * Quick exit if N = 0.
96 *
97 IF( N.LE.0 ) THEN
98 RESID = ZERO
99 RETURN
100 END IF
101 *
102 * Determine EPS and the norm of A.
103 *
104 EPS = DLAMCH( 'Epsilon' )
105 ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK )
106 *
107 * Initialize C to the identity matrix.
108 *
109 CALL ZLASET( 'Full', N, N, CZERO, CONE, C, LDC )
110 *
111 * Call ZLAVSY to form the product D * U' (or D * L' ).
112 *
113 CALL ZLAVSY( UPLO, 'Transpose', 'Non-unit', N, N, AFAC, LDAFAC,
114 $ IPIV, C, LDC, INFO )
115 *
116 * Call ZLAVSY again to multiply by U (or L ).
117 *
118 CALL ZLAVSY( UPLO, 'No transpose', 'Unit', N, N, AFAC, LDAFAC,
119 $ IPIV, C, LDC, INFO )
120 *
121 * Compute the difference C - A .
122 *
123 IF( LSAME( UPLO, 'U' ) ) THEN
124 DO 20 J = 1, N
125 DO 10 I = 1, J
126 C( I, J ) = C( I, J ) - A( I, J )
127 10 CONTINUE
128 20 CONTINUE
129 ELSE
130 DO 40 J = 1, N
131 DO 30 I = J, N
132 C( I, J ) = C( I, J ) - A( I, J )
133 30 CONTINUE
134 40 CONTINUE
135 END IF
136 *
137 * Compute norm( C - A ) / ( N * norm(A) * EPS )
138 *
139 RESID = ZLANSY( '1', UPLO, N, C, LDC, RWORK )
140 *
141 IF( ANORM.LE.ZERO ) THEN
142 IF( RESID.NE.ZERO )
143 $ RESID = ONE / EPS
144 ELSE
145 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
146 END IF
147 *
148 RETURN
149 *
150 * End of ZSYT01
151 *
152 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 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 INTEGER IPIV( * )
15 DOUBLE PRECISION RWORK( * )
16 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZSYT01 reconstructs a complex symmetric 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 transpose of L, and U' is the transpose of U.
27 *
28 * Arguments
29 * ==========
30 *
31 * UPLO (input) CHARACTER*1
32 * Specifies whether the upper or lower triangular part of the
33 * complex symmetric matrix A is stored:
34 * = 'U': Upper triangular
35 * = 'L': Lower triangular
36 *
37 * N (input) INTEGER
38 * The number of rows and columns of the matrix A. N >= 0.
39 *
40 * A (input) COMPLEX*16 array, dimension (LDA,N)
41 * The original complex symmetric matrix A.
42 *
43 * LDA (input) INTEGER
44 * The leading dimension of the array A. LDA >= max(1,N)
45 *
46 * AFAC (input) COMPLEX*16 array, dimension (LDAFAC,N)
47 * The factored form of the matrix A. AFAC contains the block
48 * diagonal matrix D and the multipliers used to obtain the
49 * factor L or U from the block L*D*L' or U*D*U' factorization
50 * as computed by ZSYTRF.
51 *
52 * LDAFAC (input) INTEGER
53 * The leading dimension of the array AFAC. LDAFAC >= max(1,N).
54 *
55 * IPIV (input) INTEGER array, dimension (N)
56 * The pivot indices from ZSYTRF.
57 *
58 * C (workspace) COMPLEX*16 array, dimension (LDC,N)
59 *
60 * LDC (integer) INTEGER
61 * The leading dimension of the array C. LDC >= max(1,N).
62 *
63 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
64 *
65 * RESID (output) DOUBLE PRECISION
66 * If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
67 * If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
68 *
69 * =====================================================================
70 *
71 * .. Parameters ..
72 DOUBLE PRECISION ZERO, ONE
73 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
74 COMPLEX*16 CZERO, CONE
75 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
76 $ CONE = ( 1.0D+0, 0.0D+0 ) )
77 * ..
78 * .. Local Scalars ..
79 INTEGER I, INFO, J
80 DOUBLE PRECISION ANORM, EPS
81 * ..
82 * .. External Functions ..
83 LOGICAL LSAME
84 DOUBLE PRECISION DLAMCH, ZLANSY
85 EXTERNAL LSAME, DLAMCH, ZLANSY
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL ZLASET, ZLAVSY
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC DBLE
92 * ..
93 * .. Executable Statements ..
94 *
95 * Quick exit if N = 0.
96 *
97 IF( N.LE.0 ) THEN
98 RESID = ZERO
99 RETURN
100 END IF
101 *
102 * Determine EPS and the norm of A.
103 *
104 EPS = DLAMCH( 'Epsilon' )
105 ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK )
106 *
107 * Initialize C to the identity matrix.
108 *
109 CALL ZLASET( 'Full', N, N, CZERO, CONE, C, LDC )
110 *
111 * Call ZLAVSY to form the product D * U' (or D * L' ).
112 *
113 CALL ZLAVSY( UPLO, 'Transpose', 'Non-unit', N, N, AFAC, LDAFAC,
114 $ IPIV, C, LDC, INFO )
115 *
116 * Call ZLAVSY again to multiply by U (or L ).
117 *
118 CALL ZLAVSY( UPLO, 'No transpose', 'Unit', N, N, AFAC, LDAFAC,
119 $ IPIV, C, LDC, INFO )
120 *
121 * Compute the difference C - A .
122 *
123 IF( LSAME( UPLO, 'U' ) ) THEN
124 DO 20 J = 1, N
125 DO 10 I = 1, J
126 C( I, J ) = C( I, J ) - A( I, J )
127 10 CONTINUE
128 20 CONTINUE
129 ELSE
130 DO 40 J = 1, N
131 DO 30 I = J, N
132 C( I, J ) = C( I, J ) - A( I, J )
133 30 CONTINUE
134 40 CONTINUE
135 END IF
136 *
137 * Compute norm( C - A ) / ( N * norm(A) * EPS )
138 *
139 RESID = ZLANSY( '1', UPLO, N, C, LDC, RWORK )
140 *
141 IF( ANORM.LE.ZERO ) THEN
142 IF( RESID.NE.ZERO )
143 $ RESID = ONE / EPS
144 ELSE
145 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
146 END IF
147 *
148 RETURN
149 *
150 * End of ZSYT01
151 *
152 END