1 SUBROUTINE SPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
2 $ PIV, RWORK, RESID, RANK )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Craig Lucas, University of Manchester / NAG Ltd.
6 * October, 2008
7 *
8 * .. Scalar Arguments ..
9 REAL RESID
10 INTEGER LDA, LDAFAC, LDPERM, N, RANK
11 CHARACTER UPLO
12 * ..
13 * .. Array Arguments ..
14 REAL A( LDA, * ), AFAC( LDAFAC, * ),
15 $ PERM( LDPERM, * ), RWORK( * )
16 INTEGER PIV( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * SPST01 reconstructs a symmetric positive semidefinite matrix A
23 * from its L or U factors and the permutation matrix P and computes
24 * the residual
25 * norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
26 * norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
27 * where EPS is the machine epsilon.
28 *
29 * Arguments
30 * ==========
31 *
32 * UPLO (input) CHARACTER*1
33 * Specifies whether the upper or lower triangular part of the
34 * symmetric matrix A is stored:
35 * = 'U': Upper triangular
36 * = 'L': Lower triangular
37 *
38 * N (input) INTEGER
39 * The number of rows and columns of the matrix A. N >= 0.
40 *
41 * A (input) REAL array, dimension (LDA,N)
42 * The original symmetric matrix A.
43 *
44 * LDA (input) INTEGER
45 * The leading dimension of the array A. LDA >= max(1,N)
46 *
47 * AFAC (input) REAL array, dimension (LDAFAC,N)
48 * The factor L or U from the L*L' or U'*U
49 * factorization of A.
50 *
51 * LDAFAC (input) INTEGER
52 * The leading dimension of the array AFAC. LDAFAC >= max(1,N).
53 *
54 * PERM (output) REAL array, dimension (LDPERM,N)
55 * Overwritten with the reconstructed matrix, and then with the
56 * difference P*L*L'*P' - A (or P*U'*U*P' - A)
57 *
58 * LDPERM (input) INTEGER
59 * The leading dimension of the array PERM.
60 * LDAPERM >= max(1,N).
61 *
62 * PIV (input) INTEGER array, dimension (N)
63 * PIV is such that the nonzero entries are
64 * P( PIV( K ), K ) = 1.
65 *
66 * RWORK (workspace) REAL array, dimension (N)
67 *
68 * RESID (output) REAL
69 * If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
70 * If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 REAL ZERO, ONE
76 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
77 * ..
78 * .. Local Scalars ..
79 REAL ANORM, EPS, T
80 INTEGER I, J, K
81 * ..
82 * .. External Functions ..
83 REAL SDOT, SLAMCH, SLANSY
84 LOGICAL LSAME
85 EXTERNAL SDOT, SLAMCH, SLANSY, LSAME
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL SSCAL, SSYR, STRMV
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC REAL
92 * ..
93 * .. Executable Statements ..
94 *
95 * Quick exit if N = 0.
96 *
97 IF( N.LE.0 ) THEN
98 RESID = ZERO
99 RETURN
100 END IF
101 *
102 * Exit with RESID = 1/EPS if ANORM = 0.
103 *
104 EPS = SLAMCH( 'Epsilon' )
105 ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK )
106 IF( ANORM.LE.ZERO ) THEN
107 RESID = ONE / EPS
108 RETURN
109 END IF
110 *
111 * Compute the product U'*U, overwriting U.
112 *
113 IF( LSAME( UPLO, 'U' ) ) THEN
114 *
115 IF( RANK.LT.N ) THEN
116 DO 110 J = RANK + 1, N
117 DO 100 I = RANK + 1, J
118 AFAC( I, J ) = ZERO
119 100 CONTINUE
120 110 CONTINUE
121 END IF
122 *
123 DO 120 K = N, 1, -1
124 *
125 * Compute the (K,K) element of the result.
126 *
127 T = SDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
128 AFAC( K, K ) = T
129 *
130 * Compute the rest of column K.
131 *
132 CALL STRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
133 $ LDAFAC, AFAC( 1, K ), 1 )
134 *
135 120 CONTINUE
136 *
137 * Compute the product L*L', overwriting L.
138 *
139 ELSE
140 *
141 IF( RANK.LT.N ) THEN
142 DO 140 J = RANK + 1, N
143 DO 130 I = J, N
144 AFAC( I, J ) = ZERO
145 130 CONTINUE
146 140 CONTINUE
147 END IF
148 *
149 DO 150 K = N, 1, -1
150 * Add a multiple of column K of the factor L to each of
151 * columns K+1 through N.
152 *
153 IF( K+1.LE.N )
154 $ CALL SSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
155 $ AFAC( K+1, K+1 ), LDAFAC )
156 *
157 * Scale column K by the diagonal element.
158 *
159 T = AFAC( K, K )
160 CALL SSCAL( N-K+1, T, AFAC( K, K ), 1 )
161 150 CONTINUE
162 *
163 END IF
164 *
165 * Form P*L*L'*P' or P*U'*U*P'
166 *
167 IF( LSAME( UPLO, 'U' ) ) THEN
168 *
169 DO 170 J = 1, N
170 DO 160 I = 1, N
171 IF( PIV( I ).LE.PIV( J ) ) THEN
172 IF( I.LE.J ) THEN
173 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
174 ELSE
175 PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
176 END IF
177 END IF
178 160 CONTINUE
179 170 CONTINUE
180 *
181 *
182 ELSE
183 *
184 DO 190 J = 1, N
185 DO 180 I = 1, N
186 IF( PIV( I ).GE.PIV( J ) ) THEN
187 IF( I.GE.J ) THEN
188 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
189 ELSE
190 PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
191 END IF
192 END IF
193 180 CONTINUE
194 190 CONTINUE
195 *
196 END IF
197 *
198 * Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A).
199 *
200 IF( LSAME( UPLO, 'U' ) ) THEN
201 DO 210 J = 1, N
202 DO 200 I = 1, J
203 PERM( I, J ) = PERM( I, J ) - A( I, J )
204 200 CONTINUE
205 210 CONTINUE
206 ELSE
207 DO 230 J = 1, N
208 DO 220 I = J, N
209 PERM( I, J ) = PERM( I, J ) - A( I, J )
210 220 CONTINUE
211 230 CONTINUE
212 END IF
213 *
214 * Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
215 * ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
216 *
217 RESID = SLANSY( '1', UPLO, N, PERM, LDAFAC, RWORK )
218 *
219 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
220 *
221 RETURN
222 *
223 * End of SPST01
224 *
225 END
2 $ PIV, RWORK, RESID, RANK )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Craig Lucas, University of Manchester / NAG Ltd.
6 * October, 2008
7 *
8 * .. Scalar Arguments ..
9 REAL RESID
10 INTEGER LDA, LDAFAC, LDPERM, N, RANK
11 CHARACTER UPLO
12 * ..
13 * .. Array Arguments ..
14 REAL A( LDA, * ), AFAC( LDAFAC, * ),
15 $ PERM( LDPERM, * ), RWORK( * )
16 INTEGER PIV( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * SPST01 reconstructs a symmetric positive semidefinite matrix A
23 * from its L or U factors and the permutation matrix P and computes
24 * the residual
25 * norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
26 * norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
27 * where EPS is the machine epsilon.
28 *
29 * Arguments
30 * ==========
31 *
32 * UPLO (input) CHARACTER*1
33 * Specifies whether the upper or lower triangular part of the
34 * symmetric matrix A is stored:
35 * = 'U': Upper triangular
36 * = 'L': Lower triangular
37 *
38 * N (input) INTEGER
39 * The number of rows and columns of the matrix A. N >= 0.
40 *
41 * A (input) REAL array, dimension (LDA,N)
42 * The original symmetric matrix A.
43 *
44 * LDA (input) INTEGER
45 * The leading dimension of the array A. LDA >= max(1,N)
46 *
47 * AFAC (input) REAL array, dimension (LDAFAC,N)
48 * The factor L or U from the L*L' or U'*U
49 * factorization of A.
50 *
51 * LDAFAC (input) INTEGER
52 * The leading dimension of the array AFAC. LDAFAC >= max(1,N).
53 *
54 * PERM (output) REAL array, dimension (LDPERM,N)
55 * Overwritten with the reconstructed matrix, and then with the
56 * difference P*L*L'*P' - A (or P*U'*U*P' - A)
57 *
58 * LDPERM (input) INTEGER
59 * The leading dimension of the array PERM.
60 * LDAPERM >= max(1,N).
61 *
62 * PIV (input) INTEGER array, dimension (N)
63 * PIV is such that the nonzero entries are
64 * P( PIV( K ), K ) = 1.
65 *
66 * RWORK (workspace) REAL array, dimension (N)
67 *
68 * RESID (output) REAL
69 * If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
70 * If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 REAL ZERO, ONE
76 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
77 * ..
78 * .. Local Scalars ..
79 REAL ANORM, EPS, T
80 INTEGER I, J, K
81 * ..
82 * .. External Functions ..
83 REAL SDOT, SLAMCH, SLANSY
84 LOGICAL LSAME
85 EXTERNAL SDOT, SLAMCH, SLANSY, LSAME
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL SSCAL, SSYR, STRMV
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC REAL
92 * ..
93 * .. Executable Statements ..
94 *
95 * Quick exit if N = 0.
96 *
97 IF( N.LE.0 ) THEN
98 RESID = ZERO
99 RETURN
100 END IF
101 *
102 * Exit with RESID = 1/EPS if ANORM = 0.
103 *
104 EPS = SLAMCH( 'Epsilon' )
105 ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK )
106 IF( ANORM.LE.ZERO ) THEN
107 RESID = ONE / EPS
108 RETURN
109 END IF
110 *
111 * Compute the product U'*U, overwriting U.
112 *
113 IF( LSAME( UPLO, 'U' ) ) THEN
114 *
115 IF( RANK.LT.N ) THEN
116 DO 110 J = RANK + 1, N
117 DO 100 I = RANK + 1, J
118 AFAC( I, J ) = ZERO
119 100 CONTINUE
120 110 CONTINUE
121 END IF
122 *
123 DO 120 K = N, 1, -1
124 *
125 * Compute the (K,K) element of the result.
126 *
127 T = SDOT( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
128 AFAC( K, K ) = T
129 *
130 * Compute the rest of column K.
131 *
132 CALL STRMV( 'Upper', 'Transpose', 'Non-unit', K-1, AFAC,
133 $ LDAFAC, AFAC( 1, K ), 1 )
134 *
135 120 CONTINUE
136 *
137 * Compute the product L*L', overwriting L.
138 *
139 ELSE
140 *
141 IF( RANK.LT.N ) THEN
142 DO 140 J = RANK + 1, N
143 DO 130 I = J, N
144 AFAC( I, J ) = ZERO
145 130 CONTINUE
146 140 CONTINUE
147 END IF
148 *
149 DO 150 K = N, 1, -1
150 * Add a multiple of column K of the factor L to each of
151 * columns K+1 through N.
152 *
153 IF( K+1.LE.N )
154 $ CALL SSYR( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
155 $ AFAC( K+1, K+1 ), LDAFAC )
156 *
157 * Scale column K by the diagonal element.
158 *
159 T = AFAC( K, K )
160 CALL SSCAL( N-K+1, T, AFAC( K, K ), 1 )
161 150 CONTINUE
162 *
163 END IF
164 *
165 * Form P*L*L'*P' or P*U'*U*P'
166 *
167 IF( LSAME( UPLO, 'U' ) ) THEN
168 *
169 DO 170 J = 1, N
170 DO 160 I = 1, N
171 IF( PIV( I ).LE.PIV( J ) ) THEN
172 IF( I.LE.J ) THEN
173 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
174 ELSE
175 PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
176 END IF
177 END IF
178 160 CONTINUE
179 170 CONTINUE
180 *
181 *
182 ELSE
183 *
184 DO 190 J = 1, N
185 DO 180 I = 1, N
186 IF( PIV( I ).GE.PIV( J ) ) THEN
187 IF( I.GE.J ) THEN
188 PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
189 ELSE
190 PERM( PIV( I ), PIV( J ) ) = AFAC( J, I )
191 END IF
192 END IF
193 180 CONTINUE
194 190 CONTINUE
195 *
196 END IF
197 *
198 * Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A).
199 *
200 IF( LSAME( UPLO, 'U' ) ) THEN
201 DO 210 J = 1, N
202 DO 200 I = 1, J
203 PERM( I, J ) = PERM( I, J ) - A( I, J )
204 200 CONTINUE
205 210 CONTINUE
206 ELSE
207 DO 230 J = 1, N
208 DO 220 I = J, N
209 PERM( I, J ) = PERM( I, J ) - A( I, J )
210 220 CONTINUE
211 230 CONTINUE
212 END IF
213 *
214 * Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
215 * ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
216 *
217 RESID = SLANSY( '1', UPLO, N, PERM, LDAFAC, RWORK )
218 *
219 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
220 *
221 RETURN
222 *
223 * End of SPST01
224 *
225 END