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