Files
lapack/TESTING/LIN/zppt02.f
Julie e18d437924 Making LAPACK's code eternal... no more version and date in source files.
GitHub is now enabling us to track accurately version and date.
No need for this anymore.
2021-03-25 10:16:58 -07:00

204 lines
5.3 KiB
FortranFixed

*> \brief \b ZPPT02
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE ZPPT02( UPLO, N, NRHS, A, X, LDX, B, LDB, RWORK,
* RESID )
*
* .. Scalar Arguments ..
* CHARACTER UPLO
* INTEGER LDB, LDX, N, NRHS
* DOUBLE PRECISION RESID
* ..
* .. Array Arguments ..
* DOUBLE PRECISION RWORK( * )
* COMPLEX*16 A( * ), B( LDB, * ), X( LDX, * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> ZPPT02 computes the residual in the solution of a Hermitian system
*> of linear equations A*x = b when packed storage is used for the
*> coefficient matrix. The ratio computed is
*>
*> RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS),
*>
*> where EPS is the machine precision.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> 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
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The number of rows and columns of the matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*> NRHS is INTEGER
*> The number of columns of B, the matrix of right hand sides.
*> NRHS >= 0.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is COMPLEX*16 array, dimension (N*(N+1)/2)
*> The original Hermitian matrix A, stored as a packed
*> triangular matrix.
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*> X is COMPLEX*16 array, dimension (LDX,NRHS)
*> The computed solution vectors for the system of linear
*> equations.
*> \endverbatim
*>
*> \param[in] LDX
*> \verbatim
*> LDX is INTEGER
*> The leading dimension of the array X. LDX >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*> B is COMPLEX*16 array, dimension (LDB,NRHS)
*> On entry, the right hand side vectors for the system of
*> linear equations.
*> On exit, B is overwritten with the difference B - A*X.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N)
*> \endverbatim
*>
*> \param[out] RESID
*> \verbatim
*> RESID is DOUBLE PRECISION
*> The maximum over the number of right hand sides of
*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complex16_lin
*
* =====================================================================
SUBROUTINE ZPPT02( UPLO, N, NRHS, A, X, LDX, B, LDB, RWORK,
$ RESID )
*
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
CHARACTER UPLO
INTEGER LDB, LDX, N, NRHS
DOUBLE PRECISION RESID
* ..
* .. Array Arguments ..
DOUBLE PRECISION RWORK( * )
COMPLEX*16 A( * ), B( LDB, * ), X( LDX, * )
* ..
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
COMPLEX*16 CONE
PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
* ..
* .. Local Scalars ..
INTEGER J
DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH, DZASUM, ZLANHP
EXTERNAL DLAMCH, DZASUM, ZLANHP
* ..
* .. External Subroutines ..
EXTERNAL ZHPMV
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* Quick exit if N = 0 or NRHS = 0.
*
IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
RESID = ZERO
RETURN
END IF
*
* Exit with RESID = 1/EPS if ANORM = 0.
*
EPS = DLAMCH( 'Epsilon' )
ANORM = ZLANHP( '1', UPLO, N, A, RWORK )
IF( ANORM.LE.ZERO ) THEN
RESID = ONE / EPS
RETURN
END IF
*
* Compute B - A*X for the matrix of right hand sides B.
*
DO 10 J = 1, NRHS
CALL ZHPMV( UPLO, N, -CONE, A, X( 1, J ), 1, CONE, B( 1, J ),
$ 1 )
10 CONTINUE
*
* Compute the maximum over the number of right hand sides of
* norm( B - A*X ) / ( norm(A) * norm(X) * EPS ) .
*
RESID = ZERO
DO 20 J = 1, NRHS
BNORM = DZASUM( N, B( 1, J ), 1 )
XNORM = DZASUM( N, X( 1, J ), 1 )
IF( XNORM.LE.ZERO ) THEN
RESID = ONE / EPS
ELSE
RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
END IF
20 CONTINUE
*
RETURN
*
* End of ZPPT02
*
END