1*> \brief \b SPPT03 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 SPPT03( 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 A( * ), AINV( * ), RWORK( * ), 21* $ WORK( LDWORK, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> SPPT03 computes the residual for a symmetric 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*> 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 REAL array, dimension (N*(N+1)/2) 57*> The original symmetric matrix A, stored as a packed 58*> triangular matrix. 59*> \endverbatim 60*> 61*> \param[in] AINV 62*> \verbatim 63*> AINV is REAL 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 REAL 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*> \ingroup single_lin 106* 107* ===================================================================== 108 SUBROUTINE SPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND, 109 $ RESID ) 110* 111* -- LAPACK test routine -- 112* -- LAPACK is a software package provided by Univ. of Tennessee, -- 113* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 114* 115* .. Scalar Arguments .. 116 CHARACTER UPLO 117 INTEGER LDWORK, N 118 REAL RCOND, RESID 119* .. 120* .. Array Arguments .. 121 REAL A( * ), AINV( * ), RWORK( * ), 122 $ WORK( LDWORK, * ) 123* .. 124* 125* ===================================================================== 126* 127* .. Parameters .. 128 REAL ZERO, ONE 129 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 130* .. 131* .. Local Scalars .. 132 INTEGER I, J, JJ 133 REAL AINVNM, ANORM, EPS 134* .. 135* .. External Functions .. 136 LOGICAL LSAME 137 REAL SLAMCH, SLANGE, SLANSP 138 EXTERNAL LSAME, SLAMCH, SLANGE, SLANSP 139* .. 140* .. Intrinsic Functions .. 141 INTRINSIC REAL 142* .. 143* .. External Subroutines .. 144 EXTERNAL SCOPY, SSPMV 145* .. 146* .. Executable Statements .. 147* 148* Quick exit if N = 0. 149* 150 IF( N.LE.0 ) THEN 151 RCOND = ONE 152 RESID = ZERO 153 RETURN 154 END IF 155* 156* Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 157* 158 EPS = SLAMCH( 'Epsilon' ) 159 ANORM = SLANSP( '1', UPLO, N, A, RWORK ) 160 AINVNM = SLANSP( '1', UPLO, N, AINV, RWORK ) 161 IF( ANORM.LE.ZERO .OR. AINVNM.EQ.ZERO ) THEN 162 RCOND = ZERO 163 RESID = ONE / EPS 164 RETURN 165 END IF 166 RCOND = ( ONE / ANORM ) / AINVNM 167* 168* UPLO = 'U': 169* Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and 170* expand it to a full matrix, then multiply by A one column at a 171* time, moving the result one column to the left. 172* 173 IF( LSAME( UPLO, 'U' ) ) THEN 174* 175* Copy AINV 176* 177 JJ = 1 178 DO 10 J = 1, N - 1 179 CALL SCOPY( J, AINV( JJ ), 1, WORK( 1, J+1 ), 1 ) 180 CALL SCOPY( J-1, AINV( JJ ), 1, WORK( J, 2 ), LDWORK ) 181 JJ = JJ + J 182 10 CONTINUE 183 JJ = ( ( N-1 )*N ) / 2 + 1 184 CALL SCOPY( N-1, AINV( JJ ), 1, WORK( N, 2 ), LDWORK ) 185* 186* Multiply by A 187* 188 DO 20 J = 1, N - 1 189 CALL SSPMV( 'Upper', N, -ONE, A, WORK( 1, J+1 ), 1, ZERO, 190 $ WORK( 1, J ), 1 ) 191 20 CONTINUE 192 CALL SSPMV( 'Upper', N, -ONE, A, AINV( JJ ), 1, ZERO, 193 $ WORK( 1, N ), 1 ) 194* 195* UPLO = 'L': 196* Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1) 197* and multiply by A, moving each column to the right. 198* 199 ELSE 200* 201* Copy AINV 202* 203 CALL SCOPY( N-1, AINV( 2 ), 1, WORK( 1, 1 ), LDWORK ) 204 JJ = N + 1 205 DO 30 J = 2, N 206 CALL SCOPY( N-J+1, AINV( JJ ), 1, WORK( J, J-1 ), 1 ) 207 CALL SCOPY( N-J, AINV( JJ+1 ), 1, WORK( J, J ), LDWORK ) 208 JJ = JJ + N - J + 1 209 30 CONTINUE 210* 211* Multiply by A 212* 213 DO 40 J = N, 2, -1 214 CALL SSPMV( 'Lower', N, -ONE, A, WORK( 1, J-1 ), 1, ZERO, 215 $ WORK( 1, J ), 1 ) 216 40 CONTINUE 217 CALL SSPMV( 'Lower', N, -ONE, A, AINV( 1 ), 1, ZERO, 218 $ WORK( 1, 1 ), 1 ) 219* 220 END IF 221* 222* Add the identity matrix to WORK . 223* 224 DO 50 I = 1, N 225 WORK( I, I ) = WORK( I, I ) + ONE 226 50 CONTINUE 227* 228* Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 229* 230 RESID = SLANGE( '1', N, N, WORK, LDWORK, RWORK ) 231* 232 RESID = ( ( RESID*RCOND ) / EPS ) / REAL( N ) 233* 234 RETURN 235* 236* End of SPPT03 237* 238 END 239