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