LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ cpst01()

subroutine cpst01 ( character  uplo,
integer  n,
complex, dimension( lda, * )  a,
integer  lda,
complex, dimension( ldafac, * )  afac,
integer  ldafac,
complex, dimension( ldperm, * )  perm,
integer  ldperm,
integer, dimension( * )  piv,
real, dimension( * )  rwork,
real  resid,
integer  rank 
)

CPST01

Purpose:
 CPST01 reconstructs an Hermitian positive semidefinite matrix A
 from its L or U factors and the permutation matrix P and computes
 the residual
    norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
    norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
 where EPS is the machine epsilon, L' is the conjugate transpose of L,
 and U' is the conjugate transpose of U.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          Hermitian matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]N
          N is INTEGER
          The number of rows and columns of the matrix A.  N >= 0.
[in]A
          A is COMPLEX array, dimension (LDA,N)
          The original Hermitian matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N)
[in]AFAC
          AFAC is COMPLEX array, dimension (LDAFAC,N)
          The factor L or U from the L*L' or U'*U
          factorization of A.
[in]LDAFAC
          LDAFAC is INTEGER
          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
[out]PERM
          PERM is COMPLEX array, dimension (LDPERM,N)
          Overwritten with the reconstructed matrix, and then with the
          difference P*L*L'*P' - A (or P*U'*U*P' - A)
[in]LDPERM
          LDPERM is INTEGER
          The leading dimension of the array PERM.
          LDAPERM >= max(1,N).
[in]PIV
          PIV is INTEGER array, dimension (N)
          PIV is such that the nonzero entries are
          P( PIV( K ), K ) = 1.
[out]RWORK
          RWORK is REAL array, dimension (N)
[out]RESID
          RESID is REAL
          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
[in]RANK
          RANK is INTEGER
          number of nonzero singular values of A.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 134 of file cpst01.f.

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*
cdotc
complex function cdotc(n, cx, incx, cy, incy)
CDOTC
Definition cdotc.f:83
cher
subroutine cher(uplo, n, alpha, x, incx, a, lda)
CHER
Definition cher.f:135
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
clanhe
real function clanhe(norm, uplo, n, a, lda, work)
CLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clanhe.f:124
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
ctrmv
subroutine ctrmv(uplo, trans, diag, n, a, lda, x, incx)
CTRMV
Definition ctrmv.f:147
Here is the call graph for this function:
Here is the caller graph for this function: