LAPACK  3.10.1
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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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: