![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CTPT02 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE CTPT02( UPLO, TRANS, DIAG, N, NRHS, AP, X, LDX, B, LDB, 00012 * WORK, RWORK, RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER DIAG, TRANS, UPLO 00016 * INTEGER LDB, LDX, N, NRHS 00017 * REAL RESID 00018 * .. 00019 * .. Array Arguments .. 00020 * REAL RWORK( * ) 00021 * COMPLEX AP( * ), B( LDB, * ), WORK( * ), X( LDX, * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> CTPT02 computes the residual for the computed solution to a 00031 *> triangular system of linear equations A*x = b, A**T *x = b, or 00032 *> A**H *x = b, when the triangular matrix A is stored in packed format. 00033 *> Here A**T denotes the transpose of A, A**H denotes the conjugate 00034 *> transpose of A, and x and b are N by NRHS matrices. The test ratio 00035 *> is the maximum over the number of right hand sides of 00036 *> the maximum over the number of right hand sides of 00037 *> norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ), 00038 *> where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon. 00039 *> \endverbatim 00040 * 00041 * Arguments: 00042 * ========== 00043 * 00044 *> \param[in] UPLO 00045 *> \verbatim 00046 *> UPLO is CHARACTER*1 00047 *> Specifies whether the matrix A is upper or lower triangular. 00048 *> = 'U': Upper triangular 00049 *> = 'L': Lower triangular 00050 *> \endverbatim 00051 *> 00052 *> \param[in] TRANS 00053 *> \verbatim 00054 *> TRANS is CHARACTER*1 00055 *> Specifies the operation applied to A. 00056 *> = 'N': A *x = b (No transpose) 00057 *> = 'T': A**T *x = b (Transpose) 00058 *> = 'C': A**H *x = b (Conjugate transpose) 00059 *> \endverbatim 00060 *> 00061 *> \param[in] DIAG 00062 *> \verbatim 00063 *> DIAG is CHARACTER*1 00064 *> Specifies whether or not the matrix A is unit triangular. 00065 *> = 'N': Non-unit triangular 00066 *> = 'U': Unit triangular 00067 *> \endverbatim 00068 *> 00069 *> \param[in] N 00070 *> \verbatim 00071 *> N is INTEGER 00072 *> The order of the matrix A. N >= 0. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] NRHS 00076 *> \verbatim 00077 *> NRHS is INTEGER 00078 *> The number of right hand sides, i.e., the number of columns 00079 *> of the matrices X and B. NRHS >= 0. 00080 *> \endverbatim 00081 *> 00082 *> \param[in] AP 00083 *> \verbatim 00084 *> AP is COMPLEX array, dimension (N*(N+1)/2) 00085 *> The upper or lower triangular matrix A, packed columnwise in 00086 *> a linear array. The j-th column of A is stored in the array 00087 *> AP as follows: 00088 *> if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j; 00089 *> if UPLO = 'L', 00090 *> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n. 00091 *> \endverbatim 00092 *> 00093 *> \param[in] X 00094 *> \verbatim 00095 *> X is COMPLEX array, dimension (LDX,NRHS) 00096 *> The computed solution vectors for the system of linear 00097 *> equations. 00098 *> \endverbatim 00099 *> 00100 *> \param[in] LDX 00101 *> \verbatim 00102 *> LDX is INTEGER 00103 *> The leading dimension of the array X. LDX >= max(1,N). 00104 *> \endverbatim 00105 *> 00106 *> \param[in] B 00107 *> \verbatim 00108 *> B is COMPLEX array, dimension (LDB,NRHS) 00109 *> The right hand side vectors for the system of linear 00110 *> equations. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] LDB 00114 *> \verbatim 00115 *> LDB is INTEGER 00116 *> The leading dimension of the array B. LDB >= max(1,N). 00117 *> \endverbatim 00118 *> 00119 *> \param[out] WORK 00120 *> \verbatim 00121 *> WORK is COMPLEX array, dimension (N) 00122 *> \endverbatim 00123 *> 00124 *> \param[out] RWORK 00125 *> \verbatim 00126 *> RWORK is REAL array, dimension (N) 00127 *> \endverbatim 00128 *> 00129 *> \param[out] RESID 00130 *> \verbatim 00131 *> RESID is REAL 00132 *> The maximum over the number of right hand sides of 00133 *> norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ). 00134 *> \endverbatim 00135 * 00136 * Authors: 00137 * ======== 00138 * 00139 *> \author Univ. of Tennessee 00140 *> \author Univ. of California Berkeley 00141 *> \author Univ. of Colorado Denver 00142 *> \author NAG Ltd. 00143 * 00144 *> \date November 2011 00145 * 00146 *> \ingroup complex_lin 00147 * 00148 * ===================================================================== 00149 SUBROUTINE CTPT02( UPLO, TRANS, DIAG, N, NRHS, AP, X, LDX, B, LDB, 00150 $ WORK, RWORK, RESID ) 00151 * 00152 * -- LAPACK test routine (version 3.4.0) -- 00153 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00154 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00155 * November 2011 00156 * 00157 * .. Scalar Arguments .. 00158 CHARACTER DIAG, TRANS, UPLO 00159 INTEGER LDB, LDX, N, NRHS 00160 REAL RESID 00161 * .. 00162 * .. Array Arguments .. 00163 REAL RWORK( * ) 00164 COMPLEX AP( * ), B( LDB, * ), WORK( * ), X( LDX, * ) 00165 * .. 00166 * 00167 * ===================================================================== 00168 * 00169 * .. Parameters .. 00170 REAL ZERO, ONE 00171 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00172 * .. 00173 * .. Local Scalars .. 00174 INTEGER J 00175 REAL ANORM, BNORM, EPS, XNORM 00176 * .. 00177 * .. External Functions .. 00178 LOGICAL LSAME 00179 REAL CLANTP, SCASUM, SLAMCH 00180 EXTERNAL LSAME, CLANTP, SCASUM, SLAMCH 00181 * .. 00182 * .. External Subroutines .. 00183 EXTERNAL CAXPY, CCOPY, CTPMV 00184 * .. 00185 * .. Intrinsic Functions .. 00186 INTRINSIC CMPLX, MAX 00187 * .. 00188 * .. Executable Statements .. 00189 * 00190 * Quick exit if N = 0 or NRHS = 0 00191 * 00192 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN 00193 RESID = ZERO 00194 RETURN 00195 END IF 00196 * 00197 * Compute the 1-norm of A or A**H. 00198 * 00199 IF( LSAME( TRANS, 'N' ) ) THEN 00200 ANORM = CLANTP( '1', UPLO, DIAG, N, AP, RWORK ) 00201 ELSE 00202 ANORM = CLANTP( 'I', UPLO, DIAG, N, AP, RWORK ) 00203 END IF 00204 * 00205 * Exit with RESID = 1/EPS if ANORM = 0. 00206 * 00207 EPS = SLAMCH( 'Epsilon' ) 00208 IF( ANORM.LE.ZERO ) THEN 00209 RESID = ONE / EPS 00210 RETURN 00211 END IF 00212 * 00213 * Compute the maximum over the number of right hand sides of 00214 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ). 00215 * 00216 RESID = ZERO 00217 DO 10 J = 1, NRHS 00218 CALL CCOPY( N, X( 1, J ), 1, WORK, 1 ) 00219 CALL CTPMV( UPLO, TRANS, DIAG, N, AP, WORK, 1 ) 00220 CALL CAXPY( N, CMPLX( -ONE ), B( 1, J ), 1, WORK, 1 ) 00221 BNORM = SCASUM( N, WORK, 1 ) 00222 XNORM = SCASUM( N, X( 1, J ), 1 ) 00223 IF( XNORM.LE.ZERO ) THEN 00224 RESID = ONE / EPS 00225 ELSE 00226 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) 00227 END IF 00228 10 CONTINUE 00229 * 00230 RETURN 00231 * 00232 * End of CTPT02 00233 * 00234 END