![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DQRT17 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * DOUBLE PRECISION FUNCTION DQRT17( TRANS, IRESID, M, N, NRHS, A, 00012 * LDA, X, LDX, B, LDB, C, WORK, LWORK ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER TRANS 00016 * INTEGER IRESID, LDA, LDB, LDX, LWORK, M, N, NRHS 00017 * .. 00018 * .. Array Arguments .. 00019 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDB, * ), 00020 * $ WORK( LWORK ), X( LDX, * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> DQRT17 computes the ratio 00030 *> 00031 *> || R'*op(A) ||/(||A||*alpha*max(M,N,NRHS)*eps) 00032 *> 00033 *> where R = op(A)*X - B, op(A) is A or A', and 00034 *> 00035 *> alpha = ||B|| if IRESID = 1 (zero-residual problem) 00036 *> alpha = ||R|| if IRESID = 2 (otherwise). 00037 *> \endverbatim 00038 * 00039 * Arguments: 00040 * ========== 00041 * 00042 *> \param[in] TRANS 00043 *> \verbatim 00044 *> TRANS is CHARACTER*1 00045 *> Specifies whether or not the transpose of A is used. 00046 *> = 'N': No transpose, op(A) = A. 00047 *> = 'T': Transpose, op(A) = A'. 00048 *> \endverbatim 00049 *> 00050 *> \param[in] IRESID 00051 *> \verbatim 00052 *> IRESID is INTEGER 00053 *> IRESID = 1 indicates zero-residual problem. 00054 *> IRESID = 2 indicates non-zero residual. 00055 *> \endverbatim 00056 *> 00057 *> \param[in] M 00058 *> \verbatim 00059 *> M is INTEGER 00060 *> The number of rows of the matrix A. 00061 *> If TRANS = 'N', the number of rows of the matrix B. 00062 *> If TRANS = 'T', the number of rows of the matrix X. 00063 *> \endverbatim 00064 *> 00065 *> \param[in] N 00066 *> \verbatim 00067 *> N is INTEGER 00068 *> The number of columns of the matrix A. 00069 *> If TRANS = 'N', the number of rows of the matrix X. 00070 *> If TRANS = 'T', the number of rows of the matrix B. 00071 *> \endverbatim 00072 *> 00073 *> \param[in] NRHS 00074 *> \verbatim 00075 *> NRHS is INTEGER 00076 *> The number of columns of the matrices X and B. 00077 *> \endverbatim 00078 *> 00079 *> \param[in] A 00080 *> \verbatim 00081 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00082 *> The m-by-n matrix A. 00083 *> \endverbatim 00084 *> 00085 *> \param[in] LDA 00086 *> \verbatim 00087 *> LDA is INTEGER 00088 *> The leading dimension of the array A. LDA >= M. 00089 *> \endverbatim 00090 *> 00091 *> \param[in] X 00092 *> \verbatim 00093 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 00094 *> If TRANS = 'N', the n-by-nrhs matrix X. 00095 *> If TRANS = 'T', the m-by-nrhs matrix X. 00096 *> \endverbatim 00097 *> 00098 *> \param[in] LDX 00099 *> \verbatim 00100 *> LDX is INTEGER 00101 *> The leading dimension of the array X. 00102 *> If TRANS = 'N', LDX >= N. 00103 *> If TRANS = 'T', LDX >= M. 00104 *> \endverbatim 00105 *> 00106 *> \param[in] B 00107 *> \verbatim 00108 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 00109 *> If TRANS = 'N', the m-by-nrhs matrix B. 00110 *> If TRANS = 'T', the n-by-nrhs matrix B. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] LDB 00114 *> \verbatim 00115 *> LDB is INTEGER 00116 *> The leading dimension of the array B. 00117 *> If TRANS = 'N', LDB >= M. 00118 *> If TRANS = 'T', LDB >= N. 00119 *> \endverbatim 00120 *> 00121 *> \param[out] C 00122 *> \verbatim 00123 *> C is DOUBLE PRECISION array, dimension (LDB,NRHS) 00124 *> \endverbatim 00125 *> 00126 *> \param[out] WORK 00127 *> \verbatim 00128 *> WORK is DOUBLE PRECISION array, dimension (LWORK) 00129 *> \endverbatim 00130 *> 00131 *> \param[in] LWORK 00132 *> \verbatim 00133 *> LWORK is INTEGER 00134 *> The length of the array WORK. LWORK >= NRHS*(M+N). 00135 *> \endverbatim 00136 * 00137 * Authors: 00138 * ======== 00139 * 00140 *> \author Univ. of Tennessee 00141 *> \author Univ. of California Berkeley 00142 *> \author Univ. of Colorado Denver 00143 *> \author NAG Ltd. 00144 * 00145 *> \date November 2011 00146 * 00147 *> \ingroup double_lin 00148 * 00149 * ===================================================================== 00150 DOUBLE PRECISION FUNCTION DQRT17( TRANS, IRESID, M, N, NRHS, A, 00151 $ LDA, X, LDX, B, LDB, C, WORK, LWORK ) 00152 * 00153 * -- LAPACK test routine (version 3.4.0) -- 00154 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00155 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00156 * November 2011 00157 * 00158 * .. Scalar Arguments .. 00159 CHARACTER TRANS 00160 INTEGER IRESID, LDA, LDB, LDX, LWORK, M, N, NRHS 00161 * .. 00162 * .. Array Arguments .. 00163 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDB, * ), 00164 $ WORK( LWORK ), X( LDX, * ) 00165 * .. 00166 * 00167 * ===================================================================== 00168 * 00169 * .. Parameters .. 00170 DOUBLE PRECISION ZERO, ONE 00171 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00172 * .. 00173 * .. Local Scalars .. 00174 INTEGER INFO, ISCL, NCOLS, NROWS 00175 DOUBLE PRECISION BIGNUM, ERR, NORMA, NORMB, NORMRS, NORMX, 00176 $ SMLNUM 00177 * .. 00178 * .. Local Arrays .. 00179 DOUBLE PRECISION RWORK( 1 ) 00180 * .. 00181 * .. External Functions .. 00182 LOGICAL LSAME 00183 DOUBLE PRECISION DLAMCH, DLANGE 00184 EXTERNAL LSAME, DLAMCH, DLANGE 00185 * .. 00186 * .. External Subroutines .. 00187 EXTERNAL DGEMM, DLACPY, DLASCL, XERBLA 00188 * .. 00189 * .. Intrinsic Functions .. 00190 INTRINSIC DBLE, MAX 00191 * .. 00192 * .. Executable Statements .. 00193 * 00194 DQRT17 = ZERO 00195 * 00196 IF( LSAME( TRANS, 'N' ) ) THEN 00197 NROWS = M 00198 NCOLS = N 00199 ELSE IF( LSAME( TRANS, 'T' ) ) THEN 00200 NROWS = N 00201 NCOLS = M 00202 ELSE 00203 CALL XERBLA( 'DQRT17', 1 ) 00204 RETURN 00205 END IF 00206 * 00207 IF( LWORK.LT.NCOLS*NRHS ) THEN 00208 CALL XERBLA( 'DQRT17', 13 ) 00209 RETURN 00210 END IF 00211 * 00212 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.LE.0 ) THEN 00213 RETURN 00214 END IF 00215 * 00216 NORMA = DLANGE( 'One-norm', M, N, A, LDA, RWORK ) 00217 SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' ) 00218 BIGNUM = ONE / SMLNUM 00219 ISCL = 0 00220 * 00221 * compute residual and scale it 00222 * 00223 CALL DLACPY( 'All', NROWS, NRHS, B, LDB, C, LDB ) 00224 CALL DGEMM( TRANS, 'No transpose', NROWS, NRHS, NCOLS, -ONE, A, 00225 $ LDA, X, LDX, ONE, C, LDB ) 00226 NORMRS = DLANGE( 'Max', NROWS, NRHS, C, LDB, RWORK ) 00227 IF( NORMRS.GT.SMLNUM ) THEN 00228 ISCL = 1 00229 CALL DLASCL( 'General', 0, 0, NORMRS, ONE, NROWS, NRHS, C, LDB, 00230 $ INFO ) 00231 END IF 00232 * 00233 * compute R'*A 00234 * 00235 CALL DGEMM( 'Transpose', TRANS, NRHS, NCOLS, NROWS, ONE, C, LDB, 00236 $ A, LDA, ZERO, WORK, NRHS ) 00237 * 00238 * compute and properly scale error 00239 * 00240 ERR = DLANGE( 'One-norm', NRHS, NCOLS, WORK, NRHS, RWORK ) 00241 IF( NORMA.NE.ZERO ) 00242 $ ERR = ERR / NORMA 00243 * 00244 IF( ISCL.EQ.1 ) 00245 $ ERR = ERR*NORMRS 00246 * 00247 IF( IRESID.EQ.1 ) THEN 00248 NORMB = DLANGE( 'One-norm', NROWS, NRHS, B, LDB, RWORK ) 00249 IF( NORMB.NE.ZERO ) 00250 $ ERR = ERR / NORMB 00251 ELSE 00252 NORMX = DLANGE( 'One-norm', NCOLS, NRHS, X, LDX, RWORK ) 00253 IF( NORMX.NE.ZERO ) 00254 $ ERR = ERR / NORMX 00255 END IF 00256 * 00257 DQRT17 = ERR / ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N, NRHS ) ) ) 00258 RETURN 00259 * 00260 * End of DQRT17 00261 * 00262 END