![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZTRT02 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 ZTRT02( UPLO, TRANS, DIAG, N, NRHS, A, LDA, X, LDX, B, 00012 * LDB, WORK, RWORK, RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER DIAG, TRANS, UPLO 00016 * INTEGER LDA, LDB, LDX, N, NRHS 00017 * DOUBLE PRECISION RESID 00018 * .. 00019 * .. Array Arguments .. 00020 * DOUBLE PRECISION RWORK( * ) 00021 * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ), 00022 * $ X( LDX, * ) 00023 * .. 00024 * 00025 * 00026 *> \par Purpose: 00027 * ============= 00028 *> 00029 *> \verbatim 00030 *> 00031 *> ZTRT02 computes the residual for the computed solution to a 00032 *> triangular system of linear equations A*x = b, A**T *x = b, 00033 *> or A**H *x = b. Here A is a triangular matrix, A**T is the transpose 00034 *> of A, A**H is the conjugate transpose of A, and x and b are N by NRHS 00035 *> matrices. The test ratio is the maximum over the number of right 00036 *> 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] A 00083 *> \verbatim 00084 *> A is COMPLEX*16 array, dimension (LDA,N) 00085 *> The triangular matrix A. If UPLO = 'U', the leading n by n 00086 *> upper triangular part of the array A contains the upper 00087 *> triangular matrix, and the strictly lower triangular part of 00088 *> A is not referenced. If UPLO = 'L', the leading n by n lower 00089 *> triangular part of the array A contains the lower triangular 00090 *> matrix, and the strictly upper triangular part of A is not 00091 *> referenced. If DIAG = 'U', the diagonal elements of A are 00092 *> also not referenced and are assumed to be 1. 00093 *> \endverbatim 00094 *> 00095 *> \param[in] LDA 00096 *> \verbatim 00097 *> LDA is INTEGER 00098 *> The leading dimension of the array A. LDA >= max(1,N). 00099 *> \endverbatim 00100 *> 00101 *> \param[in] X 00102 *> \verbatim 00103 *> X is COMPLEX*16 array, dimension (LDX,NRHS) 00104 *> The computed solution vectors for the system of linear 00105 *> equations. 00106 *> \endverbatim 00107 *> 00108 *> \param[in] LDX 00109 *> \verbatim 00110 *> LDX is INTEGER 00111 *> The leading dimension of the array X. LDX >= max(1,N). 00112 *> \endverbatim 00113 *> 00114 *> \param[in] B 00115 *> \verbatim 00116 *> B is COMPLEX*16 array, dimension (LDB,NRHS) 00117 *> The right hand side vectors for the system of linear 00118 *> equations. 00119 *> \endverbatim 00120 *> 00121 *> \param[in] LDB 00122 *> \verbatim 00123 *> LDB is INTEGER 00124 *> The leading dimension of the array B. LDB >= max(1,N). 00125 *> \endverbatim 00126 *> 00127 *> \param[out] WORK 00128 *> \verbatim 00129 *> WORK is COMPLEX*16 array, dimension (N) 00130 *> \endverbatim 00131 *> 00132 *> \param[out] RWORK 00133 *> \verbatim 00134 *> RWORK is DOUBLE PRECISION array, dimension (N) 00135 *> \endverbatim 00136 *> 00137 *> \param[out] RESID 00138 *> \verbatim 00139 *> RESID is DOUBLE PRECISION 00140 *> The maximum over the number of right hand sides of 00141 *> norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ). 00142 *> \endverbatim 00143 * 00144 * Authors: 00145 * ======== 00146 * 00147 *> \author Univ. of Tennessee 00148 *> \author Univ. of California Berkeley 00149 *> \author Univ. of Colorado Denver 00150 *> \author NAG Ltd. 00151 * 00152 *> \date November 2011 00153 * 00154 *> \ingroup complex16_lin 00155 * 00156 * ===================================================================== 00157 SUBROUTINE ZTRT02( UPLO, TRANS, DIAG, N, NRHS, A, LDA, X, LDX, B, 00158 $ LDB, WORK, RWORK, RESID ) 00159 * 00160 * -- LAPACK test routine (version 3.4.0) -- 00161 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00162 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00163 * November 2011 00164 * 00165 * .. Scalar Arguments .. 00166 CHARACTER DIAG, TRANS, UPLO 00167 INTEGER LDA, LDB, LDX, N, NRHS 00168 DOUBLE PRECISION RESID 00169 * .. 00170 * .. Array Arguments .. 00171 DOUBLE PRECISION RWORK( * ) 00172 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ), 00173 $ X( LDX, * ) 00174 * .. 00175 * 00176 * ===================================================================== 00177 * 00178 * .. Parameters .. 00179 DOUBLE PRECISION ZERO, ONE 00180 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00181 * .. 00182 * .. Local Scalars .. 00183 INTEGER J 00184 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM 00185 * .. 00186 * .. External Functions .. 00187 LOGICAL LSAME 00188 DOUBLE PRECISION DLAMCH, DZASUM, ZLANTR 00189 EXTERNAL LSAME, DLAMCH, DZASUM, ZLANTR 00190 * .. 00191 * .. External Subroutines .. 00192 EXTERNAL ZAXPY, ZCOPY, ZTRMV 00193 * .. 00194 * .. Intrinsic Functions .. 00195 INTRINSIC DCMPLX, MAX 00196 * .. 00197 * .. Executable Statements .. 00198 * 00199 * Quick exit if N = 0 or NRHS = 0 00200 * 00201 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN 00202 RESID = ZERO 00203 RETURN 00204 END IF 00205 * 00206 * Compute the 1-norm of A or A**H. 00207 * 00208 IF( LSAME( TRANS, 'N' ) ) THEN 00209 ANORM = ZLANTR( '1', UPLO, DIAG, N, N, A, LDA, RWORK ) 00210 ELSE 00211 ANORM = ZLANTR( 'I', UPLO, DIAG, N, N, A, LDA, RWORK ) 00212 END IF 00213 * 00214 * Exit with RESID = 1/EPS if ANORM = 0. 00215 * 00216 EPS = DLAMCH( 'Epsilon' ) 00217 IF( ANORM.LE.ZERO ) THEN 00218 RESID = ONE / EPS 00219 RETURN 00220 END IF 00221 * 00222 * Compute the maximum over the number of right hand sides of 00223 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ) 00224 * 00225 RESID = ZERO 00226 DO 10 J = 1, NRHS 00227 CALL ZCOPY( N, X( 1, J ), 1, WORK, 1 ) 00228 CALL ZTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 ) 00229 CALL ZAXPY( N, DCMPLX( -ONE ), B( 1, J ), 1, WORK, 1 ) 00230 BNORM = DZASUM( N, WORK, 1 ) 00231 XNORM = DZASUM( N, X( 1, J ), 1 ) 00232 IF( XNORM.LE.ZERO ) THEN 00233 RESID = ONE / EPS 00234 ELSE 00235 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) 00236 END IF 00237 10 CONTINUE 00238 * 00239 RETURN 00240 * 00241 * End of ZTRT02 00242 * 00243 END