00001 SUBROUTINE STBT02( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, X, 00002 $ LDX, B, LDB, WORK, RESID ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER DIAG, TRANS, UPLO 00010 INTEGER KD, LDAB, LDB, LDX, N, NRHS 00011 REAL RESID 00012 * .. 00013 * .. Array Arguments .. 00014 REAL AB( LDAB, * ), B( LDB, * ), WORK( * ), 00015 $ X( LDX, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * STBT02 computes the residual for the computed solution to a 00022 * triangular system of linear equations A*x = b or A' *x = b when 00023 * A is a triangular band matrix. Here A' is the transpose of A and 00024 * x and b are N by NRHS matrices. The test ratio is the maximum over 00025 * the number of right hand sides of 00026 * norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ), 00027 * where op(A) denotes A or A' and EPS is the machine epsilon. 00028 * 00029 * Arguments 00030 * ========= 00031 * 00032 * UPLO (input) CHARACTER*1 00033 * Specifies whether the matrix A is upper or lower triangular. 00034 * = 'U': Upper triangular 00035 * = 'L': Lower triangular 00036 * 00037 * TRANS (input) CHARACTER*1 00038 * Specifies the operation applied to A. 00039 * = 'N': A *x = b (No transpose) 00040 * = 'T': A'*x = b (Transpose) 00041 * = 'C': A'*x = b (Conjugate transpose = Transpose) 00042 * 00043 * DIAG (input) CHARACTER*1 00044 * Specifies whether or not the matrix A is unit triangular. 00045 * = 'N': Non-unit triangular 00046 * = 'U': Unit triangular 00047 * 00048 * N (input) INTEGER 00049 * The order of the matrix A. N >= 0. 00050 * 00051 * KD (input) INTEGER 00052 * The number of superdiagonals or subdiagonals of the 00053 * triangular band matrix A. KD >= 0. 00054 * 00055 * NRHS (input) INTEGER 00056 * The number of right hand sides, i.e., the number of columns 00057 * of the matrices X and B. NRHS >= 0. 00058 * 00059 * AB (input) REAL array, dimension (LDAB,N) 00060 * The upper or lower triangular band matrix A, stored in the 00061 * first kd+1 rows of the array. The j-th column of A is stored 00062 * in the j-th column of the array AB as follows: 00063 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 00064 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 00065 * 00066 * LDAB (input) INTEGER 00067 * The leading dimension of the array AB. LDAB >= KD+1. 00068 * 00069 * X (input) REAL array, dimension (LDX,NRHS) 00070 * The computed solution vectors for the system of linear 00071 * equations. 00072 * 00073 * LDX (input) INTEGER 00074 * The leading dimension of the array X. LDX >= max(1,N). 00075 * 00076 * B (input) REAL array, dimension (LDB,NRHS) 00077 * The right hand side vectors for the system of linear 00078 * equations. 00079 * 00080 * LDB (input) INTEGER 00081 * The leading dimension of the array B. LDB >= max(1,N). 00082 * 00083 * WORK (workspace) REAL array, dimension (N) 00084 * 00085 * RESID (output) REAL 00086 * The maximum over the number of right hand sides of 00087 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ). 00088 * 00089 * ===================================================================== 00090 * 00091 * .. Parameters .. 00092 REAL ZERO, ONE 00093 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00094 * .. 00095 * .. Local Scalars .. 00096 INTEGER J 00097 REAL ANORM, BNORM, EPS, XNORM 00098 * .. 00099 * .. External Functions .. 00100 LOGICAL LSAME 00101 REAL SASUM, SLAMCH, SLANTB 00102 EXTERNAL LSAME, SASUM, SLAMCH, SLANTB 00103 * .. 00104 * .. External Subroutines .. 00105 EXTERNAL SAXPY, SCOPY, STBMV 00106 * .. 00107 * .. Intrinsic Functions .. 00108 INTRINSIC MAX 00109 * .. 00110 * .. Executable Statements .. 00111 * 00112 * Quick exit if N = 0 or NRHS = 0 00113 * 00114 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN 00115 RESID = ZERO 00116 RETURN 00117 END IF 00118 * 00119 * Compute the 1-norm of A or A'. 00120 * 00121 IF( LSAME( TRANS, 'N' ) ) THEN 00122 ANORM = SLANTB( '1', UPLO, DIAG, N, KD, AB, LDAB, WORK ) 00123 ELSE 00124 ANORM = SLANTB( 'I', UPLO, DIAG, N, KD, AB, LDAB, WORK ) 00125 END IF 00126 * 00127 * Exit with RESID = 1/EPS if ANORM = 0. 00128 * 00129 EPS = SLAMCH( 'Epsilon' ) 00130 IF( ANORM.LE.ZERO ) THEN 00131 RESID = ONE / EPS 00132 RETURN 00133 END IF 00134 * 00135 * Compute the maximum over the number of right hand sides of 00136 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ). 00137 * 00138 RESID = ZERO 00139 DO 10 J = 1, NRHS 00140 CALL SCOPY( N, X( 1, J ), 1, WORK, 1 ) 00141 CALL STBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 ) 00142 CALL SAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 ) 00143 BNORM = SASUM( N, WORK, 1 ) 00144 XNORM = SASUM( N, X( 1, J ), 1 ) 00145 IF( XNORM.LE.ZERO ) THEN 00146 RESID = ONE / EPS 00147 ELSE 00148 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) 00149 END IF 00150 10 CONTINUE 00151 * 00152 RETURN 00153 * 00154 * End of STBT02 00155 * 00156 END
1.7.2