001:       SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
002:      $                   WORK, LWORK, 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, LWORK, 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: *  DGELSY computes the minimum-norm solution to a real linear least
022: *  squares problem:
023: *      minimize || A * X - B ||
024: *  using a complete orthogonal factorization of A.  A is an M-by-N
025: *  matrix which may be rank-deficient.
026: *
027: *  Several right hand side vectors b and solution vectors x can be
028: *  handled in a single call; they are stored as the columns of the
029: *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
030: *  matrix X.
031: *
032: *  The routine first computes a QR factorization with column pivoting:
033: *      A * P = Q * [ R11 R12 ]
034: *                  [  0  R22 ]
035: *  with R11 defined as the largest leading submatrix whose estimated
036: *  condition number is less than 1/RCOND.  The order of R11, RANK,
037: *  is the effective rank of A.
038: *
039: *  Then, R22 is considered to be negligible, and R12 is annihilated
040: *  by orthogonal transformations from the right, arriving at the
041: *  complete orthogonal factorization:
042: *     A * P = Q * [ T11 0 ] * Z
043: *                 [  0  0 ]
044: *  The minimum-norm solution is then
045: *     X = P * Z' [ inv(T11)*Q1'*B ]
046: *                [        0       ]
047: *  where Q1 consists of the first RANK columns of Q.
048: *
049: *  This routine is basically identical to the original xGELSX except
050: *  three differences:
051: *    o The call to the subroutine xGEQPF has been substituted by the
052: *      the call to the subroutine xGEQP3. This subroutine is a Blas-3
053: *      version of the QR factorization with column pivoting.
054: *    o Matrix B (the right hand side) is updated with Blas-3.
055: *    o The permutation of matrix B (the right hand side) is faster and
056: *      more simple.
057: *
058: *  Arguments
059: *  =========
060: *
061: *  M       (input) INTEGER
062: *          The number of rows of the matrix A.  M >= 0.
063: *
064: *  N       (input) INTEGER
065: *          The number of columns of the matrix A.  N >= 0.
066: *
067: *  NRHS    (input) INTEGER
068: *          The number of right hand sides, i.e., the number of
069: *          columns of matrices B and X. NRHS >= 0.
070: *
071: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
072: *          On entry, the M-by-N matrix A.
073: *          On exit, A has been overwritten by details of its
074: *          complete orthogonal factorization.
075: *
076: *  LDA     (input) INTEGER
077: *          The leading dimension of the array A.  LDA >= max(1,M).
078: *
079: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
080: *          On entry, the M-by-NRHS right hand side matrix B.
081: *          On exit, the N-by-NRHS solution matrix X.
082: *
083: *  LDB     (input) INTEGER
084: *          The leading dimension of the array B. LDB >= max(1,M,N).
085: *
086: *  JPVT    (input/output) INTEGER array, dimension (N)
087: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
088: *          to the front of AP, otherwise column i is a free column.
089: *          On exit, if JPVT(i) = k, then the i-th column of AP
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/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
104: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
105: *
106: *  LWORK   (input) INTEGER
107: *          The dimension of the array WORK.
108: *          The unblocked strategy requires that:
109: *             LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),
110: *          where MN = min( M, N ).
111: *          The block algorithm requires that:
112: *             LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
113: *          where NB is an upper bound on the blocksize returned
114: *          by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
115: *          and DORMRZ.
116: *
117: *          If LWORK = -1, then a workspace query is assumed; the routine
118: *          only calculates the optimal size of the WORK array, returns
119: *          this value as the first entry of the WORK array, and no error
120: *          message related to LWORK is issued by XERBLA.
121: *
122: *  INFO    (output) INTEGER
123: *          = 0: successful exit
124: *          < 0: If INFO = -i, the i-th argument had an illegal value.
125: *
126: *  Further Details
127: *  ===============
128: *
129: *  Based on contributions by
130: *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
131: *    E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
132: *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
133: *
134: *  =====================================================================
135: *
136: *     .. Parameters ..
137:       INTEGER            IMAX, IMIN
138:       PARAMETER          ( IMAX = 1, IMIN = 2 )
139:       DOUBLE PRECISION   ZERO, ONE
140:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
141: *     ..
142: *     .. Local Scalars ..
143:       LOGICAL            LQUERY
144:       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN,
145:      $                   LWKOPT, MN, NB, NB1, NB2, NB3, NB4
146:       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
147:      $                   SMAXPR, SMIN, SMINPR, SMLNUM, WSIZE
148: *     ..
149: *     .. External Functions ..
150:       INTEGER            ILAENV
151:       DOUBLE PRECISION   DLAMCH, DLANGE
152:       EXTERNAL           ILAENV, DLAMCH, DLANGE
153: *     ..
154: *     .. External Subroutines ..
155:       EXTERNAL           DCOPY, DGEQP3, DLABAD, DLAIC1, DLASCL, DLASET,
156:      $                   DORMQR, DORMRZ, DTRSM, DTZRZF, XERBLA
157: *     ..
158: *     .. Intrinsic Functions ..
159:       INTRINSIC          ABS, MAX, MIN
160: *     ..
161: *     .. Executable Statements ..
162: *
163:       MN = MIN( M, N )
164:       ISMIN = MN + 1
165:       ISMAX = 2*MN + 1
166: *
167: *     Test the input arguments.
168: *
169:       INFO = 0
170:       LQUERY = ( LWORK.EQ.-1 )
171:       IF( M.LT.0 ) THEN
172:          INFO = -1
173:       ELSE IF( N.LT.0 ) THEN
174:          INFO = -2
175:       ELSE IF( NRHS.LT.0 ) THEN
176:          INFO = -3
177:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
178:          INFO = -5
179:       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
180:          INFO = -7
181:       END IF
182: *
183: *     Figure out optimal block size
184: *
185:       IF( INFO.EQ.0 ) THEN
186:          IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN
187:             LWKMIN = 1
188:             LWKOPT = 1
189:          ELSE
190:             NB1 = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
191:             NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
192:             NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, NRHS, -1 )
193:             NB4 = ILAENV( 1, 'DORMRQ', ' ', M, N, NRHS, -1 )
194:             NB = MAX( NB1, NB2, NB3, NB4 )
195:             LWKMIN = MN + MAX( 2*MN, N + 1, MN + NRHS )
196:             LWKOPT = MAX( LWKMIN,
197:      $                    MN + 2*N + NB*( N + 1 ), 2*MN + NB*NRHS )
198:          END IF
199:          WORK( 1 ) = LWKOPT
200: *
201:          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
202:             INFO = -12
203:          END IF
204:       END IF
205: *
206:       IF( INFO.NE.0 ) THEN
207:          CALL XERBLA( 'DGELSY', -INFO )
208:          RETURN
209:       ELSE IF( LQUERY ) THEN
210:          RETURN
211:       END IF
212: *
213: *     Quick return if possible
214: *
215:       IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN
216:          RANK = 0
217:          RETURN
218:       END IF
219: *
220: *     Get machine parameters
221: *
222:       SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
223:       BIGNUM = ONE / SMLNUM
224:       CALL DLABAD( SMLNUM, BIGNUM )
225: *
226: *     Scale A, B if max entries outside range [SMLNUM,BIGNUM]
227: *
228:       ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
229:       IASCL = 0
230:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
231: *
232: *        Scale matrix norm up to SMLNUM
233: *
234:          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
235:          IASCL = 1
236:       ELSE IF( ANRM.GT.BIGNUM ) THEN
237: *
238: *        Scale matrix norm down to BIGNUM
239: *
240:          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
241:          IASCL = 2
242:       ELSE IF( ANRM.EQ.ZERO ) THEN
243: *
244: *        Matrix all zero. Return zero solution.
245: *
246:          CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
247:          RANK = 0
248:          GO TO 70
249:       END IF
250: *
251:       BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
252:       IBSCL = 0
253:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
254: *
255: *        Scale matrix norm up to SMLNUM
256: *
257:          CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
258:          IBSCL = 1
259:       ELSE IF( BNRM.GT.BIGNUM ) THEN
260: *
261: *        Scale matrix norm down to BIGNUM
262: *
263:          CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
264:          IBSCL = 2
265:       END IF
266: *
267: *     Compute QR factorization with column pivoting of A:
268: *        A * P = Q * R
269: *
270:       CALL DGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ),
271:      $             LWORK-MN, INFO )
272:       WSIZE = MN + WORK( MN+1 )
273: *
274: *     workspace: MN+2*N+NB*(N+1).
275: *     Details of Householder rotations stored in WORK(1:MN).
276: *
277: *     Determine RANK using incremental condition estimation
278: *
279:       WORK( ISMIN ) = ONE
280:       WORK( ISMAX ) = ONE
281:       SMAX = ABS( A( 1, 1 ) )
282:       SMIN = SMAX
283:       IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
284:          RANK = 0
285:          CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
286:          GO TO 70
287:       ELSE
288:          RANK = 1
289:       END IF
290: *
291:    10 CONTINUE
292:       IF( RANK.LT.MN ) THEN
293:          I = RANK + 1
294:          CALL DLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
295:      $                A( I, I ), SMINPR, S1, C1 )
296:          CALL DLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
297:      $                A( I, I ), SMAXPR, S2, C2 )
298: *
299:          IF( SMAXPR*RCOND.LE.SMINPR ) THEN
300:             DO 20 I = 1, RANK
301:                WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
302:                WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
303:    20       CONTINUE
304:             WORK( ISMIN+RANK ) = C1
305:             WORK( ISMAX+RANK ) = C2
306:             SMIN = SMINPR
307:             SMAX = SMAXPR
308:             RANK = RANK + 1
309:             GO TO 10
310:          END IF
311:       END IF
312: *
313: *     workspace: 3*MN.
314: *
315: *     Logically partition R = [ R11 R12 ]
316: *                             [  0  R22 ]
317: *     where R11 = R(1:RANK,1:RANK)
318: *
319: *     [R11,R12] = [ T11, 0 ] * Y
320: *
321:       IF( RANK.LT.N )
322:      $   CALL DTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ),
323:      $                LWORK-2*MN, INFO )
324: *
325: *     workspace: 2*MN.
326: *     Details of Householder rotations stored in WORK(MN+1:2*MN)
327: *
328: *     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
329: *
330:       CALL DORMQR( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ),
331:      $             B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO )
332:       WSIZE = MAX( WSIZE, 2*MN+WORK( 2*MN+1 ) )
333: *
334: *     workspace: 2*MN+NB*NRHS.
335: *
336: *     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
337: *
338:       CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
339:      $            NRHS, ONE, A, LDA, B, LDB )
340: *
341:       DO 40 J = 1, NRHS
342:          DO 30 I = RANK + 1, N
343:             B( I, J ) = ZERO
344:    30    CONTINUE
345:    40 CONTINUE
346: *
347: *     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS)
348: *
349:       IF( RANK.LT.N ) THEN
350:          CALL DORMRZ( 'Left', 'Transpose', N, NRHS, RANK, N-RANK, A,
351:      $                LDA, WORK( MN+1 ), B, LDB, WORK( 2*MN+1 ),
352:      $                LWORK-2*MN, INFO )
353:       END IF
354: *
355: *     workspace: 2*MN+NRHS.
356: *
357: *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
358: *
359:       DO 60 J = 1, NRHS
360:          DO 50 I = 1, N
361:             WORK( JPVT( I ) ) = B( I, J )
362:    50    CONTINUE
363:          CALL DCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 )
364:    60 CONTINUE
365: *
366: *     workspace: N.
367: *
368: *     Undo scaling
369: *
370:       IF( IASCL.EQ.1 ) THEN
371:          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
372:          CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
373:      $                INFO )
374:       ELSE IF( IASCL.EQ.2 ) THEN
375:          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
376:          CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
377:      $                INFO )
378:       END IF
379:       IF( IBSCL.EQ.1 ) THEN
380:          CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
381:       ELSE IF( IBSCL.EQ.2 ) THEN
382:          CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
383:       END IF
384: *
385:    70 CONTINUE
386:       WORK( 1 ) = LWKOPT
387: *
388:       RETURN
389: *
390: *     End of DGELSY
391: *
392:       END
393: