![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> SGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SGGGLM + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggglm.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggglm.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggglm.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, 00022 * INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDA, LDB, LWORK, M, N, P 00026 * .. 00027 * .. Array Arguments .. 00028 * REAL A( LDA, * ), B( LDB, * ), D( * ), WORK( * ), 00029 * $ X( * ), Y( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> SGGGLM solves a general Gauss-Markov linear model (GLM) problem: 00039 *> 00040 *> minimize || y ||_2 subject to d = A*x + B*y 00041 *> x 00042 *> 00043 *> where A is an N-by-M matrix, B is an N-by-P matrix, and d is a 00044 *> given N-vector. It is assumed that M <= N <= M+P, and 00045 *> 00046 *> rank(A) = M and rank( A B ) = N. 00047 *> 00048 *> Under these assumptions, the constrained equation is always 00049 *> consistent, and there is a unique solution x and a minimal 2-norm 00050 *> solution y, which is obtained using a generalized QR factorization 00051 *> of the matrices (A, B) given by 00052 *> 00053 *> A = Q*(R), B = Q*T*Z. 00054 *> (0) 00055 *> 00056 *> In particular, if matrix B is square nonsingular, then the problem 00057 *> GLM is equivalent to the following weighted linear least squares 00058 *> problem 00059 *> 00060 *> minimize || inv(B)*(d-A*x) ||_2 00061 *> x 00062 *> 00063 *> where inv(B) denotes the inverse of B. 00064 *> \endverbatim 00065 * 00066 * Arguments: 00067 * ========== 00068 * 00069 *> \param[in] N 00070 *> \verbatim 00071 *> N is INTEGER 00072 *> The number of rows of the matrices A and B. N >= 0. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] M 00076 *> \verbatim 00077 *> M is INTEGER 00078 *> The number of columns of the matrix A. 0 <= M <= N. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] P 00082 *> \verbatim 00083 *> P is INTEGER 00084 *> The number of columns of the matrix B. P >= N-M. 00085 *> \endverbatim 00086 *> 00087 *> \param[in,out] A 00088 *> \verbatim 00089 *> A is REAL array, dimension (LDA,M) 00090 *> On entry, the N-by-M matrix A. 00091 *> On exit, the upper triangular part of the array A contains 00092 *> the M-by-M upper triangular matrix R. 00093 *> \endverbatim 00094 *> 00095 *> \param[in] LDA 00096 *> \verbatim 00097 *> LDA is INTEGER 00098 *> The leading dimension of the array A. LDA >= max(1,N). 00099 *> \endverbatim 00100 *> 00101 *> \param[in,out] B 00102 *> \verbatim 00103 *> B is REAL array, dimension (LDB,P) 00104 *> On entry, the N-by-P matrix B. 00105 *> On exit, if N <= P, the upper triangle of the subarray 00106 *> B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T; 00107 *> if N > P, the elements on and above the (N-P)th subdiagonal 00108 *> contain the N-by-P upper trapezoidal matrix T. 00109 *> \endverbatim 00110 *> 00111 *> \param[in] LDB 00112 *> \verbatim 00113 *> LDB is INTEGER 00114 *> The leading dimension of the array B. LDB >= max(1,N). 00115 *> \endverbatim 00116 *> 00117 *> \param[in,out] D 00118 *> \verbatim 00119 *> D is REAL array, dimension (N) 00120 *> On entry, D is the left hand side of the GLM equation. 00121 *> On exit, D is destroyed. 00122 *> \endverbatim 00123 *> 00124 *> \param[out] X 00125 *> \verbatim 00126 *> X is REAL array, dimension (M) 00127 *> \endverbatim 00128 *> 00129 *> \param[out] Y 00130 *> \verbatim 00131 *> Y is REAL array, dimension (P) 00132 *> 00133 *> On exit, X and Y are the solutions of the GLM problem. 00134 *> \endverbatim 00135 *> 00136 *> \param[out] WORK 00137 *> \verbatim 00138 *> WORK is REAL array, dimension (MAX(1,LWORK)) 00139 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00140 *> \endverbatim 00141 *> 00142 *> \param[in] LWORK 00143 *> \verbatim 00144 *> LWORK is INTEGER 00145 *> The dimension of the array WORK. LWORK >= max(1,N+M+P). 00146 *> For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB, 00147 *> where NB is an upper bound for the optimal blocksizes for 00148 *> SGEQRF, SGERQF, SORMQR and SORMRQ. 00149 *> 00150 *> If LWORK = -1, then a workspace query is assumed; the routine 00151 *> only calculates the optimal size of the WORK array, returns 00152 *> this value as the first entry of the WORK array, and no error 00153 *> message related to LWORK is issued by XERBLA. 00154 *> \endverbatim 00155 *> 00156 *> \param[out] INFO 00157 *> \verbatim 00158 *> INFO is INTEGER 00159 *> = 0: successful exit. 00160 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00161 *> = 1: the upper triangular factor R associated with A in the 00162 *> generalized QR factorization of the pair (A, B) is 00163 *> singular, so that rank(A) < M; the least squares 00164 *> solution could not be computed. 00165 *> = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal 00166 *> factor T associated with B in the generalized QR 00167 *> factorization of the pair (A, B) is singular, so that 00168 *> rank( A B ) < N; the least squares solution could not 00169 *> be computed. 00170 *> \endverbatim 00171 * 00172 * Authors: 00173 * ======== 00174 * 00175 *> \author Univ. of Tennessee 00176 *> \author Univ. of California Berkeley 00177 *> \author Univ. of Colorado Denver 00178 *> \author NAG Ltd. 00179 * 00180 *> \date November 2011 00181 * 00182 *> \ingroup realOTHEReigen 00183 * 00184 * ===================================================================== 00185 SUBROUTINE SGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, 00186 $ INFO ) 00187 * 00188 * -- LAPACK driver routine (version 3.4.0) -- 00189 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00190 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00191 * November 2011 00192 * 00193 * .. Scalar Arguments .. 00194 INTEGER INFO, LDA, LDB, LWORK, M, N, P 00195 * .. 00196 * .. Array Arguments .. 00197 REAL A( LDA, * ), B( LDB, * ), D( * ), WORK( * ), 00198 $ X( * ), Y( * ) 00199 * .. 00200 * 00201 * =================================================================== 00202 * 00203 * .. Parameters .. 00204 REAL ZERO, ONE 00205 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00206 * .. 00207 * .. Local Scalars .. 00208 LOGICAL LQUERY 00209 INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3, 00210 $ NB4, NP 00211 * .. 00212 * .. External Subroutines .. 00213 EXTERNAL SCOPY, SGEMV, SGGQRF, SORMQR, SORMRQ, STRTRS, 00214 $ XERBLA 00215 * .. 00216 * .. External Functions .. 00217 INTEGER ILAENV 00218 EXTERNAL ILAENV 00219 * .. 00220 * .. Intrinsic Functions .. 00221 INTRINSIC INT, MAX, MIN 00222 * .. 00223 * .. Executable Statements .. 00224 * 00225 * Test the input parameters 00226 * 00227 INFO = 0 00228 NP = MIN( N, P ) 00229 LQUERY = ( LWORK.EQ.-1 ) 00230 IF( N.LT.0 ) THEN 00231 INFO = -1 00232 ELSE IF( M.LT.0 .OR. M.GT.N ) THEN 00233 INFO = -2 00234 ELSE IF( P.LT.0 .OR. P.LT.N-M ) THEN 00235 INFO = -3 00236 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00237 INFO = -5 00238 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00239 INFO = -7 00240 END IF 00241 * 00242 * Calculate workspace 00243 * 00244 IF( INFO.EQ.0) THEN 00245 IF( N.EQ.0 ) THEN 00246 LWKMIN = 1 00247 LWKOPT = 1 00248 ELSE 00249 NB1 = ILAENV( 1, 'SGEQRF', ' ', N, M, -1, -1 ) 00250 NB2 = ILAENV( 1, 'SGERQF', ' ', N, M, -1, -1 ) 00251 NB3 = ILAENV( 1, 'SORMQR', ' ', N, M, P, -1 ) 00252 NB4 = ILAENV( 1, 'SORMRQ', ' ', N, M, P, -1 ) 00253 NB = MAX( NB1, NB2, NB3, NB4 ) 00254 LWKMIN = M + N + P 00255 LWKOPT = M + NP + MAX( N, P )*NB 00256 END IF 00257 WORK( 1 ) = LWKOPT 00258 * 00259 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN 00260 INFO = -12 00261 END IF 00262 END IF 00263 * 00264 IF( INFO.NE.0 ) THEN 00265 CALL XERBLA( 'SGGGLM', -INFO ) 00266 RETURN 00267 ELSE IF( LQUERY ) THEN 00268 RETURN 00269 END IF 00270 * 00271 * Quick return if possible 00272 * 00273 IF( N.EQ.0 ) 00274 $ RETURN 00275 * 00276 * Compute the GQR factorization of matrices A and B: 00277 * 00278 * Q**T*A = ( R11 ) M, Q**T*B*Z**T = ( T11 T12 ) M 00279 * ( 0 ) N-M ( 0 T22 ) N-M 00280 * M M+P-N N-M 00281 * 00282 * where R11 and T22 are upper triangular, and Q and Z are 00283 * orthogonal. 00284 * 00285 CALL SGGQRF( N, M, P, A, LDA, WORK, B, LDB, WORK( M+1 ), 00286 $ WORK( M+NP+1 ), LWORK-M-NP, INFO ) 00287 LOPT = WORK( M+NP+1 ) 00288 * 00289 * Update left-hand-side vector d = Q**T*d = ( d1 ) M 00290 * ( d2 ) N-M 00291 * 00292 CALL SORMQR( 'Left', 'Transpose', N, 1, M, A, LDA, WORK, D, 00293 $ MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) 00294 LOPT = MAX( LOPT, INT( WORK( M+NP+1 ) ) ) 00295 * 00296 * Solve T22*y2 = d2 for y2 00297 * 00298 IF( N.GT.M ) THEN 00299 CALL STRTRS( 'Upper', 'No transpose', 'Non unit', N-M, 1, 00300 $ B( M+1, M+P-N+1 ), LDB, D( M+1 ), N-M, INFO ) 00301 * 00302 IF( INFO.GT.0 ) THEN 00303 INFO = 1 00304 RETURN 00305 END IF 00306 * 00307 CALL SCOPY( N-M, D( M+1 ), 1, Y( M+P-N+1 ), 1 ) 00308 END IF 00309 * 00310 * Set y1 = 0 00311 * 00312 DO 10 I = 1, M + P - N 00313 Y( I ) = ZERO 00314 10 CONTINUE 00315 * 00316 * Update d1 = d1 - T12*y2 00317 * 00318 CALL SGEMV( 'No transpose', M, N-M, -ONE, B( 1, M+P-N+1 ), LDB, 00319 $ Y( M+P-N+1 ), 1, ONE, D, 1 ) 00320 * 00321 * Solve triangular system: R11*x = d1 00322 * 00323 IF( M.GT.0 ) THEN 00324 CALL STRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, LDA, 00325 $ D, M, INFO ) 00326 * 00327 IF( INFO.GT.0 ) THEN 00328 INFO = 2 00329 RETURN 00330 END IF 00331 * 00332 * Copy D to X 00333 * 00334 CALL SCOPY( M, D, 1, X, 1 ) 00335 END IF 00336 * 00337 * Backward transformation y = Z**T *y 00338 * 00339 CALL SORMRQ( 'Left', 'Transpose', P, 1, NP, 00340 $ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y, 00341 $ MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) 00342 WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) ) 00343 * 00344 RETURN 00345 * 00346 * End of SGGGLM 00347 * 00348 END