LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 138 of file cpst01.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  REAL resid
146  INTEGER lda, ldafac, ldperm, n, rank
147  CHARACTER uplo
148 * ..
149 * .. Array Arguments ..
150  COMPLEX a( lda, * ), afac( ldafac, * ),
151  $ perm( ldperm, * )
152  REAL rwork( * )
153  INTEGER piv( * )
154 * ..
155 *
156 * =====================================================================
157 *
158 * .. Parameters ..
159  REAL zero, one
160  parameter ( zero = 0.0e+0, one = 1.0e+0 )
161  COMPLEX czero
162  parameter ( czero = ( 0.0e+0, 0.0e+0 ) )
163 * ..
164 * .. Local Scalars ..
165  COMPLEX tc
166  REAL anorm, eps, tr
167  INTEGER i, j, k
168 * ..
169 * .. External Functions ..
170  COMPLEX cdotc
171  REAL clanhe, slamch
172  LOGICAL lsame
173  EXTERNAL cdotc, clanhe, slamch, lsame
174 * ..
175 * .. External Subroutines ..
176  EXTERNAL cher, cscal, ctrmv
177 * ..
178 * .. Intrinsic Functions ..
179  INTRINSIC aimag, conjg, real
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 = slamch( 'Epsilon' )
193  anorm = clanhe( '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( aimag( 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 = cdotc( 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 ctrmv( '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 cher( '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 cscal( 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 ) ) = conjg( 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 ) ) = conjg( 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 ) - REAL( A( J, J ) )
304  220 CONTINUE
305  ELSE
306  DO 240 j = 1, n
307  perm( j, j ) = perm( j, j ) - REAL( 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 = clanhe( '1', uplo, n, perm, ldafac, rwork )
318 *
319  resid = ( ( resid / REAL( N ) ) / anorm ) / eps
320 *
321  RETURN
322 *
323 * End of CPST01
324 *
cher
subroutine cher(UPLO, N, ALPHA, X, INCX, A, LDA)
CHER
Definition: cher.f:137
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, or the element of largest absolute value of a complex Hermitian matrix.
Definition: clanhe.f:126
cscal
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
cdotc
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:54
ctrmv
subroutine ctrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
CTRMV
Definition: ctrmv.f:149
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: