![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGLMTS 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE DGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, 00012 * WORK, LWORK, RWORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER LDA, LDB, LWORK, M, N, P 00016 * DOUBLE PRECISION RESULT 00017 * .. 00018 * .. Array Arguments .. 00019 * 00020 * 00021 *> \par Purpose: 00022 * ============= 00023 *> 00024 *> \verbatim 00025 *> 00026 *> DGLMTS tests DGGGLM - a subroutine for solving the generalized 00027 *> linear model problem. 00028 *> \endverbatim 00029 * 00030 * Arguments: 00031 * ========== 00032 * 00033 *> \param[in] N 00034 *> \verbatim 00035 *> N is INTEGER 00036 *> The number of rows of the matrices A and B. N >= 0. 00037 *> \endverbatim 00038 *> 00039 *> \param[in] M 00040 *> \verbatim 00041 *> M is INTEGER 00042 *> The number of columns of the matrix A. M >= 0. 00043 *> \endverbatim 00044 *> 00045 *> \param[in] P 00046 *> \verbatim 00047 *> P is INTEGER 00048 *> The number of columns of the matrix B. P >= 0. 00049 *> \endverbatim 00050 *> 00051 *> \param[in] A 00052 *> \verbatim 00053 *> A is DOUBLE PRECISION array, dimension (LDA,M) 00054 *> The N-by-M matrix A. 00055 *> \endverbatim 00056 *> 00057 *> \param[out] AF 00058 *> \verbatim 00059 *> AF is DOUBLE PRECISION array, dimension (LDA,M) 00060 *> \endverbatim 00061 *> 00062 *> \param[in] LDA 00063 *> \verbatim 00064 *> LDA is INTEGER 00065 *> The leading dimension of the arrays A, AF. LDA >= max(M,N). 00066 *> \endverbatim 00067 *> 00068 *> \param[in] B 00069 *> \verbatim 00070 *> B is DOUBLE PRECISION array, dimension (LDB,P) 00071 *> The N-by-P matrix A. 00072 *> \endverbatim 00073 *> 00074 *> \param[out] BF 00075 *> \verbatim 00076 *> BF is DOUBLE PRECISION array, dimension (LDB,P) 00077 *> \endverbatim 00078 *> 00079 *> \param[in] LDB 00080 *> \verbatim 00081 *> LDB is INTEGER 00082 *> The leading dimension of the arrays B, BF. LDB >= max(P,N). 00083 *> \endverbatim 00084 *> 00085 *> \param[in] D 00086 *> \verbatim 00087 *> D is DOUBLE PRECISION array, dimension( N ) 00088 *> On input, the left hand side of the GLM. 00089 *> \endverbatim 00090 *> 00091 *> \param[out] DF 00092 *> \verbatim 00093 *> DF is DOUBLE PRECISION array, dimension( N ) 00094 *> \endverbatim 00095 *> 00096 *> \param[out] X 00097 *> \verbatim 00098 *> X is DOUBLE PRECISION array, dimension( M ) 00099 *> solution vector X in the GLM problem. 00100 *> \endverbatim 00101 *> 00102 *> \param[out] U 00103 *> \verbatim 00104 *> U is DOUBLE PRECISION array, dimension( P ) 00105 *> solution vector U in the GLM problem. 00106 *> \endverbatim 00107 *> 00108 *> \param[out] WORK 00109 *> \verbatim 00110 *> WORK is DOUBLE PRECISION array, dimension (LWORK) 00111 *> \endverbatim 00112 *> 00113 *> \param[in] LWORK 00114 *> \verbatim 00115 *> LWORK is INTEGER 00116 *> The dimension of the array WORK. 00117 *> \endverbatim 00118 *> 00119 *> \param[out] RWORK 00120 *> \verbatim 00121 *> RWORK is DOUBLE PRECISION array, dimension (M) 00122 *> \endverbatim 00123 *> 00124 *> \param[out] RESULT 00125 *> \verbatim 00126 *> RESULT is DOUBLE PRECISION 00127 *> The test ratio: 00128 *> norm( d - A*x - B*u ) 00129 *> RESULT = ----------------------------------------- 00130 *> (norm(A)+norm(B))*(norm(x)+norm(u))*EPS 00131 *> \endverbatim 00132 * 00133 * Authors: 00134 * ======== 00135 * 00136 *> \author Univ. of Tennessee 00137 *> \author Univ. of California Berkeley 00138 *> \author Univ. of Colorado Denver 00139 *> \author NAG Ltd. 00140 * 00141 *> \date November 2011 00142 * 00143 *> \ingroup double_eig 00144 * 00145 * ===================================================================== 00146 SUBROUTINE DGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, 00147 $ WORK, LWORK, RWORK, RESULT ) 00148 * 00149 * -- LAPACK test routine (version 3.4.0) -- 00150 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00151 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00152 * November 2011 00153 * 00154 * .. Scalar Arguments .. 00155 INTEGER LDA, LDB, LWORK, M, N, P 00156 DOUBLE PRECISION RESULT 00157 * .. 00158 * .. Array Arguments .. 00159 * 00160 * ==================================================================== 00161 * 00162 DOUBLE PRECISION A( LDA, * ), AF( LDA, * ), B( LDB, * ), 00163 $ BF( LDB, * ), D( * ), DF( * ), RWORK( * ), 00164 $ U( * ), WORK( LWORK ), X( * ) 00165 * .. 00166 * .. Parameters .. 00167 DOUBLE PRECISION ZERO, ONE 00168 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00169 * .. 00170 * .. Local Scalars .. 00171 INTEGER INFO 00172 DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM 00173 * .. 00174 * .. External Functions .. 00175 DOUBLE PRECISION DASUM, DLAMCH, DLANGE 00176 EXTERNAL DASUM, DLAMCH, DLANGE 00177 * .. 00178 * .. External Subroutines .. 00179 * 00180 EXTERNAL DCOPY, DGEMV, DGGGLM, DLACPY 00181 * .. 00182 * .. Intrinsic Functions .. 00183 INTRINSIC MAX 00184 * .. 00185 * .. Executable Statements .. 00186 * 00187 EPS = DLAMCH( 'Epsilon' ) 00188 UNFL = DLAMCH( 'Safe minimum' ) 00189 ANORM = MAX( DLANGE( '1', N, M, A, LDA, RWORK ), UNFL ) 00190 BNORM = MAX( DLANGE( '1', N, P, B, LDB, RWORK ), UNFL ) 00191 * 00192 * Copy the matrices A and B to the arrays AF and BF, 00193 * and the vector D the array DF. 00194 * 00195 CALL DLACPY( 'Full', N, M, A, LDA, AF, LDA ) 00196 CALL DLACPY( 'Full', N, P, B, LDB, BF, LDB ) 00197 CALL DCOPY( N, D, 1, DF, 1 ) 00198 * 00199 * Solve GLM problem 00200 * 00201 CALL DGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK, 00202 $ INFO ) 00203 * 00204 * Test the residual for the solution of LSE 00205 * 00206 * norm( d - A*x - B*u ) 00207 * RESULT = ----------------------------------------- 00208 * (norm(A)+norm(B))*(norm(x)+norm(u))*EPS 00209 * 00210 CALL DCOPY( N, D, 1, DF, 1 ) 00211 CALL DGEMV( 'No transpose', N, M, -ONE, A, LDA, X, 1, ONE, DF, 1 ) 00212 * 00213 CALL DGEMV( 'No transpose', N, P, -ONE, B, LDB, U, 1, ONE, DF, 1 ) 00214 * 00215 DNORM = DASUM( N, DF, 1 ) 00216 XNORM = DASUM( M, X, 1 ) + DASUM( P, U, 1 ) 00217 YNORM = ANORM + BNORM 00218 * 00219 IF( XNORM.LE.ZERO ) THEN 00220 RESULT = ZERO 00221 ELSE 00222 RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS 00223 END IF 00224 * 00225 RETURN 00226 * 00227 * End of DGLMTS 00228 * 00229 END