1*> \brief \b CSPT03 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 CSPT03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND, 12* RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER UPLO 16* INTEGER LDW, N 17* REAL RCOND, RESID 18* .. 19* .. Array Arguments .. 20* REAL RWORK( * ) 21* COMPLEX A( * ), AINV( * ), WORK( LDW, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> CSPT03 computes the residual for a complex symmetric packed matrix 31*> times its 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*> complex symmetric 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 complex symmetric 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 (symmetric) 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 (LDW,N) 71*> \endverbatim 72*> 73*> \param[in] LDW 74*> \verbatim 75*> LDW is INTEGER 76*> The leading dimension of the array WORK. LDW >= 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 November 2011 106* 107*> \ingroup complex_lin 108* 109* ===================================================================== 110 SUBROUTINE CSPT03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND, 111 $ RESID ) 112* 113* -- LAPACK test routine (version 3.4.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* November 2011 117* 118* .. Scalar Arguments .. 119 CHARACTER UPLO 120 INTEGER LDW, N 121 REAL RCOND, RESID 122* .. 123* .. Array Arguments .. 124 REAL RWORK( * ) 125 COMPLEX A( * ), AINV( * ), WORK( LDW, * ) 126* .. 127* 128* ===================================================================== 129* 130* .. Parameters .. 131 REAL ZERO, ONE 132 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 133* .. 134* .. Local Scalars .. 135 INTEGER I, ICOL, J, JCOL, K, KCOL, NALL 136 REAL AINVNM, ANORM, EPS 137 COMPLEX T 138* .. 139* .. External Functions .. 140 LOGICAL LSAME 141 REAL CLANGE, CLANSP, SLAMCH 142 COMPLEX CDOTU 143 EXTERNAL LSAME, CLANGE, CLANSP, SLAMCH, CDOTU 144* .. 145* .. Intrinsic Functions .. 146 INTRINSIC REAL 147* .. 148* .. Executable Statements .. 149* 150* Quick exit if N = 0. 151* 152 IF( N.LE.0 ) THEN 153 RCOND = ONE 154 RESID = ZERO 155 RETURN 156 END IF 157* 158* Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 159* 160 EPS = SLAMCH( 'Epsilon' ) 161 ANORM = CLANSP( '1', UPLO, N, A, RWORK ) 162 AINVNM = CLANSP( '1', UPLO, N, AINV, RWORK ) 163 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 164 RCOND = ZERO 165 RESID = ONE / EPS 166 RETURN 167 END IF 168 RCOND = ( ONE/ANORM ) / AINVNM 169* 170* Case where both A and AINV are upper triangular: 171* Each element of - A * AINV is computed by taking the dot product 172* of a row of A with a column of AINV. 173* 174 IF( LSAME( UPLO, 'U' ) ) THEN 175 DO 70 I = 1, N 176 ICOL = ( ( I-1 )*I ) / 2 + 1 177* 178* Code when J <= I 179* 180 DO 30 J = 1, I 181 JCOL = ( ( J-1 )*J ) / 2 + 1 182 T = CDOTU( J, A( ICOL ), 1, AINV( JCOL ), 1 ) 183 JCOL = JCOL + 2*J - 1 184 KCOL = ICOL - 1 185 DO 10 K = J + 1, I 186 T = T + A( KCOL+K )*AINV( JCOL ) 187 JCOL = JCOL + K 188 10 CONTINUE 189 KCOL = KCOL + 2*I 190 DO 20 K = I + 1, N 191 T = T + A( KCOL )*AINV( JCOL ) 192 KCOL = KCOL + K 193 JCOL = JCOL + K 194 20 CONTINUE 195 WORK( I, J ) = -T 196 30 CONTINUE 197* 198* Code when J > I 199* 200 DO 60 J = I + 1, N 201 JCOL = ( ( J-1 )*J ) / 2 + 1 202 T = CDOTU( I, A( ICOL ), 1, AINV( JCOL ), 1 ) 203 JCOL = JCOL - 1 204 KCOL = ICOL + 2*I - 1 205 DO 40 K = I + 1, J 206 T = T + A( KCOL )*AINV( JCOL+K ) 207 KCOL = KCOL + K 208 40 CONTINUE 209 JCOL = JCOL + 2*J 210 DO 50 K = J + 1, N 211 T = T + A( KCOL )*AINV( JCOL ) 212 KCOL = KCOL + K 213 JCOL = JCOL + K 214 50 CONTINUE 215 WORK( I, J ) = -T 216 60 CONTINUE 217 70 CONTINUE 218 ELSE 219* 220* Case where both A and AINV are lower triangular 221* 222 NALL = ( N*( N+1 ) ) / 2 223 DO 140 I = 1, N 224* 225* Code when J <= I 226* 227 ICOL = NALL - ( ( N-I+1 )*( N-I+2 ) ) / 2 + 1 228 DO 100 J = 1, I 229 JCOL = NALL - ( ( N-J )*( N-J+1 ) ) / 2 - ( N-I ) 230 T = CDOTU( N-I+1, A( ICOL ), 1, AINV( JCOL ), 1 ) 231 KCOL = I 232 JCOL = J 233 DO 80 K = 1, J - 1 234 T = T + A( KCOL )*AINV( JCOL ) 235 JCOL = JCOL + N - K 236 KCOL = KCOL + N - K 237 80 CONTINUE 238 JCOL = JCOL - J 239 DO 90 K = J, I - 1 240 T = T + A( KCOL )*AINV( JCOL+K ) 241 KCOL = KCOL + N - K 242 90 CONTINUE 243 WORK( I, J ) = -T 244 100 CONTINUE 245* 246* Code when J > I 247* 248 ICOL = NALL - ( ( N-I )*( N-I+1 ) ) / 2 249 DO 130 J = I + 1, N 250 JCOL = NALL - ( ( N-J+1 )*( N-J+2 ) ) / 2 + 1 251 T = CDOTU( N-J+1, A( ICOL-N+J ), 1, AINV( JCOL ), 1 ) 252 KCOL = I 253 JCOL = J 254 DO 110 K = 1, I - 1 255 T = T + A( KCOL )*AINV( JCOL ) 256 JCOL = JCOL + N - K 257 KCOL = KCOL + N - K 258 110 CONTINUE 259 KCOL = KCOL - I 260 DO 120 K = I, J - 1 261 T = T + A( KCOL+K )*AINV( JCOL ) 262 JCOL = JCOL + N - K 263 120 CONTINUE 264 WORK( I, J ) = -T 265 130 CONTINUE 266 140 CONTINUE 267 END IF 268* 269* Add the identity matrix to WORK . 270* 271 DO 150 I = 1, N 272 WORK( I, I ) = WORK( I, I ) + ONE 273 150 CONTINUE 274* 275* Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 276* 277 RESID = CLANGE( '1', N, N, WORK, LDW, RWORK ) 278* 279 RESID = ( ( RESID*RCOND )/EPS ) / REAL( N ) 280* 281 RETURN 282* 283* End of CSPT03 284* 285 END 286