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