1*> \brief \b ZTPT01 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE ZTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID ) 12* 13* .. Scalar Arguments .. 14* CHARACTER DIAG, UPLO 15* INTEGER N 16* DOUBLE PRECISION RCOND, RESID 17* .. 18* .. Array Arguments .. 19* DOUBLE PRECISION RWORK( * ) 20* COMPLEX*16 AINVP( * ), AP( * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> ZTPT01 computes the residual for a triangular matrix A times its 30*> inverse when A is stored in packed format: 31*> RESID = norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ), 32*> where EPS is the machine epsilon. 33*> \endverbatim 34* 35* Arguments: 36* ========== 37* 38*> \param[in] UPLO 39*> \verbatim 40*> UPLO is CHARACTER*1 41*> Specifies whether the matrix A is upper or lower triangular. 42*> = 'U': Upper triangular 43*> = 'L': Lower triangular 44*> \endverbatim 45*> 46*> \param[in] DIAG 47*> \verbatim 48*> DIAG is CHARACTER*1 49*> Specifies whether or not the matrix A is unit triangular. 50*> = 'N': Non-unit triangular 51*> = 'U': Unit triangular 52*> \endverbatim 53*> 54*> \param[in] N 55*> \verbatim 56*> N is INTEGER 57*> The order of the matrix A. N >= 0. 58*> \endverbatim 59*> 60*> \param[in] AP 61*> \verbatim 62*> AP is COMPLEX*16 array, dimension (N*(N+1)/2) 63*> The original upper or lower triangular matrix A, packed 64*> columnwise in a linear array. The j-th column of A is stored 65*> in the array AP as follows: 66*> if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j; 67*> if UPLO = 'L', 68*> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n. 69*> \endverbatim 70*> 71*> \param[in] AINVP 72*> \verbatim 73*> AINVP is COMPLEX*16 array, dimension (N*(N+1)/2) 74*> On entry, the (triangular) inverse of the matrix A, packed 75*> columnwise in a linear array as in AP. 76*> On exit, the contents of AINVP are destroyed. 77*> \endverbatim 78*> 79*> \param[out] RCOND 80*> \verbatim 81*> RCOND is DOUBLE PRECISION 82*> The reciprocal condition number of A, computed as 83*> 1/(norm(A) * norm(AINV)). 84*> \endverbatim 85*> 86*> \param[out] RWORK 87*> \verbatim 88*> RWORK is DOUBLE PRECISION array, dimension (N) 89*> \endverbatim 90*> 91*> \param[out] RESID 92*> \verbatim 93*> RESID is DOUBLE PRECISION 94*> norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ) 95*> \endverbatim 96* 97* Authors: 98* ======== 99* 100*> \author Univ. of Tennessee 101*> \author Univ. of California Berkeley 102*> \author Univ. of Colorado Denver 103*> \author NAG Ltd. 104* 105*> \date November 2011 106* 107*> \ingroup complex16_lin 108* 109* ===================================================================== 110 SUBROUTINE ZTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID ) 111* 112* -- LAPACK test routine (version 3.4.0) -- 113* -- LAPACK is a software package provided by Univ. of Tennessee, -- 114* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 115* November 2011 116* 117* .. Scalar Arguments .. 118 CHARACTER DIAG, UPLO 119 INTEGER N 120 DOUBLE PRECISION RCOND, RESID 121* .. 122* .. Array Arguments .. 123 DOUBLE PRECISION RWORK( * ) 124 COMPLEX*16 AINVP( * ), AP( * ) 125* .. 126* 127* ===================================================================== 128* 129* .. Parameters .. 130 DOUBLE PRECISION ZERO, ONE 131 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 132* .. 133* .. Local Scalars .. 134 LOGICAL UNITD 135 INTEGER J, JC 136 DOUBLE PRECISION AINVNM, ANORM, EPS 137* .. 138* .. External Functions .. 139 LOGICAL LSAME 140 DOUBLE PRECISION DLAMCH, ZLANTP 141 EXTERNAL LSAME, DLAMCH, ZLANTP 142* .. 143* .. External Subroutines .. 144 EXTERNAL ZTPMV 145* .. 146* .. Intrinsic Functions .. 147 INTRINSIC DBLE 148* .. 149* .. Executable Statements .. 150* 151* Quick exit if N = 0. 152* 153 IF( N.LE.0 ) THEN 154 RCOND = ONE 155 RESID = ZERO 156 RETURN 157 END IF 158* 159* Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 160* 161 EPS = DLAMCH( 'Epsilon' ) 162 ANORM = ZLANTP( '1', UPLO, DIAG, N, AP, RWORK ) 163 AINVNM = ZLANTP( '1', UPLO, DIAG, N, AINVP, RWORK ) 164 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 165 RCOND = ZERO 166 RESID = ONE / EPS 167 RETURN 168 END IF 169 RCOND = ( ONE / ANORM ) / AINVNM 170* 171* Compute A * AINV, overwriting AINV. 172* 173 UNITD = LSAME( DIAG, 'U' ) 174 IF( LSAME( UPLO, 'U' ) ) THEN 175 JC = 1 176 DO 10 J = 1, N 177 IF( UNITD ) 178 $ AINVP( JC+J-1 ) = ONE 179* 180* Form the j-th column of A*AINV. 181* 182 CALL ZTPMV( 'Upper', 'No transpose', DIAG, J, AP, 183 $ AINVP( JC ), 1 ) 184* 185* Subtract 1 from the diagonal to form A*AINV - I. 186* 187 AINVP( JC+J-1 ) = AINVP( JC+J-1 ) - ONE 188 JC = JC + J 189 10 CONTINUE 190 ELSE 191 JC = 1 192 DO 20 J = 1, N 193 IF( UNITD ) 194 $ AINVP( JC ) = ONE 195* 196* Form the j-th column of A*AINV. 197* 198 CALL ZTPMV( 'Lower', 'No transpose', DIAG, N-J+1, AP( JC ), 199 $ AINVP( JC ), 1 ) 200* 201* Subtract 1 from the diagonal to form A*AINV - I. 202* 203 AINVP( JC ) = AINVP( JC ) - ONE 204 JC = JC + N - J + 1 205 20 CONTINUE 206 END IF 207* 208* Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS) 209* 210 RESID = ZLANTP( '1', UPLO, 'Non-unit', N, AINVP, RWORK ) 211* 212 RESID = ( ( RESID*RCOND ) / DBLE( N ) ) / EPS 213* 214 RETURN 215* 216* End of ZTPT01 217* 218 END 219