LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 138 of file zpst01.f.

138 *
139 * -- LAPACK test routine (version 3.4.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 * November 2011
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 *
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
zlanhe
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
zdotc
complex *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
Definition: zdotc.f:54
zher
subroutine zher(UPLO, N, ALPHA, X, INCX, A, LDA)
ZHER
Definition: zher.f:137
ztrmv
subroutine ztrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRMV
Definition: ztrmv.f:149
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: