LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cpst01()

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.

Definition at line 134 of file cpst01.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  REAL RESID
143  INTEGER LDA, LDAFAC, LDPERM, N, RANK
144  CHARACTER UPLO
145 * ..
146 * .. Array Arguments ..
147  COMPLEX A( LDA, * ), AFAC( LDAFAC, * ),
148  $ PERM( LDPERM, * )
149  REAL RWORK( * )
150  INTEGER PIV( * )
151 * ..
152 *
153 * =====================================================================
154 *
155 * .. Parameters ..
156  REAL ZERO, ONE
157  parameter( zero = 0.0e+0, one = 1.0e+0 )
158  COMPLEX CZERO
159  parameter( czero = ( 0.0e+0, 0.0e+0 ) )
160 * ..
161 * .. Local Scalars ..
162  COMPLEX TC
163  REAL ANORM, EPS, TR
164  INTEGER I, J, K
165 * ..
166 * .. External Functions ..
167  COMPLEX CDOTC
168  REAL CLANHE, SLAMCH
169  LOGICAL LSAME
170  EXTERNAL cdotc, clanhe, slamch, lsame
171 * ..
172 * .. External Subroutines ..
173  EXTERNAL cher, cscal, ctrmv
174 * ..
175 * .. Intrinsic Functions ..
176  INTRINSIC aimag, conjg, real
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 = slamch( 'Epsilon' )
190  anorm = clanhe( '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( aimag( 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 = cdotc( 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 ctrmv( '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 cher( '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 cscal( 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 ) ) = conjg( 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 ) ) = conjg( 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 ) - real( a( j, j ) )
301  220 CONTINUE
302  ELSE
303  DO 240 j = 1, n
304  perm( j, j ) = perm( j, j ) - real( 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 = clanhe( '1', uplo, n, perm, ldafac, rwork )
315 *
316  resid = ( ( resid / real( n ) ) / anorm ) / eps
317 *
318  RETURN
319 *
320 * End of CPST01
321 *
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:83
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine ctrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
CTRMV
Definition: ctrmv.f:147
subroutine cher(UPLO, N, ALPHA, X, INCX, A, LDA)
CHER
Definition: cher.f:135
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,...
Definition: clanhe.f:124
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: