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