LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cpst01.f
Go to the documentation of this file.
1*> \brief \b CPST01
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE CPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
12* PIV, RWORK, RESID, RANK )
13*
14* .. Scalar Arguments ..
15* REAL RESID
16* INTEGER LDA, LDAFAC, LDPERM, N, RANK
17* CHARACTER UPLO
18* ..
19* .. Array Arguments ..
20* COMPLEX A( LDA, * ), AFAC( LDAFAC, * ),
21* $ PERM( LDPERM, * )
22* REAL RWORK( * )
23* INTEGER PIV( * )
24* ..
25*
26*
27*> \par Purpose:
28* =============
29*>
30*> \verbatim
31*>
32*> CPST01 reconstructs an Hermitian positive semidefinite matrix A
33*> from its L or U factors and the permutation matrix P and computes
34*> the residual
35*> norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
36*> norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
37*> where EPS is the machine epsilon, L' is the conjugate transpose of L,
38*> and U' is the conjugate transpose of U.
39*> \endverbatim
40*
41* Arguments:
42* ==========
43*
44*> \param[in] UPLO
45*> \verbatim
46*> UPLO is CHARACTER*1
47*> Specifies whether the upper or lower triangular part of the
48*> Hermitian matrix A is stored:
49*> = 'U': Upper triangular
50*> = 'L': Lower triangular
51*> \endverbatim
52*>
53*> \param[in] N
54*> \verbatim
55*> N is INTEGER
56*> The number of rows and columns of the matrix A. N >= 0.
57*> \endverbatim
58*>
59*> \param[in] A
60*> \verbatim
61*> A is COMPLEX array, dimension (LDA,N)
62*> The original Hermitian matrix A.
63*> \endverbatim
64*>
65*> \param[in] LDA
66*> \verbatim
67*> LDA is INTEGER
68*> The leading dimension of the array A. LDA >= max(1,N)
69*> \endverbatim
70*>
71*> \param[in] AFAC
72*> \verbatim
73*> AFAC is COMPLEX array, dimension (LDAFAC,N)
74*> The factor L or U from the L*L' or U'*U
75*> factorization of A.
76*> \endverbatim
77*>
78*> \param[in] LDAFAC
79*> \verbatim
80*> LDAFAC is INTEGER
81*> The leading dimension of the array AFAC. LDAFAC >= max(1,N).
82*> \endverbatim
83*>
84*> \param[out] PERM
85*> \verbatim
86*> PERM is COMPLEX array, dimension (LDPERM,N)
87*> Overwritten with the reconstructed matrix, and then with the
88*> difference P*L*L'*P' - A (or P*U'*U*P' - A)
89*> \endverbatim
90*>
91*> \param[in] LDPERM
92*> \verbatim
93*> LDPERM is INTEGER
94*> The leading dimension of the array PERM.
95*> LDAPERM >= max(1,N).
96*> \endverbatim
97*>
98*> \param[in] PIV
99*> \verbatim
100*> PIV is INTEGER array, dimension (N)
101*> PIV is such that the nonzero entries are
102*> P( PIV( K ), K ) = 1.
103*> \endverbatim
104*>
105*> \param[out] RWORK
106*> \verbatim
107*> RWORK is REAL array, dimension (N)
108*> \endverbatim
109*>
110*> \param[out] RESID
111*> \verbatim
112*> RESID is REAL
113*> If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
114*> If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
115*> \endverbatim
116*>
117*> \param[in] RANK
118*> \verbatim
119*> RANK is INTEGER
120*> number of nonzero singular values of A.
121*> \endverbatim
122*
123* Authors:
124* ========
125*
126*> \author Univ. of Tennessee
127*> \author Univ. of California Berkeley
128*> \author Univ. of Colorado Denver
129*> \author NAG Ltd.
130*
131*> \ingroup complex_lin
132*
133* =====================================================================
134 SUBROUTINE cpst01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
135 $ PIV, RWORK, RESID, RANK )
136*
137* -- LAPACK test routine --
138* -- LAPACK is a software package provided by Univ. of Tennessee, --
139* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140*
141* .. Scalar Arguments ..
142 REAL RESID
143 INTEGER LDA, LDAFAC, LDPERM, N, RANK
144 CHARACTER UPLO
145* ..
146* .. Array Arguments ..
147 COMPLEX A( LDA, * ), AFAC( LDAFAC, * ),
148 $ perm( ldperm, * )
149 REAL RWORK( * )
150 INTEGER PIV( * )
151* ..
152*
153* =====================================================================
154*
155* .. Parameters ..
156 REAL ZERO, ONE
157 parameter( zero = 0.0e+0, one = 1.0e+0 )
158 COMPLEX CZERO
159 parameter( czero = ( 0.0e+0, 0.0e+0 ) )
160* ..
161* .. Local Scalars ..
162 COMPLEX TC
163 REAL ANORM, EPS, TR
164 INTEGER I, J, K
165* ..
166* .. External Functions ..
167 COMPLEX CDOTC
168 REAL CLANHE, SLAMCH
169 LOGICAL LSAME
170 EXTERNAL cdotc, clanhe, slamch, lsame
171* ..
172* .. External Subroutines ..
173 EXTERNAL cher, cscal, ctrmv
174* ..
175* .. Intrinsic Functions ..
176 INTRINSIC aimag, conjg, real
177* ..
178* .. Executable Statements ..
179*
180* Quick exit if N = 0.
181*
182 IF( n.LE.0 ) THEN
183 resid = zero
184 RETURN
185 END IF
186*
187* Exit with RESID = 1/EPS if ANORM = 0.
188*
189 eps = slamch( 'Epsilon' )
190 anorm = clanhe( '1', uplo, n, a, lda, rwork )
191 IF( anorm.LE.zero ) THEN
192 resid = one / eps
193 RETURN
194 END IF
195*
196* Check the imaginary parts of the diagonal elements and return with
197* an error code if any are nonzero.
198*
199 DO 100 j = 1, n
200 IF( aimag( afac( j, j ) ).NE.zero ) THEN
201 resid = one / eps
202 RETURN
203 END IF
204 100 CONTINUE
205*
206* Compute the product U'*U, overwriting U.
207*
208 IF( lsame( uplo, 'U' ) ) THEN
209*
210 IF( rank.LT.n ) THEN
211 DO 120 j = rank + 1, n
212 DO 110 i = rank + 1, j
213 afac( i, j ) = czero
214 110 CONTINUE
215 120 CONTINUE
216 END IF
217*
218 DO 130 k = n, 1, -1
219*
220* Compute the (K,K) element of the result.
221*
222 tr = real( cdotc( k, afac( 1, k ), 1, afac( 1, k ), 1 ) )
223 afac( k, k ) = tr
224*
225* Compute the rest of column K.
226*
227 CALL ctrmv( 'Upper', 'Conjugate', 'Non-unit', k-1, afac,
228 $ ldafac, afac( 1, k ), 1 )
229*
230 130 CONTINUE
231*
232* Compute the product L*L', overwriting L.
233*
234 ELSE
235*
236 IF( rank.LT.n ) THEN
237 DO 150 j = rank + 1, n
238 DO 140 i = j, n
239 afac( i, j ) = czero
240 140 CONTINUE
241 150 CONTINUE
242 END IF
243*
244 DO 160 k = n, 1, -1
245* Add a multiple of column K of the factor L to each of
246* columns K+1 through N.
247*
248 IF( k+1.LE.n )
249 $ CALL cher( 'Lower', n-k, one, afac( k+1, k ), 1,
250 $ afac( k+1, k+1 ), ldafac )
251*
252* Scale column K by the diagonal element.
253*
254 tc = afac( k, k )
255 CALL cscal( n-k+1, tc, afac( k, k ), 1 )
256 160 CONTINUE
257*
258 END IF
259*
260* Form P*L*L'*P' or P*U'*U*P'
261*
262 IF( lsame( uplo, 'U' ) ) THEN
263*
264 DO 180 j = 1, n
265 DO 170 i = 1, n
266 IF( piv( i ).LE.piv( j ) ) THEN
267 IF( i.LE.j ) THEN
268 perm( piv( i ), piv( j ) ) = afac( i, j )
269 ELSE
270 perm( piv( i ), piv( j ) ) = conjg( afac( j, i ) )
271 END IF
272 END IF
273 170 CONTINUE
274 180 CONTINUE
275*
276*
277 ELSE
278*
279 DO 200 j = 1, n
280 DO 190 i = 1, n
281 IF( piv( i ).GE.piv( j ) ) THEN
282 IF( i.GE.j ) THEN
283 perm( piv( i ), piv( j ) ) = afac( i, j )
284 ELSE
285 perm( piv( i ), piv( j ) ) = conjg( afac( j, i ) )
286 END IF
287 END IF
288 190 CONTINUE
289 200 CONTINUE
290*
291 END IF
292*
293* Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A).
294*
295 IF( lsame( uplo, 'U' ) ) THEN
296 DO 220 j = 1, n
297 DO 210 i = 1, j - 1
298 perm( i, j ) = perm( i, j ) - a( i, j )
299 210 CONTINUE
300 perm( j, j ) = perm( j, j ) - real( a( j, j ) )
301 220 CONTINUE
302 ELSE
303 DO 240 j = 1, n
304 perm( j, j ) = perm( j, j ) - real( a( j, j ) )
305 DO 230 i = j + 1, n
306 perm( i, j ) = perm( i, j ) - a( i, j )
307 230 CONTINUE
308 240 CONTINUE
309 END IF
310*
311* Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
312* ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
313*
314 resid = clanhe( '1', uplo, n, perm, ldafac, rwork )
315*
316 resid = ( ( resid / real( n ) ) / anorm ) / eps
317*
318 RETURN
319*
320* End of CPST01
321*
322 END
subroutine cpst01(uplo, n, a, lda, afac, ldafac, perm, ldperm, piv, rwork, resid, rank)
CPST01
Definition cpst01.f:136
subroutine cher(uplo, n, alpha, x, incx, a, lda)
CHER
Definition cher.f:135
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
subroutine ctrmv(uplo, trans, diag, n, a, lda, x, incx)
CTRMV
Definition ctrmv.f:147