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