1*> \brief <b> DGELS solves overdetermined or underdetermined systems for GE matrices</b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DGELS + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgels.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgels.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgels.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, 22* INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER TRANS 26* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS 27* .. 28* .. Array Arguments .. 29* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> DGELS solves overdetermined or underdetermined real linear systems 39*> involving an M-by-N matrix A, or its transpose, using a QR or LQ 40*> factorization of A. It is assumed that A has full rank. 41*> 42*> The following options are provided: 43*> 44*> 1. If TRANS = 'N' and m >= n: find the least squares solution of 45*> an overdetermined system, i.e., solve the least squares problem 46*> minimize || B - A*X ||. 47*> 48*> 2. If TRANS = 'N' and m < n: find the minimum norm solution of 49*> an underdetermined system A * X = B. 50*> 51*> 3. If TRANS = 'T' and m >= n: find the minimum norm solution of 52*> an underdetermined system A**T * X = B. 53*> 54*> 4. If TRANS = 'T' and m < n: find the least squares solution of 55*> an overdetermined system, i.e., solve the least squares problem 56*> minimize || B - A**T * X ||. 57*> 58*> Several right hand side vectors b and solution vectors x can be 59*> handled in a single call; they are stored as the columns of the 60*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution 61*> matrix X. 62*> \endverbatim 63* 64* Arguments: 65* ========== 66* 67*> \param[in] TRANS 68*> \verbatim 69*> TRANS is CHARACTER*1 70*> = 'N': the linear system involves A; 71*> = 'T': the linear system involves A**T. 72*> \endverbatim 73*> 74*> \param[in] M 75*> \verbatim 76*> M is INTEGER 77*> The number of rows of the matrix A. M >= 0. 78*> \endverbatim 79*> 80*> \param[in] N 81*> \verbatim 82*> N is INTEGER 83*> The number of columns of the matrix A. N >= 0. 84*> \endverbatim 85*> 86*> \param[in] NRHS 87*> \verbatim 88*> NRHS is INTEGER 89*> The number of right hand sides, i.e., the number of 90*> columns of the matrices B and X. NRHS >=0. 91*> \endverbatim 92*> 93*> \param[in,out] A 94*> \verbatim 95*> A is DOUBLE PRECISION array, dimension (LDA,N) 96*> On entry, the M-by-N matrix A. 97*> On exit, 98*> if M >= N, A is overwritten by details of its QR 99*> factorization as returned by DGEQRF; 100*> if M < N, A is overwritten by details of its LQ 101*> factorization as returned by DGELQF. 102*> \endverbatim 103*> 104*> \param[in] LDA 105*> \verbatim 106*> LDA is INTEGER 107*> The leading dimension of the array A. LDA >= max(1,M). 108*> \endverbatim 109*> 110*> \param[in,out] B 111*> \verbatim 112*> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 113*> On entry, the matrix B of right hand side vectors, stored 114*> columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS 115*> if TRANS = 'T'. 116*> On exit, if INFO = 0, B is overwritten by the solution 117*> vectors, stored columnwise: 118*> if TRANS = 'N' and m >= n, rows 1 to n of B contain the least 119*> squares solution vectors; the residual sum of squares for the 120*> solution in each column is given by the sum of squares of 121*> elements N+1 to M in that column; 122*> if TRANS = 'N' and m < n, rows 1 to N of B contain the 123*> minimum norm solution vectors; 124*> if TRANS = 'T' and m >= n, rows 1 to M of B contain the 125*> minimum norm solution vectors; 126*> if TRANS = 'T' and m < n, rows 1 to M of B contain the 127*> least squares solution vectors; the residual sum of squares 128*> for the solution in each column is given by the sum of 129*> squares of elements M+1 to N in that column. 130*> \endverbatim 131*> 132*> \param[in] LDB 133*> \verbatim 134*> LDB is INTEGER 135*> The leading dimension of the array B. LDB >= MAX(1,M,N). 136*> \endverbatim 137*> 138*> \param[out] WORK 139*> \verbatim 140*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 141*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 142*> \endverbatim 143*> 144*> \param[in] LWORK 145*> \verbatim 146*> LWORK is INTEGER 147*> The dimension of the array WORK. 148*> LWORK >= max( 1, MN + max( MN, NRHS ) ). 149*> For optimal performance, 150*> LWORK >= max( 1, MN + max( MN, NRHS )*NB ). 151*> where MN = min(M,N) and NB is the optimum block size. 152*> 153*> If LWORK = -1, then a workspace query is assumed; the routine 154*> only calculates the optimal size of the WORK array, returns 155*> this value as the first entry of the WORK array, and no error 156*> message related to LWORK is issued by XERBLA. 157*> \endverbatim 158*> 159*> \param[out] INFO 160*> \verbatim 161*> INFO is INTEGER 162*> = 0: successful exit 163*> < 0: if INFO = -i, the i-th argument had an illegal value 164*> > 0: if INFO = i, the i-th diagonal element of the 165*> triangular factor of A is zero, so that A does not have 166*> full rank; the least squares solution could not be 167*> computed. 168*> \endverbatim 169* 170* Authors: 171* ======== 172* 173*> \author Univ. of Tennessee 174*> \author Univ. of California Berkeley 175*> \author Univ. of Colorado Denver 176*> \author NAG Ltd. 177* 178*> \ingroup doubleGEsolve 179* 180* ===================================================================== 181 SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, 182 $ INFO ) 183* 184* -- LAPACK driver routine -- 185* -- LAPACK is a software package provided by Univ. of Tennessee, -- 186* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 187* 188* .. Scalar Arguments .. 189 CHARACTER TRANS 190 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS 191* .. 192* .. Array Arguments .. 193 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) 194* .. 195* 196* ===================================================================== 197* 198* .. Parameters .. 199 DOUBLE PRECISION ZERO, ONE 200 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 201* .. 202* .. Local Scalars .. 203 LOGICAL LQUERY, TPSD 204 INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE 205 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM 206* .. 207* .. Local Arrays .. 208 DOUBLE PRECISION RWORK( 1 ) 209* .. 210* .. External Functions .. 211 LOGICAL LSAME 212 INTEGER ILAENV 213 DOUBLE PRECISION DLAMCH, DLANGE 214 EXTERNAL LSAME, ILAENV, DLABAD, DLAMCH, DLANGE 215* .. 216* .. External Subroutines .. 217 EXTERNAL DGELQF, DGEQRF, DLASCL, DLASET, DORMLQ, DORMQR, 218 $ DTRTRS, XERBLA 219* .. 220* .. Intrinsic Functions .. 221 INTRINSIC DBLE, MAX, MIN 222* .. 223* .. Executable Statements .. 224* 225* Test the input arguments. 226* 227 INFO = 0 228 MN = MIN( M, N ) 229 LQUERY = ( LWORK.EQ.-1 ) 230 IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'T' ) ) ) THEN 231 INFO = -1 232 ELSE IF( M.LT.0 ) THEN 233 INFO = -2 234 ELSE IF( N.LT.0 ) THEN 235 INFO = -3 236 ELSE IF( NRHS.LT.0 ) THEN 237 INFO = -4 238 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 239 INFO = -6 240 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN 241 INFO = -8 242 ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY ) 243 $ THEN 244 INFO = -10 245 END IF 246* 247* Figure out optimal block size 248* 249 IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN 250* 251 TPSD = .TRUE. 252 IF( LSAME( TRANS, 'N' ) ) 253 $ TPSD = .FALSE. 254* 255 IF( M.GE.N ) THEN 256 NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) 257 IF( TPSD ) THEN 258 NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LN', M, NRHS, N, 259 $ -1 ) ) 260 ELSE 261 NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N, 262 $ -1 ) ) 263 END IF 264 ELSE 265 NB = ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) 266 IF( TPSD ) THEN 267 NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M, 268 $ -1 ) ) 269 ELSE 270 NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LN', N, NRHS, M, 271 $ -1 ) ) 272 END IF 273 END IF 274* 275 WSIZE = MAX( 1, MN+MAX( MN, NRHS )*NB ) 276 WORK( 1 ) = DBLE( WSIZE ) 277* 278 END IF 279* 280 IF( INFO.NE.0 ) THEN 281 CALL XERBLA( 'DGELS ', -INFO ) 282 RETURN 283 ELSE IF( LQUERY ) THEN 284 RETURN 285 END IF 286* 287* Quick return if possible 288* 289 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 290 CALL DLASET( 'Full', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 291 RETURN 292 END IF 293* 294* Get machine parameters 295* 296 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) 297 BIGNUM = ONE / SMLNUM 298 CALL DLABAD( SMLNUM, BIGNUM ) 299* 300* Scale A, B if max element outside range [SMLNUM,BIGNUM] 301* 302 ANRM = DLANGE( 'M', M, N, A, LDA, RWORK ) 303 IASCL = 0 304 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 305* 306* Scale matrix norm up to SMLNUM 307* 308 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 309 IASCL = 1 310 ELSE IF( ANRM.GT.BIGNUM ) THEN 311* 312* Scale matrix norm down to BIGNUM 313* 314 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 315 IASCL = 2 316 ELSE IF( ANRM.EQ.ZERO ) THEN 317* 318* Matrix all zero. Return zero solution. 319* 320 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 321 GO TO 50 322 END IF 323* 324 BROW = M 325 IF( TPSD ) 326 $ BROW = N 327 BNRM = DLANGE( 'M', BROW, NRHS, B, LDB, RWORK ) 328 IBSCL = 0 329 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 330* 331* Scale matrix norm up to SMLNUM 332* 333 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB, 334 $ INFO ) 335 IBSCL = 1 336 ELSE IF( BNRM.GT.BIGNUM ) THEN 337* 338* Scale matrix norm down to BIGNUM 339* 340 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB, 341 $ INFO ) 342 IBSCL = 2 343 END IF 344* 345 IF( M.GE.N ) THEN 346* 347* compute QR factorization of A 348* 349 CALL DGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN, 350 $ INFO ) 351* 352* workspace at least N, optimally N*NB 353* 354 IF( .NOT.TPSD ) THEN 355* 356* Least-Squares Problem min || A * X - B || 357* 358* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS) 359* 360 CALL DORMQR( 'Left', 'Transpose', M, NRHS, N, A, LDA, 361 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, 362 $ INFO ) 363* 364* workspace at least NRHS, optimally NRHS*NB 365* 366* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS) 367* 368 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS, 369 $ A, LDA, B, LDB, INFO ) 370* 371 IF( INFO.GT.0 ) THEN 372 RETURN 373 END IF 374* 375 SCLLEN = N 376* 377 ELSE 378* 379* Underdetermined system of equations A**T * X = B 380* 381* B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS) 382* 383 CALL DTRTRS( 'Upper', 'Transpose', 'Non-unit', N, NRHS, 384 $ A, LDA, B, LDB, INFO ) 385* 386 IF( INFO.GT.0 ) THEN 387 RETURN 388 END IF 389* 390* B(N+1:M,1:NRHS) = ZERO 391* 392 DO 20 J = 1, NRHS 393 DO 10 I = N + 1, M 394 B( I, J ) = ZERO 395 10 CONTINUE 396 20 CONTINUE 397* 398* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS) 399* 400 CALL DORMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA, 401 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, 402 $ INFO ) 403* 404* workspace at least NRHS, optimally NRHS*NB 405* 406 SCLLEN = M 407* 408 END IF 409* 410 ELSE 411* 412* Compute LQ factorization of A 413* 414 CALL DGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN, 415 $ INFO ) 416* 417* workspace at least M, optimally M*NB. 418* 419 IF( .NOT.TPSD ) THEN 420* 421* underdetermined system of equations A * X = B 422* 423* B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS) 424* 425 CALL DTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS, 426 $ A, LDA, B, LDB, INFO ) 427* 428 IF( INFO.GT.0 ) THEN 429 RETURN 430 END IF 431* 432* B(M+1:N,1:NRHS) = 0 433* 434 DO 40 J = 1, NRHS 435 DO 30 I = M + 1, N 436 B( I, J ) = ZERO 437 30 CONTINUE 438 40 CONTINUE 439* 440* B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS) 441* 442 CALL DORMLQ( 'Left', 'Transpose', N, NRHS, M, A, LDA, 443 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, 444 $ INFO ) 445* 446* workspace at least NRHS, optimally NRHS*NB 447* 448 SCLLEN = N 449* 450 ELSE 451* 452* overdetermined system min || A**T * X - B || 453* 454* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS) 455* 456 CALL DORMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA, 457 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, 458 $ INFO ) 459* 460* workspace at least NRHS, optimally NRHS*NB 461* 462* B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS) 463* 464 CALL DTRTRS( 'Lower', 'Transpose', 'Non-unit', M, NRHS, 465 $ A, LDA, B, LDB, INFO ) 466* 467 IF( INFO.GT.0 ) THEN 468 RETURN 469 END IF 470* 471 SCLLEN = M 472* 473 END IF 474* 475 END IF 476* 477* Undo scaling 478* 479 IF( IASCL.EQ.1 ) THEN 480 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB, 481 $ INFO ) 482 ELSE IF( IASCL.EQ.2 ) THEN 483 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB, 484 $ INFO ) 485 END IF 486 IF( IBSCL.EQ.1 ) THEN 487 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB, 488 $ INFO ) 489 ELSE IF( IBSCL.EQ.2 ) THEN 490 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB, 491 $ INFO ) 492 END IF 493* 494 50 CONTINUE 495 WORK( 1 ) = DBLE( WSIZE ) 496* 497 RETURN 498* 499* End of DGELS 500* 501 END 502