00001 SUBROUTINE DTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, WORK, RESID ) 00002 * 00003 * -- LAPACK test routine (version 3.1) -- 00004 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00005 * November 2006 00006 * 00007 * .. Scalar Arguments .. 00008 CHARACTER DIAG, UPLO 00009 INTEGER N 00010 DOUBLE PRECISION RCOND, RESID 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION AINVP( * ), AP( * ), WORK( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * DTPT01 computes the residual for a triangular matrix A times its 00020 * inverse when A is stored in packed format: 00021 * RESID = norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ), 00022 * where EPS is the machine epsilon. 00023 * 00024 * Arguments 00025 * ========== 00026 * 00027 * UPLO (input) CHARACTER*1 00028 * Specifies whether the matrix A is upper or lower triangular. 00029 * = 'U': Upper triangular 00030 * = 'L': Lower triangular 00031 * 00032 * DIAG (input) CHARACTER*1 00033 * Specifies whether or not the matrix A is unit triangular. 00034 * = 'N': Non-unit triangular 00035 * = 'U': Unit triangular 00036 * 00037 * N (input) INTEGER 00038 * The order of the matrix A. N >= 0. 00039 * 00040 * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2) 00041 * The original upper or lower triangular matrix A, packed 00042 * columnwise in a linear array. The j-th column of A is stored 00043 * in the array AP as follows: 00044 * if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j; 00045 * if UPLO = 'L', 00046 * AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n. 00047 * 00048 * AINVP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) 00049 * On entry, the (triangular) inverse of the matrix A, packed 00050 * columnwise in a linear array as in AP. 00051 * On exit, the contents of AINVP are destroyed. 00052 * 00053 * RCOND (output) DOUBLE PRECISION 00054 * The reciprocal condition number of A, computed as 00055 * 1/(norm(A) * norm(AINV)). 00056 * 00057 * WORK (workspace) DOUBLE PRECISION array, dimension (N) 00058 * 00059 * RESID (output) DOUBLE PRECISION 00060 * norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ) 00061 * 00062 * ===================================================================== 00063 * 00064 * .. Parameters .. 00065 DOUBLE PRECISION ZERO, ONE 00066 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00067 * .. 00068 * .. Local Scalars .. 00069 LOGICAL UNITD 00070 INTEGER J, JC 00071 DOUBLE PRECISION AINVNM, ANORM, EPS 00072 * .. 00073 * .. External Functions .. 00074 LOGICAL LSAME 00075 DOUBLE PRECISION DLAMCH, DLANTP 00076 EXTERNAL LSAME, DLAMCH, DLANTP 00077 * .. 00078 * .. External Subroutines .. 00079 EXTERNAL DTPMV 00080 * .. 00081 * .. Intrinsic Functions .. 00082 INTRINSIC DBLE 00083 * .. 00084 * .. Executable Statements .. 00085 * 00086 * Quick exit if N = 0. 00087 * 00088 IF( N.LE.0 ) THEN 00089 RCOND = ONE 00090 RESID = ZERO 00091 RETURN 00092 END IF 00093 * 00094 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00095 * 00096 EPS = DLAMCH( 'Epsilon' ) 00097 ANORM = DLANTP( '1', UPLO, DIAG, N, AP, WORK ) 00098 AINVNM = DLANTP( '1', UPLO, DIAG, N, AINVP, WORK ) 00099 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00100 RCOND = ZERO 00101 RESID = ONE / EPS 00102 RETURN 00103 END IF 00104 RCOND = ( ONE / ANORM ) / AINVNM 00105 * 00106 * Compute A * AINV, overwriting AINV. 00107 * 00108 UNITD = LSAME( DIAG, 'U' ) 00109 IF( LSAME( UPLO, 'U' ) ) THEN 00110 JC = 1 00111 DO 10 J = 1, N 00112 IF( UNITD ) 00113 $ AINVP( JC+J-1 ) = ONE 00114 * 00115 * Form the j-th column of A*AINV 00116 * 00117 CALL DTPMV( 'Upper', 'No transpose', DIAG, J, AP, 00118 $ AINVP( JC ), 1 ) 00119 * 00120 * Subtract 1 from the diagonal 00121 * 00122 AINVP( JC+J-1 ) = AINVP( JC+J-1 ) - ONE 00123 JC = JC + J 00124 10 CONTINUE 00125 ELSE 00126 JC = 1 00127 DO 20 J = 1, N 00128 IF( UNITD ) 00129 $ AINVP( JC ) = ONE 00130 * 00131 * Form the j-th column of A*AINV 00132 * 00133 CALL DTPMV( 'Lower', 'No transpose', DIAG, N-J+1, AP( JC ), 00134 $ AINVP( JC ), 1 ) 00135 * 00136 * Subtract 1 from the diagonal 00137 * 00138 AINVP( JC ) = AINVP( JC ) - ONE 00139 JC = JC + N - J + 1 00140 20 CONTINUE 00141 END IF 00142 * 00143 * Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS) 00144 * 00145 RESID = DLANTP( '1', UPLO, 'Non-unit', N, AINVP, WORK ) 00146 * 00147 RESID = ( ( RESID*RCOND ) / DBLE( N ) ) / EPS 00148 * 00149 RETURN 00150 * 00151 * End of DTPT01 00152 * 00153 END
1.7.2