3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
11 * SUBROUTINE ZTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID )
13 * .. Scalar Arguments ..
14 * CHARACTER DIAG, UPLO
16 * DOUBLE PRECISION RCOND, RESID
18 * .. Array Arguments ..
19 * DOUBLE PRECISION RWORK( * )
20 * COMPLEX*16 AINVP( * ), AP( * )
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.
40 *> UPLO is CHARACTER*1
41 *> Specifies whether the matrix A is upper or lower triangular.
42 *> = 'U': Upper triangular
43 *> = 'L': Lower triangular
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
57 *> The order of the matrix A. N >= 0.
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;
68 *> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
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.
81 *> RCOND is DOUBLE PRECISION
82 *> The reciprocal condition number of A, computed as
83 *> 1/(norm(A) * norm(AINV)).
88 *> RWORK is DOUBLE PRECISION array, dimension (N)
93 *> RESID is DOUBLE PRECISION
94 *> norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS )
100 *> \author Univ. of Tennessee
101 *> \author Univ. of California Berkeley
102 *> \author Univ. of Colorado Denver
105 *> \date November 2011
107 *> \ingroup complex16_lin
109 * =====================================================================
110 SUBROUTINE ZTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID )
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..--
117 * .. Scalar Arguments ..
120 DOUBLE PRECISION RCOND, RESID
122 * .. Array Arguments ..
123 DOUBLE PRECISION RWORK( * )
124 COMPLEX*16 AINVP( * ), AP( * )
127 * =====================================================================
130 DOUBLE PRECISION ZERO, ONE
131 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
133 * .. Local Scalars ..
136 DOUBLE PRECISION AINVNM, ANORM, EPS
138 * .. External Functions ..
140 DOUBLE PRECISION DLAMCH, ZLANTP
141 EXTERNAL LSAME, DLAMCH, ZLANTP
143 * .. External Subroutines ..
146 * .. Intrinsic Functions ..
149 * .. Executable Statements ..
151 * Quick exit if N = 0.
159 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
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
169 RCOND = ( ONE / ANORM ) / AINVNM
171 * Compute A * AINV, overwriting AINV.
173 UNITD = LSAME( DIAG, 'U' )
174 IF( LSAME( UPLO, 'U' ) ) THEN
178 $ AINVP( JC+J-1 ) = ONE
180 * Form the j-th column of A*AINV.
182 CALL ZTPMV( 'Upper', 'No transpose', DIAG, J, AP,
185 * Subtract 1 from the diagonal to form A*AINV - I.
187 AINVP( JC+J-1 ) = AINVP( JC+J-1 ) - ONE
196 * Form the j-th column of A*AINV.
198 CALL ZTPMV( 'Lower', 'No transpose', DIAG, N-J+1, AP( JC ),
201 * Subtract 1 from the diagonal to form A*AINV - I.
203 AINVP( JC ) = AINVP( JC ) - ONE
208 * Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
210 RESID = ZLANTP( '1', UPLO, 'Non-unit', N, AINVP, RWORK )
212 RESID = ( ( RESID*RCOND ) / DBLE( N ) ) / EPS