![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CPPT03 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 CPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND, 00012 * RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER UPLO 00016 * INTEGER LDWORK, N 00017 * REAL RCOND, RESID 00018 * .. 00019 * .. Array Arguments .. 00020 * REAL RWORK( * ) 00021 * COMPLEX A( * ), AINV( * ), WORK( LDWORK, * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> CPPT03 computes the residual for a Hermitian packed matrix times its 00031 *> inverse: 00032 *> norm( I - A*AINV ) / ( 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 upper or lower triangular part of the 00043 *> Hermitian matrix A is stored: 00044 *> = 'U': Upper triangular 00045 *> = 'L': Lower triangular 00046 *> \endverbatim 00047 *> 00048 *> \param[in] N 00049 *> \verbatim 00050 *> N is INTEGER 00051 *> The number of rows and columns of the matrix A. N >= 0. 00052 *> \endverbatim 00053 *> 00054 *> \param[in] A 00055 *> \verbatim 00056 *> A is COMPLEX array, dimension (N*(N+1)/2) 00057 *> The original Hermitian matrix A, stored as a packed 00058 *> triangular matrix. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] AINV 00062 *> \verbatim 00063 *> AINV is COMPLEX array, dimension (N*(N+1)/2) 00064 *> The (Hermitian) inverse of the matrix A, stored as a packed 00065 *> triangular matrix. 00066 *> \endverbatim 00067 *> 00068 *> \param[out] WORK 00069 *> \verbatim 00070 *> WORK is COMPLEX array, dimension (LDWORK,N) 00071 *> \endverbatim 00072 *> 00073 *> \param[in] LDWORK 00074 *> \verbatim 00075 *> LDWORK is INTEGER 00076 *> The leading dimension of the array WORK. LDWORK >= max(1,N). 00077 *> \endverbatim 00078 *> 00079 *> \param[out] RWORK 00080 *> \verbatim 00081 *> RWORK is REAL array, dimension (N) 00082 *> \endverbatim 00083 *> 00084 *> \param[out] RCOND 00085 *> \verbatim 00086 *> RCOND is REAL 00087 *> The reciprocal of the condition number of A, computed as 00088 *> ( 1/norm(A) ) / norm(AINV). 00089 *> \endverbatim 00090 *> 00091 *> \param[out] RESID 00092 *> \verbatim 00093 *> RESID is REAL 00094 *> norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS ) 00095 *> \endverbatim 00096 * 00097 * Authors: 00098 * ======== 00099 * 00100 *> \author Univ. of Tennessee 00101 *> \author Univ. of California Berkeley 00102 *> \author Univ. of Colorado Denver 00103 *> \author NAG Ltd. 00104 * 00105 *> \date November 2011 00106 * 00107 *> \ingroup complex_lin 00108 * 00109 * ===================================================================== 00110 SUBROUTINE CPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND, 00111 $ RESID ) 00112 * 00113 * -- LAPACK test routine (version 3.4.0) -- 00114 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00115 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00116 * November 2011 00117 * 00118 * .. Scalar Arguments .. 00119 CHARACTER UPLO 00120 INTEGER LDWORK, N 00121 REAL RCOND, RESID 00122 * .. 00123 * .. Array Arguments .. 00124 REAL RWORK( * ) 00125 COMPLEX A( * ), AINV( * ), WORK( LDWORK, * ) 00126 * .. 00127 * 00128 * ===================================================================== 00129 * 00130 * .. Parameters .. 00131 REAL ZERO, ONE 00132 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00133 COMPLEX CZERO, CONE 00134 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00135 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00136 * .. 00137 * .. Local Scalars .. 00138 INTEGER I, J, JJ 00139 REAL AINVNM, ANORM, EPS 00140 * .. 00141 * .. External Functions .. 00142 LOGICAL LSAME 00143 REAL CLANGE, CLANHP, SLAMCH 00144 EXTERNAL LSAME, CLANGE, CLANHP, SLAMCH 00145 * .. 00146 * .. Intrinsic Functions .. 00147 INTRINSIC CONJG, REAL 00148 * .. 00149 * .. External Subroutines .. 00150 EXTERNAL CCOPY, CHPMV 00151 * .. 00152 * .. Executable Statements .. 00153 * 00154 * Quick exit if N = 0. 00155 * 00156 IF( N.LE.0 ) THEN 00157 RCOND = ONE 00158 RESID = ZERO 00159 RETURN 00160 END IF 00161 * 00162 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00163 * 00164 EPS = SLAMCH( 'Epsilon' ) 00165 ANORM = CLANHP( '1', UPLO, N, A, RWORK ) 00166 AINVNM = CLANHP( '1', UPLO, N, AINV, RWORK ) 00167 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00168 RCOND = ZERO 00169 RESID = ONE / EPS 00170 RETURN 00171 END IF 00172 RCOND = ( ONE/ANORM ) / AINVNM 00173 * 00174 * UPLO = 'U': 00175 * Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and 00176 * expand it to a full matrix, then multiply by A one column at a 00177 * time, moving the result one column to the left. 00178 * 00179 IF( LSAME( UPLO, 'U' ) ) THEN 00180 * 00181 * Copy AINV 00182 * 00183 JJ = 1 00184 DO 20 J = 1, N - 1 00185 CALL CCOPY( J, AINV( JJ ), 1, WORK( 1, J+1 ), 1 ) 00186 DO 10 I = 1, J - 1 00187 WORK( J, I+1 ) = CONJG( AINV( JJ+I-1 ) ) 00188 10 CONTINUE 00189 JJ = JJ + J 00190 20 CONTINUE 00191 JJ = ( ( N-1 )*N ) / 2 + 1 00192 DO 30 I = 1, N - 1 00193 WORK( N, I+1 ) = CONJG( AINV( JJ+I-1 ) ) 00194 30 CONTINUE 00195 * 00196 * Multiply by A 00197 * 00198 DO 40 J = 1, N - 1 00199 CALL CHPMV( 'Upper', N, -CONE, A, WORK( 1, J+1 ), 1, CZERO, 00200 $ WORK( 1, J ), 1 ) 00201 40 CONTINUE 00202 CALL CHPMV( 'Upper', N, -CONE, A, AINV( JJ ), 1, CZERO, 00203 $ WORK( 1, N ), 1 ) 00204 * 00205 * UPLO = 'L': 00206 * Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1) 00207 * and multiply by A, moving each column to the right. 00208 * 00209 ELSE 00210 * 00211 * Copy AINV 00212 * 00213 DO 50 I = 1, N - 1 00214 WORK( 1, I ) = CONJG( AINV( I+1 ) ) 00215 50 CONTINUE 00216 JJ = N + 1 00217 DO 70 J = 2, N 00218 CALL CCOPY( N-J+1, AINV( JJ ), 1, WORK( J, J-1 ), 1 ) 00219 DO 60 I = 1, N - J 00220 WORK( J, J+I-1 ) = CONJG( AINV( JJ+I ) ) 00221 60 CONTINUE 00222 JJ = JJ + N - J + 1 00223 70 CONTINUE 00224 * 00225 * Multiply by A 00226 * 00227 DO 80 J = N, 2, -1 00228 CALL CHPMV( 'Lower', N, -CONE, A, WORK( 1, J-1 ), 1, CZERO, 00229 $ WORK( 1, J ), 1 ) 00230 80 CONTINUE 00231 CALL CHPMV( 'Lower', N, -CONE, A, AINV( 1 ), 1, CZERO, 00232 $ WORK( 1, 1 ), 1 ) 00233 * 00234 END IF 00235 * 00236 * Add the identity matrix to WORK . 00237 * 00238 DO 90 I = 1, N 00239 WORK( I, I ) = WORK( I, I ) + CONE 00240 90 CONTINUE 00241 * 00242 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 00243 * 00244 RESID = CLANGE( '1', N, N, WORK, LDWORK, RWORK ) 00245 * 00246 RESID = ( ( RESID*RCOND )/EPS ) / REAL( N ) 00247 * 00248 RETURN 00249 * 00250 * End of CPPT03 00251 * 00252 END