 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ zpst01()

 subroutine zpst01 ( character UPLO, integer N, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldafac, * ) AFAC, integer LDAFAC, complex*16, dimension( ldperm, * ) PERM, integer LDPERM, integer, dimension( * ) PIV, double precision, dimension( * ) RWORK, double precision RESID, integer RANK )

ZPST01

Purpose:
ZPST01 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*16 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*16 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*16 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 DOUBLE PRECISION array, dimension (N) [out] RESID RESID is DOUBLE PRECISION 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.

Definition at line 134 of file zpst01.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  DOUBLE PRECISION RESID
143  INTEGER LDA, LDAFAC, LDPERM, N, RANK
144  CHARACTER UPLO
145 * ..
146 * .. Array Arguments ..
147  COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ),
148  \$ PERM( LDPERM, * )
149  DOUBLE PRECISION RWORK( * )
150  INTEGER PIV( * )
151 * ..
152 *
153 * =====================================================================
154 *
155 * .. Parameters ..
156  DOUBLE PRECISION ZERO, ONE
157  parameter( zero = 0.0d+0, one = 1.0d+0 )
158  COMPLEX*16 CZERO
159  parameter( czero = ( 0.0d+0, 0.0d+0 ) )
160 * ..
161 * .. Local Scalars ..
162  COMPLEX*16 TC
163  DOUBLE PRECISION ANORM, EPS, TR
164  INTEGER I, J, K
165 * ..
166 * .. External Functions ..
167  COMPLEX*16 ZDOTC
168  DOUBLE PRECISION DLAMCH, ZLANHE
169  LOGICAL LSAME
170  EXTERNAL zdotc, dlamch, zlanhe, lsame
171 * ..
172 * .. External Subroutines ..
173  EXTERNAL zher, zscal, ztrmv
174 * ..
175 * .. Intrinsic Functions ..
176  INTRINSIC dble, dconjg, dimag
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 = dlamch( 'Epsilon' )
190  anorm = zlanhe( '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( dimag( 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 = zdotc( 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 ztrmv( '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 zher( '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 zscal( 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 ) ) = dconjg( 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 ) ) = dconjg( 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 ) - dble( a( j, j ) )
301  220 CONTINUE
302  ELSE
303  DO 240 j = 1, n
304  perm( j, j ) = perm( j, j ) - dble( 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 = zlanhe( '1', uplo, n, perm, ldafac, rwork )
315 *
316  resid = ( ( resid / dble( n ) ) / anorm ) / eps
317 *
318  RETURN
319 *
320 * End of ZPST01
321 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
complex *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
Definition: zdotc.f:83
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine ztrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRMV
Definition: ztrmv.f:147
subroutine zher(UPLO, N, ALPHA, X, INCX, A, LDA)
ZHER
Definition: zher.f:135
double precision function zlanhe(NORM, UPLO, N, A, LDA, WORK)
ZLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: zlanhe.f:124
Here is the call graph for this function:
Here is the caller graph for this function: