1*> \brief \b CPPT03 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE CPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND, 12* RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER UPLO 16* INTEGER LDWORK, N 17* REAL RCOND, RESID 18* .. 19* .. Array Arguments .. 20* REAL RWORK( * ) 21* COMPLEX A( * ), AINV( * ), WORK( LDWORK, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> CPPT03 computes the residual for a Hermitian packed matrix times its 31*> inverse: 32*> norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ), 33*> where EPS is the machine epsilon. 34*> \endverbatim 35* 36* Arguments: 37* ========== 38* 39*> \param[in] UPLO 40*> \verbatim 41*> UPLO is CHARACTER*1 42*> Specifies whether the upper or lower triangular part of the 43*> Hermitian matrix A is stored: 44*> = 'U': Upper triangular 45*> = 'L': Lower triangular 46*> \endverbatim 47*> 48*> \param[in] N 49*> \verbatim 50*> N is INTEGER 51*> The number of rows and columns of the matrix A. N >= 0. 52*> \endverbatim 53*> 54*> \param[in] A 55*> \verbatim 56*> A is COMPLEX array, dimension (N*(N+1)/2) 57*> The original Hermitian matrix A, stored as a packed 58*> triangular matrix. 59*> \endverbatim 60*> 61*> \param[in] AINV 62*> \verbatim 63*> AINV is COMPLEX array, dimension (N*(N+1)/2) 64*> The (Hermitian) inverse of the matrix A, stored as a packed 65*> triangular matrix. 66*> \endverbatim 67*> 68*> \param[out] WORK 69*> \verbatim 70*> WORK is COMPLEX array, dimension (LDWORK,N) 71*> \endverbatim 72*> 73*> \param[in] LDWORK 74*> \verbatim 75*> LDWORK is INTEGER 76*> The leading dimension of the array WORK. LDWORK >= max(1,N). 77*> \endverbatim 78*> 79*> \param[out] RWORK 80*> \verbatim 81*> RWORK is REAL array, dimension (N) 82*> \endverbatim 83*> 84*> \param[out] RCOND 85*> \verbatim 86*> RCOND is REAL 87*> The reciprocal of the condition number of A, computed as 88*> ( 1/norm(A) ) / norm(AINV). 89*> \endverbatim 90*> 91*> \param[out] RESID 92*> \verbatim 93*> RESID is REAL 94*> norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS ) 95*> \endverbatim 96* 97* Authors: 98* ======== 99* 100*> \author Univ. of Tennessee 101*> \author Univ. of California Berkeley 102*> \author Univ. of Colorado Denver 103*> \author NAG Ltd. 104* 105*> \date December 2016 106* 107*> \ingroup complex_lin 108* 109* ===================================================================== 110 SUBROUTINE CPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND, 111 $ RESID ) 112* 113* -- LAPACK test routine (version 3.7.0) -- 114* -- LAPACK is a software package provided by Univ. of Tennessee, -- 115* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 116* December 2016 117* 118* .. Scalar Arguments .. 119 CHARACTER UPLO 120 INTEGER LDWORK, N 121 REAL RCOND, RESID 122* .. 123* .. Array Arguments .. 124 REAL RWORK( * ) 125 COMPLEX A( * ), AINV( * ), WORK( LDWORK, * ) 126* .. 127* 128* ===================================================================== 129* 130* .. Parameters .. 131 REAL ZERO, ONE 132 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 133 COMPLEX CZERO, CONE 134 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 135 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 136* .. 137* .. Local Scalars .. 138 INTEGER I, J, JJ 139 REAL AINVNM, ANORM, EPS 140* .. 141* .. External Functions .. 142 LOGICAL LSAME 143 REAL CLANGE, CLANHP, SLAMCH 144 EXTERNAL LSAME, CLANGE, CLANHP, SLAMCH 145* .. 146* .. Intrinsic Functions .. 147 INTRINSIC CONJG, REAL 148* .. 149* .. External Subroutines .. 150 EXTERNAL CCOPY, CHPMV 151* .. 152* .. Executable Statements .. 153* 154* Quick exit if N = 0. 155* 156 IF( N.LE.0 ) THEN 157 RCOND = ONE 158 RESID = ZERO 159 RETURN 160 END IF 161* 162* Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 163* 164 EPS = SLAMCH( 'Epsilon' ) 165 ANORM = CLANHP( '1', UPLO, N, A, RWORK ) 166 AINVNM = CLANHP( '1', UPLO, N, AINV, RWORK ) 167 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 168 RCOND = ZERO 169 RESID = ONE / EPS 170 RETURN 171 END IF 172 RCOND = ( ONE/ANORM ) / AINVNM 173* 174* UPLO = 'U': 175* Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and 176* expand it to a full matrix, then multiply by A one column at a 177* time, moving the result one column to the left. 178* 179 IF( LSAME( UPLO, 'U' ) ) THEN 180* 181* Copy AINV 182* 183 JJ = 1 184 DO 20 J = 1, N - 1 185 CALL CCOPY( J, AINV( JJ ), 1, WORK( 1, J+1 ), 1 ) 186 DO 10 I = 1, J - 1 187 WORK( J, I+1 ) = CONJG( AINV( JJ+I-1 ) ) 188 10 CONTINUE 189 JJ = JJ + J 190 20 CONTINUE 191 JJ = ( ( N-1 )*N ) / 2 + 1 192 DO 30 I = 1, N - 1 193 WORK( N, I+1 ) = CONJG( AINV( JJ+I-1 ) ) 194 30 CONTINUE 195* 196* Multiply by A 197* 198 DO 40 J = 1, N - 1 199 CALL CHPMV( 'Upper', N, -CONE, A, WORK( 1, J+1 ), 1, CZERO, 200 $ WORK( 1, J ), 1 ) 201 40 CONTINUE 202 CALL CHPMV( 'Upper', N, -CONE, A, AINV( JJ ), 1, CZERO, 203 $ WORK( 1, N ), 1 ) 204* 205* UPLO = 'L': 206* Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1) 207* and multiply by A, moving each column to the right. 208* 209 ELSE 210* 211* Copy AINV 212* 213 DO 50 I = 1, N - 1 214 WORK( 1, I ) = CONJG( AINV( I+1 ) ) 215 50 CONTINUE 216 JJ = N + 1 217 DO 70 J = 2, N 218 CALL CCOPY( N-J+1, AINV( JJ ), 1, WORK( J, J-1 ), 1 ) 219 DO 60 I = 1, N - J 220 WORK( J, J+I-1 ) = CONJG( AINV( JJ+I ) ) 221 60 CONTINUE 222 JJ = JJ + N - J + 1 223 70 CONTINUE 224* 225* Multiply by A 226* 227 DO 80 J = N, 2, -1 228 CALL CHPMV( 'Lower', N, -CONE, A, WORK( 1, J-1 ), 1, CZERO, 229 $ WORK( 1, J ), 1 ) 230 80 CONTINUE 231 CALL CHPMV( 'Lower', N, -CONE, A, AINV( 1 ), 1, CZERO, 232 $ WORK( 1, 1 ), 1 ) 233* 234 END IF 235* 236* Add the identity matrix to WORK . 237* 238 DO 90 I = 1, N 239 WORK( I, I ) = WORK( I, I ) + CONE 240 90 CONTINUE 241* 242* Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 243* 244 RESID = CLANGE( '1', N, N, WORK, LDWORK, RWORK ) 245* 246 RESID = ( ( RESID*RCOND )/EPS ) / REAL( N ) 247* 248 RETURN 249* 250* End of CPPT03 251* 252 END 253