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