1*> \brief \b SQRT16 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 SQRT16( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, 12* RWORK, RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER TRANS 16* INTEGER LDA, LDB, LDX, M, N, NRHS 17* REAL RESID 18* .. 19* .. Array Arguments .. 20* REAL A( LDA, * ), B( LDB, * ), RWORK( * ), 21* $ X( LDX, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> SQRT16 computes the residual for a solution of a system of linear 31*> equations A*x = b or A'*x = b: 32*> RESID = norm(B - A*X) / ( max(m,n) * norm(A) * norm(X) * EPS ), 33*> where EPS is the machine epsilon. 34*> \endverbatim 35* 36* Arguments: 37* ========== 38* 39*> \param[in] TRANS 40*> \verbatim 41*> TRANS is CHARACTER*1 42*> Specifies the form of the system of equations: 43*> = 'N': A *x = b 44*> = 'T': A'*x = b, where A' is the transpose of A 45*> = 'C': A'*x = b, where A' is the transpose of A 46*> \endverbatim 47*> 48*> \param[in] M 49*> \verbatim 50*> M is INTEGER 51*> The number of rows of the matrix A. M >= 0. 52*> \endverbatim 53*> 54*> \param[in] N 55*> \verbatim 56*> N is INTEGER 57*> The number of columns of the matrix A. N >= 0. 58*> \endverbatim 59*> 60*> \param[in] NRHS 61*> \verbatim 62*> NRHS is INTEGER 63*> The number of columns of B, the matrix of right hand sides. 64*> NRHS >= 0. 65*> \endverbatim 66*> 67*> \param[in] A 68*> \verbatim 69*> A is REAL array, dimension (LDA,N) 70*> The original M x N matrix A. 71*> \endverbatim 72*> 73*> \param[in] LDA 74*> \verbatim 75*> LDA is INTEGER 76*> The leading dimension of the array A. LDA >= max(1,M). 77*> \endverbatim 78*> 79*> \param[in] X 80*> \verbatim 81*> X is REAL array, dimension (LDX,NRHS) 82*> The computed solution vectors for the system of linear 83*> equations. 84*> \endverbatim 85*> 86*> \param[in] LDX 87*> \verbatim 88*> LDX is INTEGER 89*> The leading dimension of the array X. If TRANS = 'N', 90*> LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M). 91*> \endverbatim 92*> 93*> \param[in,out] B 94*> \verbatim 95*> B is REAL array, dimension (LDB,NRHS) 96*> On entry, the right hand side vectors for the system of 97*> linear equations. 98*> On exit, B is overwritten with the difference B - A*X. 99*> \endverbatim 100*> 101*> \param[in] LDB 102*> \verbatim 103*> LDB is INTEGER 104*> The leading dimension of the array B. IF TRANS = 'N', 105*> LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N). 106*> \endverbatim 107*> 108*> \param[out] RWORK 109*> \verbatim 110*> RWORK is REAL array, dimension (M) 111*> \endverbatim 112*> 113*> \param[out] RESID 114*> \verbatim 115*> RESID is REAL 116*> The maximum over the number of right hand sides of 117*> norm(B - A*X) / ( max(m,n) * norm(A) * norm(X) * EPS ). 118*> \endverbatim 119* 120* Authors: 121* ======== 122* 123*> \author Univ. of Tennessee 124*> \author Univ. of California Berkeley 125*> \author Univ. of Colorado Denver 126*> \author NAG Ltd. 127* 128*> \ingroup single_lin 129* 130* ===================================================================== 131 SUBROUTINE SQRT16( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, 132 $ RWORK, RESID ) 133* 134* -- LAPACK test routine -- 135* -- LAPACK is a software package provided by Univ. of Tennessee, -- 136* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 137* 138* .. Scalar Arguments .. 139 CHARACTER TRANS 140 INTEGER LDA, LDB, LDX, M, N, NRHS 141 REAL RESID 142* .. 143* .. Array Arguments .. 144 REAL A( LDA, * ), B( LDB, * ), RWORK( * ), 145 $ X( LDX, * ) 146* .. 147* 148* ===================================================================== 149* 150* .. Parameters .. 151 REAL ZERO, ONE 152 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 153* .. 154* .. Local Scalars .. 155 INTEGER J, N1, N2 156 REAL ANORM, BNORM, EPS, XNORM 157* .. 158* .. External Functions .. 159 LOGICAL LSAME 160 REAL SASUM, SLAMCH, SLANGE 161 EXTERNAL LSAME, SASUM, SLAMCH, SLANGE 162* .. 163* .. External Subroutines .. 164 EXTERNAL SGEMM 165* .. 166* .. Intrinsic Functions .. 167 INTRINSIC MAX 168* .. 169* .. Executable Statements .. 170* 171* Quick exit if M = 0 or N = 0 or NRHS = 0 172* 173 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN 174 RESID = ZERO 175 RETURN 176 END IF 177* 178 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN 179 ANORM = SLANGE( 'I', M, N, A, LDA, RWORK ) 180 N1 = N 181 N2 = M 182 ELSE 183 ANORM = SLANGE( '1', M, N, A, LDA, RWORK ) 184 N1 = M 185 N2 = N 186 END IF 187* 188 EPS = SLAMCH( 'Epsilon' ) 189* 190* Compute B - A*X (or B - A'*X ) and store in B. 191* 192 CALL SGEMM( TRANS, 'No transpose', N1, NRHS, N2, -ONE, A, LDA, X, 193 $ LDX, ONE, B, LDB ) 194* 195* Compute the maximum over the number of right hand sides of 196* norm(B - A*X) / ( max(m,n) * norm(A) * norm(X) * EPS ) . 197* 198 RESID = ZERO 199 DO 10 J = 1, NRHS 200 BNORM = SASUM( N1, B( 1, J ), 1 ) 201 XNORM = SASUM( N2, X( 1, J ), 1 ) 202 IF( ANORM.EQ.ZERO .AND. BNORM.EQ.ZERO ) THEN 203 RESID = ZERO 204 ELSE IF( ANORM.LE.ZERO .OR. XNORM.LE.ZERO ) THEN 205 RESID = ONE / EPS 206 ELSE 207 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / 208 $ ( MAX( M, N )*EPS ) ) 209 END IF 210 10 CONTINUE 211* 212 RETURN 213* 214* End of SQRT16 215* 216 END 217