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