1 SUBROUTINE DGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, 2 $ WORK, 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* March 31, 1993 8* 9* .. Scalar Arguments .. 10 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK 11 DOUBLE PRECISION RCOND 12* .. 13* .. Array Arguments .. 14 INTEGER JPVT( * ) 15 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* This routine is deprecated and has been replaced by routine DGELSY. 22* 23* DGELSX computes the minimum-norm solution to a real linear least 24* squares problem: 25* minimize || A * X - B || 26* using a complete orthogonal factorization of A. A is an M-by-N 27* matrix which may be rank-deficient. 28* 29* Several right hand side vectors b and solution vectors x can be 30* handled in a single call; they are stored as the columns of the 31* M-by-NRHS right hand side matrix B and the N-by-NRHS solution 32* matrix X. 33* 34* The routine first computes a QR factorization with column pivoting: 35* A * P = Q * [ R11 R12 ] 36* [ 0 R22 ] 37* with R11 defined as the largest leading submatrix whose estimated 38* condition number is less than 1/RCOND. The order of R11, RANK, 39* is the effective rank of A. 40* 41* Then, R22 is considered to be negligible, and R12 is annihilated 42* by orthogonal transformations from the right, arriving at the 43* complete orthogonal factorization: 44* A * P = Q * [ T11 0 ] * Z 45* [ 0 0 ] 46* The minimum-norm solution is then 47* X = P * Z' [ inv(T11)*Q1'*B ] 48* [ 0 ] 49* where Q1 consists of the first RANK columns of Q. 50* 51* Arguments 52* ========= 53* 54* M (input) INTEGER 55* The number of rows of the matrix A. M >= 0. 56* 57* N (input) INTEGER 58* The number of columns of the matrix A. N >= 0. 59* 60* NRHS (input) INTEGER 61* The number of right hand sides, i.e., the number of 62* columns of matrices B and X. NRHS >= 0. 63* 64* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 65* On entry, the M-by-N matrix A. 66* On exit, A has been overwritten by details of its 67* complete orthogonal factorization. 68* 69* LDA (input) INTEGER 70* The leading dimension of the array A. LDA >= max(1,M). 71* 72* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 73* On entry, the M-by-NRHS right hand side matrix B. 74* On exit, the N-by-NRHS solution matrix X. 75* If m >= n and RANK = n, the residual sum-of-squares for 76* the solution in the i-th column is given by the sum of 77* squares of elements N+1:M in that column. 78* 79* LDB (input) INTEGER 80* The leading dimension of the array B. LDB >= max(1,M,N). 81* 82* JPVT (input/output) INTEGER array, dimension (N) 83* On entry, if JPVT(i) .ne. 0, the i-th column of A is an 84* initial column, otherwise it is a free column. Before 85* the QR factorization of A, all initial columns are 86* permuted to the leading positions; only the remaining 87* free columns are moved as a result of column pivoting 88* during the factorization. 89* On exit, if JPVT(i) = k, then the i-th column of A*P 90* was the k-th column of A. 91* 92* RCOND (input) DOUBLE PRECISION 93* RCOND is used to determine the effective rank of A, which 94* is defined as the order of the largest leading triangular 95* submatrix R11 in the QR factorization with pivoting of A, 96* whose estimated condition number < 1/RCOND. 97* 98* RANK (output) INTEGER 99* The effective rank of A, i.e., the order of the submatrix 100* R11. This is the same as the order of the submatrix T11 101* in the complete orthogonal factorization of A. 102* 103* WORK (workspace) DOUBLE PRECISION array, dimension 104* (max( min(M,N)+3*N, 2*min(M,N)+NRHS )), 105* 106* INFO (output) INTEGER 107* = 0: successful exit 108* < 0: if INFO = -i, the i-th argument had an illegal value 109* 110* ===================================================================== 111* 112* .. Parameters .. 113 INTEGER IMAX, IMIN 114 PARAMETER ( IMAX = 1, IMIN = 2 ) 115 DOUBLE PRECISION ZERO, ONE, DONE, NTDONE 116 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, DONE = ZERO, 117 $ NTDONE = ONE ) 118* .. 119* .. Local Scalars .. 120 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN 121 DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX, 122 $ SMAXPR, SMIN, SMINPR, SMLNUM, T1, T2 123* .. 124* .. External Functions .. 125 DOUBLE PRECISION DLAMCH, DLANGE 126 EXTERNAL DLAMCH, DLANGE 127* .. 128* .. External Subroutines .. 129 EXTERNAL DGEQPF, DLAIC1, DLASCL, DLASET, DLATZM, DORM2R, 130 $ DTRSM, DTZRQF, XERBLA 131* .. 132* .. Intrinsic Functions .. 133 INTRINSIC ABS, MAX, MIN 134* .. 135* .. Executable Statements .. 136* 137 MN = MIN( M, N ) 138 ISMIN = MN + 1 139 ISMAX = 2*MN + 1 140* 141* Test the input arguments. 142* 143 INFO = 0 144 IF( M.LT.0 ) THEN 145 INFO = -1 146 ELSE IF( N.LT.0 ) THEN 147 INFO = -2 148 ELSE IF( NRHS.LT.0 ) THEN 149 INFO = -3 150 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 151 INFO = -5 152 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN 153 INFO = -7 154 END IF 155* 156 IF( INFO.NE.0 ) THEN 157 CALL XERBLA( 'DGELSX', -INFO ) 158 RETURN 159 END IF 160* 161* Quick return if possible 162* 163 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 164 RANK = 0 165 RETURN 166 END IF 167* 168* Get machine parameters 169* 170 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) 171 BIGNUM = ONE / SMLNUM 172 CALL DLABAD( SMLNUM, BIGNUM ) 173* 174* Scale A, B if max elements outside range [SMLNUM,BIGNUM] 175* 176 ANRM = DLANGE( 'M', M, N, A, LDA, WORK ) 177 IASCL = 0 178 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 179* 180* Scale matrix norm up to SMLNUM 181* 182 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 183 IASCL = 1 184 ELSE IF( ANRM.GT.BIGNUM ) THEN 185* 186* Scale matrix norm down to BIGNUM 187* 188 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 189 IASCL = 2 190 ELSE IF( ANRM.EQ.ZERO ) THEN 191* 192* Matrix all zero. Return zero solution. 193* 194 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 195 RANK = 0 196 GO TO 100 197 END IF 198* 199 BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK ) 200 IBSCL = 0 201 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 202* 203* Scale matrix norm up to SMLNUM 204* 205 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 206 IBSCL = 1 207 ELSE IF( BNRM.GT.BIGNUM ) THEN 208* 209* Scale matrix norm down to BIGNUM 210* 211 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 212 IBSCL = 2 213 END IF 214* 215* Compute QR factorization with column pivoting of A: 216* A * P = Q * R 217* 218 CALL DGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), INFO ) 219* 220* workspace 3*N. Details of Householder rotations stored 221* in WORK(1:MN). 222* 223* Determine RANK using incremental condition estimation 224* 225 WORK( ISMIN ) = ONE 226 WORK( ISMAX ) = ONE 227 SMAX = ABS( A( 1, 1 ) ) 228 SMIN = SMAX 229 IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN 230 RANK = 0 231 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) 232 GO TO 100 233 ELSE 234 RANK = 1 235 END IF 236* 237 10 CONTINUE 238 IF( RANK.LT.MN ) THEN 239 I = RANK + 1 240 CALL DLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ), 241 $ A( I, I ), SMINPR, S1, C1 ) 242 CALL DLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ), 243 $ A( I, I ), SMAXPR, S2, C2 ) 244* 245 IF( SMAXPR*RCOND.LE.SMINPR ) THEN 246 DO 20 I = 1, RANK 247 WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 ) 248 WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 ) 249 20 CONTINUE 250 WORK( ISMIN+RANK ) = C1 251 WORK( ISMAX+RANK ) = C2 252 SMIN = SMINPR 253 SMAX = SMAXPR 254 RANK = RANK + 1 255 GO TO 10 256 END IF 257 END IF 258* 259* Logically partition R = [ R11 R12 ] 260* [ 0 R22 ] 261* where R11 = R(1:RANK,1:RANK) 262* 263* [R11,R12] = [ T11, 0 ] * Y 264* 265 IF( RANK.LT.N ) 266 $ CALL DTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO ) 267* 268* Details of Householder rotations stored in WORK(MN+1:2*MN) 269* 270* B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) 271* 272 CALL DORM2R( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ), 273 $ B, LDB, WORK( 2*MN+1 ), INFO ) 274* 275* workspace NRHS 276* 277* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) 278* 279 CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK, 280 $ NRHS, ONE, A, LDA, B, LDB ) 281* 282 DO 40 I = RANK + 1, N 283 DO 30 J = 1, NRHS 284 B( I, J ) = ZERO 285 30 CONTINUE 286 40 CONTINUE 287* 288* B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) 289* 290 IF( RANK.LT.N ) THEN 291 DO 50 I = 1, RANK 292 CALL DLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA, 293 $ WORK( MN+I ), B( I, 1 ), B( RANK+1, 1 ), LDB, 294 $ WORK( 2*MN+1 ) ) 295 50 CONTINUE 296 END IF 297* 298* workspace NRHS 299* 300* B(1:N,1:NRHS) := P * B(1:N,1:NRHS) 301* 302 DO 90 J = 1, NRHS 303 DO 60 I = 1, N 304 WORK( 2*MN+I ) = NTDONE 305 60 CONTINUE 306 DO 80 I = 1, N 307 IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN 308 IF( JPVT( I ).NE.I ) THEN 309 K = I 310 T1 = B( K, J ) 311 T2 = B( JPVT( K ), J ) 312 70 CONTINUE 313 B( JPVT( K ), J ) = T1 314 WORK( 2*MN+K ) = DONE 315 T1 = T2 316 K = JPVT( K ) 317 T2 = B( JPVT( K ), J ) 318 IF( JPVT( K ).NE.I ) 319 $ GO TO 70 320 B( I, J ) = T1 321 WORK( 2*MN+K ) = DONE 322 END IF 323 END IF 324 80 CONTINUE 325 90 CONTINUE 326* 327* Undo scaling 328* 329 IF( IASCL.EQ.1 ) THEN 330 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 331 CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA, 332 $ INFO ) 333 ELSE IF( IASCL.EQ.2 ) THEN 334 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 335 CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA, 336 $ INFO ) 337 END IF 338 IF( IBSCL.EQ.1 ) THEN 339 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 340 ELSE IF( IBSCL.EQ.2 ) THEN 341 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 342 END IF 343* 344 100 CONTINUE 345* 346 RETURN 347* 348* End of DGELSX 349* 350 END 351