![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZTPT01 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 ZTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID ) 00012 * 00013 * .. Scalar Arguments .. 00014 * CHARACTER DIAG, UPLO 00015 * INTEGER N 00016 * DOUBLE PRECISION RCOND, RESID 00017 * .. 00018 * .. Array Arguments .. 00019 * DOUBLE PRECISION RWORK( * ) 00020 * COMPLEX*16 AINVP( * ), AP( * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> ZTPT01 computes the residual for a triangular matrix A times its 00030 *> inverse when A is stored in packed format: 00031 *> RESID = norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ), 00032 *> where EPS is the machine epsilon. 00033 *> \endverbatim 00034 * 00035 * Arguments: 00036 * ========== 00037 * 00038 *> \param[in] UPLO 00039 *> \verbatim 00040 *> UPLO is CHARACTER*1 00041 *> Specifies whether the matrix A is upper or lower triangular. 00042 *> = 'U': Upper triangular 00043 *> = 'L': Lower triangular 00044 *> \endverbatim 00045 *> 00046 *> \param[in] DIAG 00047 *> \verbatim 00048 *> DIAG is CHARACTER*1 00049 *> Specifies whether or not the matrix A is unit triangular. 00050 *> = 'N': Non-unit triangular 00051 *> = 'U': Unit triangular 00052 *> \endverbatim 00053 *> 00054 *> \param[in] N 00055 *> \verbatim 00056 *> N is INTEGER 00057 *> The order of the matrix A. N >= 0. 00058 *> \endverbatim 00059 *> 00060 *> \param[in] AP 00061 *> \verbatim 00062 *> AP is COMPLEX*16 array, dimension (N*(N+1)/2) 00063 *> The original upper or lower triangular matrix A, packed 00064 *> columnwise in a linear array. The j-th column of A is stored 00065 *> in the array AP as follows: 00066 *> if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j; 00067 *> if UPLO = 'L', 00068 *> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] AINVP 00072 *> \verbatim 00073 *> AINVP is COMPLEX*16 array, dimension (N*(N+1)/2) 00074 *> On entry, the (triangular) inverse of the matrix A, packed 00075 *> columnwise in a linear array as in AP. 00076 *> On exit, the contents of AINVP are destroyed. 00077 *> \endverbatim 00078 *> 00079 *> \param[out] RCOND 00080 *> \verbatim 00081 *> RCOND is DOUBLE PRECISION 00082 *> The reciprocal condition number of A, computed as 00083 *> 1/(norm(A) * norm(AINV)). 00084 *> \endverbatim 00085 *> 00086 *> \param[out] RWORK 00087 *> \verbatim 00088 *> RWORK is DOUBLE PRECISION array, dimension (N) 00089 *> \endverbatim 00090 *> 00091 *> \param[out] RESID 00092 *> \verbatim 00093 *> RESID is DOUBLE PRECISION 00094 *> norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ) 00095 *> \endverbatim 00096 * 00097 * Authors: 00098 * ======== 00099 * 00100 *> \author Univ. of Tennessee 00101 *> \author Univ. of California Berkeley 00102 *> \author Univ. of Colorado Denver 00103 *> \author NAG Ltd. 00104 * 00105 *> \date November 2011 00106 * 00107 *> \ingroup complex16_lin 00108 * 00109 * ===================================================================== 00110 SUBROUTINE ZTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID ) 00111 * 00112 * -- LAPACK test routine (version 3.4.0) -- 00113 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00114 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00115 * November 2011 00116 * 00117 * .. Scalar Arguments .. 00118 CHARACTER DIAG, UPLO 00119 INTEGER N 00120 DOUBLE PRECISION RCOND, RESID 00121 * .. 00122 * .. Array Arguments .. 00123 DOUBLE PRECISION RWORK( * ) 00124 COMPLEX*16 AINVP( * ), AP( * ) 00125 * .. 00126 * 00127 * ===================================================================== 00128 * 00129 * .. Parameters .. 00130 DOUBLE PRECISION ZERO, ONE 00131 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00132 * .. 00133 * .. Local Scalars .. 00134 LOGICAL UNITD 00135 INTEGER J, JC 00136 DOUBLE PRECISION AINVNM, ANORM, EPS 00137 * .. 00138 * .. External Functions .. 00139 LOGICAL LSAME 00140 DOUBLE PRECISION DLAMCH, ZLANTP 00141 EXTERNAL LSAME, DLAMCH, ZLANTP 00142 * .. 00143 * .. External Subroutines .. 00144 EXTERNAL ZTPMV 00145 * .. 00146 * .. Intrinsic Functions .. 00147 INTRINSIC DBLE 00148 * .. 00149 * .. Executable Statements .. 00150 * 00151 * Quick exit if N = 0. 00152 * 00153 IF( N.LE.0 ) THEN 00154 RCOND = ONE 00155 RESID = ZERO 00156 RETURN 00157 END IF 00158 * 00159 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00160 * 00161 EPS = DLAMCH( 'Epsilon' ) 00162 ANORM = ZLANTP( '1', UPLO, DIAG, N, AP, RWORK ) 00163 AINVNM = ZLANTP( '1', UPLO, DIAG, N, AINVP, RWORK ) 00164 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00165 RCOND = ZERO 00166 RESID = ONE / EPS 00167 RETURN 00168 END IF 00169 RCOND = ( ONE / ANORM ) / AINVNM 00170 * 00171 * Compute A * AINV, overwriting AINV. 00172 * 00173 UNITD = LSAME( DIAG, 'U' ) 00174 IF( LSAME( UPLO, 'U' ) ) THEN 00175 JC = 1 00176 DO 10 J = 1, N 00177 IF( UNITD ) 00178 $ AINVP( JC+J-1 ) = ONE 00179 * 00180 * Form the j-th column of A*AINV. 00181 * 00182 CALL ZTPMV( 'Upper', 'No transpose', DIAG, J, AP, 00183 $ AINVP( JC ), 1 ) 00184 * 00185 * Subtract 1 from the diagonal to form A*AINV - I. 00186 * 00187 AINVP( JC+J-1 ) = AINVP( JC+J-1 ) - ONE 00188 JC = JC + J 00189 10 CONTINUE 00190 ELSE 00191 JC = 1 00192 DO 20 J = 1, N 00193 IF( UNITD ) 00194 $ AINVP( JC ) = ONE 00195 * 00196 * Form the j-th column of A*AINV. 00197 * 00198 CALL ZTPMV( 'Lower', 'No transpose', DIAG, N-J+1, AP( JC ), 00199 $ AINVP( JC ), 1 ) 00200 * 00201 * Subtract 1 from the diagonal to form A*AINV - I. 00202 * 00203 AINVP( JC ) = AINVP( JC ) - ONE 00204 JC = JC + N - J + 1 00205 20 CONTINUE 00206 END IF 00207 * 00208 * Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS) 00209 * 00210 RESID = ZLANTP( '1', UPLO, 'Non-unit', N, AINVP, RWORK ) 00211 * 00212 RESID = ( ( RESID*RCOND ) / DBLE( N ) ) / EPS 00213 * 00214 RETURN 00215 * 00216 * End of ZTPT01 00217 * 00218 END