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