1*> \brief \b DPPT03 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 DPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND, 12* RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER UPLO 16* INTEGER LDWORK, N 17* DOUBLE PRECISION RCOND, RESID 18* .. 19* .. Array Arguments .. 20* DOUBLE PRECISION A( * ), AINV( * ), RWORK( * ), 21* $ WORK( LDWORK, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> DPPT03 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N) 82*> \endverbatim 83*> 84*> \param[out] RCOND 85*> \verbatim 86*> RCOND is DOUBLE PRECISION 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 DOUBLE PRECISION 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 double_lin 108* 109* ===================================================================== 110 SUBROUTINE DPPT03( 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 DOUBLE PRECISION RCOND, RESID 122* .. 123* .. Array Arguments .. 124 DOUBLE PRECISION A( * ), AINV( * ), RWORK( * ), 125 $ WORK( LDWORK, * ) 126* .. 127* 128* ===================================================================== 129* 130* .. Parameters .. 131 DOUBLE PRECISION ZERO, ONE 132 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 133* .. 134* .. Local Scalars .. 135 INTEGER I, J, JJ 136 DOUBLE PRECISION AINVNM, ANORM, EPS 137* .. 138* .. External Functions .. 139 LOGICAL LSAME 140 DOUBLE PRECISION DLAMCH, DLANGE, DLANSP 141 EXTERNAL LSAME, DLAMCH, DLANGE, DLANSP 142* .. 143* .. Intrinsic Functions .. 144 INTRINSIC DBLE 145* .. 146* .. External Subroutines .. 147 EXTERNAL DCOPY, DSPMV 148* .. 149* .. Executable Statements .. 150* 151* Quick exit if N = 0. 152* 153 IF( N.LE.0 ) THEN 154 RCOND = ONE 155 RESID = ZERO 156 RETURN 157 END IF 158* 159* Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 160* 161 EPS = DLAMCH( 'Epsilon' ) 162 ANORM = DLANSP( '1', UPLO, N, A, RWORK ) 163 AINVNM = DLANSP( '1', UPLO, N, AINV, RWORK ) 164 IF( ANORM.LE.ZERO .OR. AINVNM.EQ.ZERO ) THEN 165 RCOND = ZERO 166 RESID = ONE / EPS 167 RETURN 168 END IF 169 RCOND = ( ONE / ANORM ) / AINVNM 170* 171* UPLO = 'U': 172* Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and 173* expand it to a full matrix, then multiply by A one column at a 174* time, moving the result one column to the left. 175* 176 IF( LSAME( UPLO, 'U' ) ) THEN 177* 178* Copy AINV 179* 180 JJ = 1 181 DO 10 J = 1, N - 1 182 CALL DCOPY( J, AINV( JJ ), 1, WORK( 1, J+1 ), 1 ) 183 CALL DCOPY( J-1, AINV( JJ ), 1, WORK( J, 2 ), LDWORK ) 184 JJ = JJ + J 185 10 CONTINUE 186 JJ = ( ( N-1 )*N ) / 2 + 1 187 CALL DCOPY( N-1, AINV( JJ ), 1, WORK( N, 2 ), LDWORK ) 188* 189* Multiply by A 190* 191 DO 20 J = 1, N - 1 192 CALL DSPMV( 'Upper', N, -ONE, A, WORK( 1, J+1 ), 1, ZERO, 193 $ WORK( 1, J ), 1 ) 194 20 CONTINUE 195 CALL DSPMV( 'Upper', N, -ONE, A, AINV( JJ ), 1, ZERO, 196 $ WORK( 1, N ), 1 ) 197* 198* UPLO = 'L': 199* Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1) 200* and multiply by A, moving each column to the right. 201* 202 ELSE 203* 204* Copy AINV 205* 206 CALL DCOPY( N-1, AINV( 2 ), 1, WORK( 1, 1 ), LDWORK ) 207 JJ = N + 1 208 DO 30 J = 2, N 209 CALL DCOPY( N-J+1, AINV( JJ ), 1, WORK( J, J-1 ), 1 ) 210 CALL DCOPY( N-J, AINV( JJ+1 ), 1, WORK( J, J ), LDWORK ) 211 JJ = JJ + N - J + 1 212 30 CONTINUE 213* 214* Multiply by A 215* 216 DO 40 J = N, 2, -1 217 CALL DSPMV( 'Lower', N, -ONE, A, WORK( 1, J-1 ), 1, ZERO, 218 $ WORK( 1, J ), 1 ) 219 40 CONTINUE 220 CALL DSPMV( 'Lower', N, -ONE, A, AINV( 1 ), 1, ZERO, 221 $ WORK( 1, 1 ), 1 ) 222* 223 END IF 224* 225* Add the identity matrix to WORK . 226* 227 DO 50 I = 1, N 228 WORK( I, I ) = WORK( I, I ) + ONE 229 50 CONTINUE 230* 231* Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 232* 233 RESID = DLANGE( '1', N, N, WORK, LDWORK, RWORK ) 234* 235 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N ) 236* 237 RETURN 238* 239* End of DPPT03 240* 241 END 242