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