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