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