![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGLMTS 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 ZGLMTS( 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 *> ZGLMTS tests ZGGGLM - 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 COMPLEX*16 array, dimension (LDA,M) 00054 *> The N-by-M matrix A. 00055 *> \endverbatim 00056 *> 00057 *> \param[out] AF 00058 *> \verbatim 00059 *> AF is COMPLEX*16 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 COMPLEX*16 array, dimension (LDB,P) 00071 *> The N-by-P matrix A. 00072 *> \endverbatim 00073 *> 00074 *> \param[out] BF 00075 *> \verbatim 00076 *> BF is COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension( N ) 00094 *> \endverbatim 00095 *> 00096 *> \param[out] X 00097 *> \verbatim 00098 *> X is COMPLEX*16 array, dimension( M ) 00099 *> solution vector X in the GLM problem. 00100 *> \endverbatim 00101 *> 00102 *> \param[out] U 00103 *> \verbatim 00104 *> U is COMPLEX*16 array, dimension( P ) 00105 *> solution vector U in the GLM problem. 00106 *> \endverbatim 00107 *> 00108 *> \param[out] WORK 00109 *> \verbatim 00110 *> WORK is COMPLEX*16 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 complex16_eig 00144 * 00145 * ===================================================================== 00146 SUBROUTINE ZGLMTS( 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 RWORK( * ) 00163 COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ), 00164 $ BF( LDB, * ), D( * ), DF( * ), U( * ), 00165 $ WORK( LWORK ), X( * ) 00166 * .. 00167 * .. Parameters .. 00168 DOUBLE PRECISION ZERO 00169 PARAMETER ( ZERO = 0.0D+0 ) 00170 COMPLEX*16 CONE 00171 PARAMETER ( CONE = 1.0D+0 ) 00172 * .. 00173 * .. Local Scalars .. 00174 INTEGER INFO 00175 DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM 00176 * .. 00177 * .. External Functions .. 00178 DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE 00179 EXTERNAL DLAMCH, DZASUM, ZLANGE 00180 * .. 00181 * .. External Subroutines .. 00182 * 00183 EXTERNAL ZCOPY, ZGEMV, ZGGGLM, ZLACPY 00184 * .. 00185 * .. Intrinsic Functions .. 00186 INTRINSIC MAX 00187 * .. 00188 * .. Executable Statements .. 00189 * 00190 EPS = DLAMCH( 'Epsilon' ) 00191 UNFL = DLAMCH( 'Safe minimum' ) 00192 ANORM = MAX( ZLANGE( '1', N, M, A, LDA, RWORK ), UNFL ) 00193 BNORM = MAX( ZLANGE( '1', N, P, B, LDB, RWORK ), UNFL ) 00194 * 00195 * Copy the matrices A and B to the arrays AF and BF, 00196 * and the vector D the array DF. 00197 * 00198 CALL ZLACPY( 'Full', N, M, A, LDA, AF, LDA ) 00199 CALL ZLACPY( 'Full', N, P, B, LDB, BF, LDB ) 00200 CALL ZCOPY( N, D, 1, DF, 1 ) 00201 * 00202 * Solve GLM problem 00203 * 00204 CALL ZGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK, 00205 $ INFO ) 00206 * 00207 * Test the residual for the solution of LSE 00208 * 00209 * norm( d - A*x - B*u ) 00210 * RESULT = ----------------------------------------- 00211 * (norm(A)+norm(B))*(norm(x)+norm(u))*EPS 00212 * 00213 CALL ZCOPY( N, D, 1, DF, 1 ) 00214 CALL ZGEMV( 'No transpose', N, M, -CONE, A, LDA, X, 1, CONE, DF, 00215 $ 1 ) 00216 * 00217 CALL ZGEMV( 'No transpose', N, P, -CONE, B, LDB, U, 1, CONE, DF, 00218 $ 1 ) 00219 * 00220 DNORM = DZASUM( N, DF, 1 ) 00221 XNORM = DZASUM( M, X, 1 ) + DZASUM( P, U, 1 ) 00222 YNORM = ANORM + BNORM 00223 * 00224 IF( XNORM.LE.ZERO ) THEN 00225 RESULT = ZERO 00226 ELSE 00227 RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS 00228 END IF 00229 * 00230 RETURN 00231 * 00232 * End of ZGLMTS 00233 * 00234 END