1*> \brief \b DPOT03 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 DPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK, 12* RWORK, RCOND, RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER UPLO 16* INTEGER LDA, LDAINV, LDWORK, N 17* DOUBLE PRECISION RCOND, RESID 18* .. 19* .. Array Arguments .. 20* DOUBLE PRECISION A( LDA, * ), AINV( LDAINV, * ), RWORK( * ), 21* $ WORK( LDWORK, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> DPOT03 computes the residual for a symmetric 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 (LDA,N) 57*> The original symmetric matrix A. 58*> \endverbatim 59*> 60*> \param[in] LDA 61*> \verbatim 62*> LDA is INTEGER 63*> The leading dimension of the array A. LDA >= max(1,N) 64*> \endverbatim 65*> 66*> \param[in,out] AINV 67*> \verbatim 68*> AINV is DOUBLE PRECISION array, dimension (LDAINV,N) 69*> On entry, the inverse of the matrix A, stored as a symmetric 70*> matrix in the same format as A. 71*> In this version, AINV is expanded into a full matrix and 72*> multiplied by A, so the opposing triangle of AINV will be 73*> changed; i.e., if the upper triangular part of AINV is 74*> stored, the lower triangular part will be used as work space. 75*> \endverbatim 76*> 77*> \param[in] LDAINV 78*> \verbatim 79*> LDAINV is INTEGER 80*> The leading dimension of the array AINV. LDAINV >= max(1,N). 81*> \endverbatim 82*> 83*> \param[out] WORK 84*> \verbatim 85*> WORK is DOUBLE PRECISION array, dimension (LDWORK,N) 86*> \endverbatim 87*> 88*> \param[in] LDWORK 89*> \verbatim 90*> LDWORK is INTEGER 91*> The leading dimension of the array WORK. LDWORK >= max(1,N). 92*> \endverbatim 93*> 94*> \param[out] RWORK 95*> \verbatim 96*> RWORK is DOUBLE PRECISION array, dimension (N) 97*> \endverbatim 98*> 99*> \param[out] RCOND 100*> \verbatim 101*> RCOND is DOUBLE PRECISION 102*> The reciprocal of the condition number of A, computed as 103*> ( 1/norm(A) ) / norm(AINV). 104*> \endverbatim 105*> 106*> \param[out] RESID 107*> \verbatim 108*> RESID is DOUBLE PRECISION 109*> norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS ) 110*> \endverbatim 111* 112* Authors: 113* ======== 114* 115*> \author Univ. of Tennessee 116*> \author Univ. of California Berkeley 117*> \author Univ. of Colorado Denver 118*> \author NAG Ltd. 119* 120*> \ingroup double_lin 121* 122* ===================================================================== 123 SUBROUTINE DPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK, 124 $ RWORK, RCOND, RESID ) 125* 126* -- LAPACK test routine -- 127* -- LAPACK is a software package provided by Univ. of Tennessee, -- 128* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 129* 130* .. Scalar Arguments .. 131 CHARACTER UPLO 132 INTEGER LDA, LDAINV, LDWORK, N 133 DOUBLE PRECISION RCOND, RESID 134* .. 135* .. Array Arguments .. 136 DOUBLE PRECISION A( LDA, * ), AINV( LDAINV, * ), RWORK( * ), 137 $ WORK( LDWORK, * ) 138* .. 139* 140* ===================================================================== 141* 142* .. Parameters .. 143 DOUBLE PRECISION ZERO, ONE 144 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 145* .. 146* .. Local Scalars .. 147 INTEGER I, J 148 DOUBLE PRECISION AINVNM, ANORM, EPS 149* .. 150* .. External Functions .. 151 LOGICAL LSAME 152 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 153 EXTERNAL LSAME, DLAMCH, DLANGE, DLANSY 154* .. 155* .. External Subroutines .. 156 EXTERNAL DSYMM 157* .. 158* .. Intrinsic Functions .. 159 INTRINSIC DBLE 160* .. 161* .. Executable Statements .. 162* 163* Quick exit if N = 0. 164* 165 IF( N.LE.0 ) THEN 166 RCOND = ONE 167 RESID = ZERO 168 RETURN 169 END IF 170* 171* Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 172* 173 EPS = DLAMCH( 'Epsilon' ) 174 ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK ) 175 AINVNM = DLANSY( '1', UPLO, N, AINV, LDAINV, RWORK ) 176 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 177 RCOND = ZERO 178 RESID = ONE / EPS 179 RETURN 180 END IF 181 RCOND = ( ONE / ANORM ) / AINVNM 182* 183* Expand AINV into a full matrix and call DSYMM to multiply 184* AINV on the left by A. 185* 186 IF( LSAME( UPLO, 'U' ) ) THEN 187 DO 20 J = 1, N 188 DO 10 I = 1, J - 1 189 AINV( J, I ) = AINV( I, J ) 190 10 CONTINUE 191 20 CONTINUE 192 ELSE 193 DO 40 J = 1, N 194 DO 30 I = J + 1, N 195 AINV( J, I ) = AINV( I, J ) 196 30 CONTINUE 197 40 CONTINUE 198 END IF 199 CALL DSYMM( 'Left', UPLO, N, N, -ONE, A, LDA, AINV, LDAINV, ZERO, 200 $ WORK, LDWORK ) 201* 202* Add the identity matrix to WORK . 203* 204 DO 50 I = 1, N 205 WORK( I, I ) = WORK( I, I ) + ONE 206 50 CONTINUE 207* 208* Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 209* 210 RESID = DLANGE( '1', N, N, WORK, LDWORK, RWORK ) 211* 212 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N ) 213* 214 RETURN 215* 216* End of DPOT03 217* 218 END 219