 LAPACK  3.8.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.```
Date
December 2016

Definition at line 138 of file zpst01.f.

138 *
139 * -- LAPACK test routine (version 3.7.0) --
140 * -- LAPACK is a software package provided by Univ. of Tennessee, --
141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142 * December 2016
143 *
144 * .. Scalar Arguments ..
145  DOUBLE PRECISION resid
146  INTEGER lda, ldafac, ldperm, n, rank
147  CHARACTER uplo
148 * ..
149 * .. Array Arguments ..
150  COMPLEX*16 a( lda, * ), afac( ldafac, * ),
151  \$ perm( ldperm, * )
152  DOUBLE PRECISION rwork( * )
153  INTEGER piv( * )
154 * ..
155 *
156 * =====================================================================
157 *
158 * .. Parameters ..
159  DOUBLE PRECISION zero, one
160  parameter( zero = 0.0d+0, one = 1.0d+0 )
161  COMPLEX*16 czero
162  parameter( czero = ( 0.0d+0, 0.0d+0 ) )
163 * ..
164 * .. Local Scalars ..
165  COMPLEX*16 tc
166  DOUBLE PRECISION anorm, eps, tr
167  INTEGER i, j, k
168 * ..
169 * .. External Functions ..
170  COMPLEX*16 zdotc
171  DOUBLE PRECISION dlamch, zlanhe
172  LOGICAL lsame
173  EXTERNAL zdotc, dlamch, zlanhe, lsame
174 * ..
175 * .. External Subroutines ..
176  EXTERNAL zher, zscal, ztrmv
177 * ..
178 * .. Intrinsic Functions ..
179  INTRINSIC dble, dconjg, dimag
180 * ..
181 * .. Executable Statements ..
182 *
183 * Quick exit if N = 0.
184 *
185  IF( n.LE.0 ) THEN
186  resid = zero
187  RETURN
188  END IF
189 *
190 * Exit with RESID = 1/EPS if ANORM = 0.
191 *
192  eps = dlamch( 'Epsilon' )
193  anorm = zlanhe( '1', uplo, n, a, lda, rwork )
194  IF( anorm.LE.zero ) THEN
195  resid = one / eps
196  RETURN
197  END IF
198 *
199 * Check the imaginary parts of the diagonal elements and return with
200 * an error code if any are nonzero.
201 *
202  DO 100 j = 1, n
203  IF( dimag( afac( j, j ) ).NE.zero ) THEN
204  resid = one / eps
205  RETURN
206  END IF
207  100 CONTINUE
208 *
209 * Compute the product U'*U, overwriting U.
210 *
211  IF( lsame( uplo, 'U' ) ) THEN
212 *
213  IF( rank.LT.n ) THEN
214  DO 120 j = rank + 1, n
215  DO 110 i = rank + 1, j
216  afac( i, j ) = czero
217  110 CONTINUE
218  120 CONTINUE
219  END IF
220 *
221  DO 130 k = n, 1, -1
222 *
223 * Compute the (K,K) element of the result.
224 *
225  tr = zdotc( k, afac( 1, k ), 1, afac( 1, k ), 1 )
226  afac( k, k ) = tr
227 *
228 * Compute the rest of column K.
229 *
230  CALL ztrmv( 'Upper', 'Conjugate', 'Non-unit', k-1, afac,
231  \$ ldafac, afac( 1, k ), 1 )
232 *
233  130 CONTINUE
234 *
235 * Compute the product L*L', overwriting L.
236 *
237  ELSE
238 *
239  IF( rank.LT.n ) THEN
240  DO 150 j = rank + 1, n
241  DO 140 i = j, n
242  afac( i, j ) = czero
243  140 CONTINUE
244  150 CONTINUE
245  END IF
246 *
247  DO 160 k = n, 1, -1
248 * Add a multiple of column K of the factor L to each of
249 * columns K+1 through N.
250 *
251  IF( k+1.LE.n )
252  \$ CALL zher( 'Lower', n-k, one, afac( k+1, k ), 1,
253  \$ afac( k+1, k+1 ), ldafac )
254 *
255 * Scale column K by the diagonal element.
256 *
257  tc = afac( k, k )
258  CALL zscal( n-k+1, tc, afac( k, k ), 1 )
259  160 CONTINUE
260 *
261  END IF
262 *
263 * Form P*L*L'*P' or P*U'*U*P'
264 *
265  IF( lsame( uplo, 'U' ) ) THEN
266 *
267  DO 180 j = 1, n
268  DO 170 i = 1, n
269  IF( piv( i ).LE.piv( j ) ) THEN
270  IF( i.LE.j ) THEN
271  perm( piv( i ), piv( j ) ) = afac( i, j )
272  ELSE
273  perm( piv( i ), piv( j ) ) = dconjg( afac( j, i ) )
274  END IF
275  END IF
276  170 CONTINUE
277  180 CONTINUE
278 *
279 *
280  ELSE
281 *
282  DO 200 j = 1, n
283  DO 190 i = 1, n
284  IF( piv( i ).GE.piv( j ) ) THEN
285  IF( i.GE.j ) THEN
286  perm( piv( i ), piv( j ) ) = afac( i, j )
287  ELSE
288  perm( piv( i ), piv( j ) ) = dconjg( afac( j, i ) )
289  END IF
290  END IF
291  190 CONTINUE
292  200 CONTINUE
293 *
294  END IF
295 *
296 * Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A).
297 *
298  IF( lsame( uplo, 'U' ) ) THEN
299  DO 220 j = 1, n
300  DO 210 i = 1, j - 1
301  perm( i, j ) = perm( i, j ) - a( i, j )
302  210 CONTINUE
303  perm( j, j ) = perm( j, j ) - dble( a( j, j ) )
304  220 CONTINUE
305  ELSE
306  DO 240 j = 1, n
307  perm( j, j ) = perm( j, j ) - dble( a( j, j ) )
308  DO 230 i = j + 1, n
309  perm( i, j ) = perm( i, j ) - a( i, j )
310  230 CONTINUE
311  240 CONTINUE
312  END IF
313 *
314 * Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
315 * ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
316 *
317  resid = zlanhe( '1', uplo, n, perm, ldafac, rwork )
318 *
319  resid = ( ( resid / dble( n ) ) / anorm ) / eps
320 *
321  RETURN
322 *
323 * End of ZPST01
324 *
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, or the element of largest absolute value of a complex Hermitian matrix.
Definition: zlanhe.f:126
subroutine zher(UPLO, N, ALPHA, X, INCX, A, LDA)
ZHER
Definition: zher.f:137
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine ztrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRMV
Definition: ztrmv.f:149
complex *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
Definition: zdotc.f:85
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
Here is the call graph for this function:
Here is the caller graph for this function: