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