001:       SUBROUTINE DGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
002:      $                   WORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            INFO, LDA, LDB, M, N, NRHS, RANK
011:       DOUBLE PRECISION   RCOND
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            JPVT( * )
015:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  This routine is deprecated and has been replaced by routine DGELSY.
022: *
023: *  DGELSX computes the minimum-norm solution to a real linear least
024: *  squares problem:
025: *      minimize || A * X - B ||
026: *  using a complete orthogonal factorization of A.  A is an M-by-N
027: *  matrix which may be rank-deficient.
028: *
029: *  Several right hand side vectors b and solution vectors x can be
030: *  handled in a single call; they are stored as the columns of the
031: *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
032: *  matrix X.
033: *
034: *  The routine first computes a QR factorization with column pivoting:
035: *      A * P = Q * [ R11 R12 ]
036: *                  [  0  R22 ]
037: *  with R11 defined as the largest leading submatrix whose estimated
038: *  condition number is less than 1/RCOND.  The order of R11, RANK,
039: *  is the effective rank of A.
040: *
041: *  Then, R22 is considered to be negligible, and R12 is annihilated
042: *  by orthogonal transformations from the right, arriving at the
043: *  complete orthogonal factorization:
044: *     A * P = Q * [ T11 0 ] * Z
045: *                 [  0  0 ]
046: *  The minimum-norm solution is then
047: *     X = P * Z' [ inv(T11)*Q1'*B ]
048: *                [        0       ]
049: *  where Q1 consists of the first RANK columns of Q.
050: *
051: *  Arguments
052: *  =========
053: *
054: *  M       (input) INTEGER
055: *          The number of rows of the matrix A.  M >= 0.
056: *
057: *  N       (input) INTEGER
058: *          The number of columns of the matrix A.  N >= 0.
059: *
060: *  NRHS    (input) INTEGER
061: *          The number of right hand sides, i.e., the number of
062: *          columns of matrices B and X. NRHS >= 0.
063: *
064: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
065: *          On entry, the M-by-N matrix A.
066: *          On exit, A has been overwritten by details of its
067: *          complete orthogonal factorization.
068: *
069: *  LDA     (input) INTEGER
070: *          The leading dimension of the array A.  LDA >= max(1,M).
071: *
072: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
073: *          On entry, the M-by-NRHS right hand side matrix B.
074: *          On exit, the N-by-NRHS solution matrix X.
075: *          If m >= n and RANK = n, the residual sum-of-squares for
076: *          the solution in the i-th column is given by the sum of
077: *          squares of elements N+1:M in that column.
078: *
079: *  LDB     (input) INTEGER
080: *          The leading dimension of the array B. LDB >= max(1,M,N).
081: *
082: *  JPVT    (input/output) INTEGER array, dimension (N)
083: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is an
084: *          initial column, otherwise it is a free column.  Before
085: *          the QR factorization of A, all initial columns are
086: *          permuted to the leading positions; only the remaining
087: *          free columns are moved as a result of column pivoting
088: *          during the factorization.
089: *          On exit, if JPVT(i) = k, then the i-th column of A*P
090: *          was the k-th column of A.
091: *
092: *  RCOND   (input) DOUBLE PRECISION
093: *          RCOND is used to determine the effective rank of A, which
094: *          is defined as the order of the largest leading triangular
095: *          submatrix R11 in the QR factorization with pivoting of A,
096: *          whose estimated condition number < 1/RCOND.
097: *
098: *  RANK    (output) INTEGER
099: *          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: