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