1 SUBROUTINE SPPT01( UPLO, N, A, AFAC, 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 N
10 REAL RESID
11 * ..
12 * .. Array Arguments ..
13 REAL A( * ), AFAC( * ), RWORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * SPPT01 reconstructs a symmetric positive definite packed matrix A
20 * from its L*L' or U'*U factorization and computes the residual
21 * norm( L*L' - A ) / ( N * norm(A) * EPS ) or
22 * norm( U'*U - A ) / ( N * norm(A) * EPS ),
23 * where 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/output) REAL array, dimension (N*(N+1)/2)
42 * On entry, the factor L or U from the L*L' or U'*U
43 * factorization of A, stored as a packed triangular matrix.
44 * Overwritten with the reconstructed matrix, and then with the
45 * difference L*L' - A (or U'*U - A).
46 *
47 * RWORK (workspace) REAL array, dimension (N)
48 *
49 * RESID (output) REAL
50 * If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
51 * If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
52 *
53 * =====================================================================
54 *
55 * .. Parameters ..
56 REAL ZERO, ONE
57 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
58 * ..
59 * .. Local Scalars ..
60 INTEGER I, K, KC, NPP
61 REAL ANORM, EPS, T
62 * ..
63 * .. External Functions ..
64 LOGICAL LSAME
65 REAL SDOT, SLAMCH, SLANSP
66 EXTERNAL LSAME, SDOT, SLAMCH, SLANSP
67 * ..
68 * .. External Subroutines ..
69 EXTERNAL SSCAL, SSPR, STPMV
70 * ..
71 * .. Intrinsic Functions ..
72 INTRINSIC REAL
73 * ..
74 * .. Executable Statements ..
75 *
76 * Quick exit if N = 0
77 *
78 IF( N.LE.0 ) THEN
79 RESID = ZERO
80 RETURN
81 END IF
82 *
83 * Exit with RESID = 1/EPS if ANORM = 0.
84 *
85 EPS = SLAMCH( 'Epsilon' )
86 ANORM = SLANSP( '1', UPLO, N, A, RWORK )
87 IF( ANORM.LE.ZERO ) THEN
88 RESID = ONE / EPS
89 RETURN
90 END IF
91 *
92 * Compute the product U'*U, overwriting U.
93 *
94 IF( LSAME( UPLO, 'U' ) ) THEN
95 KC = ( N*( N-1 ) ) / 2 + 1
96 DO 10 K = N, 1, -1
97 *
98 * Compute the (K,K) element of the result.
99 *
100 T = SDOT( K, AFAC( KC ), 1, AFAC( KC ), 1 )
101 AFAC( KC+K-1 ) = T
102 *
103 * Compute the rest of column K.
104 *
105 IF( K.GT.1 ) THEN
106 CALL STPMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
107 $ AFAC( KC ), 1 )
108 KC = KC - ( K-1 )
109 END IF
110 10 CONTINUE
111 *
112 * Compute the product L*L', overwriting L.
113 *
114 ELSE
115 KC = ( N*( N+1 ) ) / 2
116 DO 20 K = N, 1, -1
117 *
118 * Add a multiple of column K of the factor L to each of
119 * columns K+1 through N.
120 *
121 IF( K.LT.N )
122 $ CALL SSPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1,
123 $ AFAC( KC+N-K+1 ) )
124 *
125 * Scale column K by the diagonal element.
126 *
127 T = AFAC( KC )
128 CALL SSCAL( N-K+1, T, AFAC( KC ), 1 )
129 *
130 KC = KC - ( N-K+2 )
131 20 CONTINUE
132 END IF
133 *
134 * Compute the difference L*L' - A (or U'*U - A).
135 *
136 NPP = N*( N+1 ) / 2
137 DO 30 I = 1, NPP
138 AFAC( I ) = AFAC( I ) - A( I )
139 30 CONTINUE
140 *
141 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
142 *
143 RESID = SLANSP( '1', UPLO, N, AFAC, RWORK )
144 *
145 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
146 *
147 RETURN
148 *
149 * End of SPPT01
150 *
151 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 N
10 REAL RESID
11 * ..
12 * .. Array Arguments ..
13 REAL A( * ), AFAC( * ), RWORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * SPPT01 reconstructs a symmetric positive definite packed matrix A
20 * from its L*L' or U'*U factorization and computes the residual
21 * norm( L*L' - A ) / ( N * norm(A) * EPS ) or
22 * norm( U'*U - A ) / ( N * norm(A) * EPS ),
23 * where 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/output) REAL array, dimension (N*(N+1)/2)
42 * On entry, the factor L or U from the L*L' or U'*U
43 * factorization of A, stored as a packed triangular matrix.
44 * Overwritten with the reconstructed matrix, and then with the
45 * difference L*L' - A (or U'*U - A).
46 *
47 * RWORK (workspace) REAL array, dimension (N)
48 *
49 * RESID (output) REAL
50 * If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
51 * If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
52 *
53 * =====================================================================
54 *
55 * .. Parameters ..
56 REAL ZERO, ONE
57 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
58 * ..
59 * .. Local Scalars ..
60 INTEGER I, K, KC, NPP
61 REAL ANORM, EPS, T
62 * ..
63 * .. External Functions ..
64 LOGICAL LSAME
65 REAL SDOT, SLAMCH, SLANSP
66 EXTERNAL LSAME, SDOT, SLAMCH, SLANSP
67 * ..
68 * .. External Subroutines ..
69 EXTERNAL SSCAL, SSPR, STPMV
70 * ..
71 * .. Intrinsic Functions ..
72 INTRINSIC REAL
73 * ..
74 * .. Executable Statements ..
75 *
76 * Quick exit if N = 0
77 *
78 IF( N.LE.0 ) THEN
79 RESID = ZERO
80 RETURN
81 END IF
82 *
83 * Exit with RESID = 1/EPS if ANORM = 0.
84 *
85 EPS = SLAMCH( 'Epsilon' )
86 ANORM = SLANSP( '1', UPLO, N, A, RWORK )
87 IF( ANORM.LE.ZERO ) THEN
88 RESID = ONE / EPS
89 RETURN
90 END IF
91 *
92 * Compute the product U'*U, overwriting U.
93 *
94 IF( LSAME( UPLO, 'U' ) ) THEN
95 KC = ( N*( N-1 ) ) / 2 + 1
96 DO 10 K = N, 1, -1
97 *
98 * Compute the (K,K) element of the result.
99 *
100 T = SDOT( K, AFAC( KC ), 1, AFAC( KC ), 1 )
101 AFAC( KC+K-1 ) = T
102 *
103 * Compute the rest of column K.
104 *
105 IF( K.GT.1 ) THEN
106 CALL STPMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
107 $ AFAC( KC ), 1 )
108 KC = KC - ( K-1 )
109 END IF
110 10 CONTINUE
111 *
112 * Compute the product L*L', overwriting L.
113 *
114 ELSE
115 KC = ( N*( N+1 ) ) / 2
116 DO 20 K = N, 1, -1
117 *
118 * Add a multiple of column K of the factor L to each of
119 * columns K+1 through N.
120 *
121 IF( K.LT.N )
122 $ CALL SSPR( 'Lower', N-K, ONE, AFAC( KC+1 ), 1,
123 $ AFAC( KC+N-K+1 ) )
124 *
125 * Scale column K by the diagonal element.
126 *
127 T = AFAC( KC )
128 CALL SSCAL( N-K+1, T, AFAC( KC ), 1 )
129 *
130 KC = KC - ( N-K+2 )
131 20 CONTINUE
132 END IF
133 *
134 * Compute the difference L*L' - A (or U'*U - A).
135 *
136 NPP = N*( N+1 ) / 2
137 DO 30 I = 1, NPP
138 AFAC( I ) = AFAC( I ) - A( I )
139 30 CONTINUE
140 *
141 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
142 *
143 RESID = SLANSP( '1', UPLO, N, AFAC, RWORK )
144 *
145 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
146 *
147 RETURN
148 *
149 * End of SPPT01
150 *
151 END