1*> \brief \b DGETSLS
2*
3*  Definition:
4*  ===========
5*
6*       SUBROUTINE DGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
7*     $                     WORK, LWORK, INFO )
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*
18*> \par Purpose:
19*  =============
20*>
21*> \verbatim
22*>
23*> DGETSLS solves overdetermined or underdetermined real linear systems
24*> involving an M-by-N matrix A, using a tall skinny QR or short wide LQ
25*> factorization of A.  It is assumed that A has full rank.
26*>
27*>
28*>
29*> The following options are provided:
30*>
31*> 1. If TRANS = 'N' and m >= n:  find the least squares solution of
32*>    an overdetermined system, i.e., solve the least squares problem
33*>                 minimize || B - A*X ||.
34*>
35*> 2. If TRANS = 'N' and m < n:  find the minimum norm solution of
36*>    an underdetermined system A * X = B.
37*>
38*> 3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
39*>    an undetermined system A**T * X = B.
40*>
41*> 4. If TRANS = 'T' and m < n:  find the least squares solution of
42*>    an overdetermined system, i.e., solve the least squares problem
43*>                 minimize || B - A**T * X ||.
44*>
45*> Several right hand side vectors b and solution vectors x can be
46*> handled in a single call; they are stored as the columns of the
47*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
48*> matrix X.
49*> \endverbatim
50*
51*  Arguments:
52*  ==========
53*
54*> \param[in] TRANS
55*> \verbatim
56*>          TRANS is CHARACTER*1
57*>          = 'N': the linear system involves A;
58*>          = 'T': the linear system involves A**T.
59*> \endverbatim
60*>
61*> \param[in] M
62*> \verbatim
63*>          M is INTEGER
64*>          The number of rows of the matrix A.  M >= 0.
65*> \endverbatim
66*>
67*> \param[in] N
68*> \verbatim
69*>          N is INTEGER
70*>          The number of columns of the matrix A.  N >= 0.
71*> \endverbatim
72*>
73*> \param[in] NRHS
74*> \verbatim
75*>          NRHS is INTEGER
76*>          The number of right hand sides, i.e., the number of
77*>          columns of the matrices B and X. NRHS >=0.
78*> \endverbatim
79*>
80*> \param[in,out] A
81*> \verbatim
82*>          A is DOUBLE PRECISION array, dimension (LDA,N)
83*>          On entry, the M-by-N matrix A.
84*>          On exit,
85*>          A is overwritten by details of its QR or LQ
86*>          factorization as returned by DGEQR or DGELQ.
87*> \endverbatim
88*>
89*> \param[in] LDA
90*> \verbatim
91*>          LDA is INTEGER
92*>          The leading dimension of the array A.  LDA >= max(1,M).
93*> \endverbatim
94*>
95*> \param[in,out] B
96*> \verbatim
97*>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
98*>          On entry, the matrix B of right hand side vectors, stored
99*>          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
100*>          if TRANS = 'T'.
101*>          On exit, if INFO = 0, B is overwritten by the solution
102*>          vectors, stored columnwise:
103*>          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
104*>          squares solution vectors.
105*>          if TRANS = 'N' and m < n, rows 1 to N of B contain the
106*>          minimum norm solution vectors;
107*>          if TRANS = 'T' and m >= n, rows 1 to M of B contain the
108*>          minimum norm solution vectors;
109*>          if TRANS = 'T' and m < n, rows 1 to M of B contain the
110*>          least squares solution vectors.
111*> \endverbatim
112*>
113*> \param[in] LDB
114*> \verbatim
115*>          LDB is INTEGER
116*>          The leading dimension of the array B. LDB >= MAX(1,M,N).
117*> \endverbatim
118*>
119*> \param[out] WORK
120*> \verbatim
121*>          (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
122*>          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
123*>          or optimal, if query was assumed) LWORK.
124*>          See LWORK for details.
125*> \endverbatim
126*>
127*> \param[in] LWORK
128*> \verbatim
129*>          LWORK is INTEGER
130*>          The dimension of the array WORK.
131*>          If LWORK = -1 or -2, then a workspace query is assumed.
132*>          If LWORK = -1, the routine calculates optimal size of WORK for the
133*>          optimal performance and returns this value in WORK(1).
134*>          If LWORK = -2, the routine calculates minimal size of WORK and
135*>          returns this value in WORK(1).
136*> \endverbatim
137*>
138*> \param[out] INFO
139*> \verbatim
140*>          INFO is INTEGER
141*>          = 0:  successful exit
142*>          < 0:  if INFO = -i, the i-th argument had an illegal value
143*>          > 0:  if INFO =  i, the i-th diagonal element of the
144*>                triangular factor of A is zero, so that A does not have
145*>                full rank; the least squares solution could not be
146*>                computed.
147*> \endverbatim
148*
149*  Authors:
150*  ========
151*
152*> \author Univ. of Tennessee
153*> \author Univ. of California Berkeley
154*> \author Univ. of Colorado Denver
155*> \author NAG Ltd.
156*
157*> \ingroup doubleGEsolve
158*
159*  =====================================================================
160      SUBROUTINE DGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
161     $                    WORK, LWORK, INFO )
162*
163*  -- LAPACK driver routine --
164*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
165*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166*
167*     .. Scalar Arguments ..
168      CHARACTER          TRANS
169      INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
170*     ..
171*     .. Array Arguments ..
172      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
173*
174*     ..
175*
176*  =====================================================================
177*
178*     .. Parameters ..
179      DOUBLE PRECISION   ZERO, ONE
180      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
181*     ..
182*     .. Local Scalars ..
183      LOGICAL            LQUERY, TRAN
184      INTEGER            I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
185     $                   SCLLEN, MNK, TSZO, TSZM, LWO, LWM, LW1, LW2,
186     $                   WSIZEO, WSIZEM, INFO2
187      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, SMLNUM, TQ( 5 ), WORKQ( 1 )
188*     ..
189*     .. External Functions ..
190      LOGICAL            LSAME
191      INTEGER            ILAENV
192      DOUBLE PRECISION   DLAMCH, DLANGE
193      EXTERNAL           LSAME, ILAENV, DLABAD, DLAMCH, DLANGE
194*     ..
195*     .. External Subroutines ..
196      EXTERNAL           DGEQR, DGEMQR, DLASCL, DLASET,
197     $                   DTRTRS, XERBLA, DGELQ, DGEMLQ
198*     ..
199*     .. Intrinsic Functions ..
200      INTRINSIC          DBLE, MAX, MIN, INT
201*     ..
202*     .. Executable Statements ..
203*
204*     Test the input arguments.
205*
206      INFO = 0
207      MINMN = MIN( M, N )
208      MAXMN = MAX( M, N )
209      MNK   = MAX( MINMN, NRHS )
210      TRAN  = LSAME( TRANS, 'T' )
211*
212      LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
213      IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
214     $    LSAME( TRANS, 'T' ) ) ) THEN
215         INFO = -1
216      ELSE IF( M.LT.0 ) THEN
217         INFO = -2
218      ELSE IF( N.LT.0 ) THEN
219         INFO = -3
220      ELSE IF( NRHS.LT.0 ) THEN
221         INFO = -4
222      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
223         INFO = -6
224      ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
225         INFO = -8
226      END IF
227*
228      IF( INFO.EQ.0 ) THEN
229*
230*     Determine the block size and minimum LWORK
231*
232       IF( M.GE.N ) THEN
233         CALL DGEQR( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
234         TSZO = INT( TQ( 1 ) )
235         LWO  = INT( WORKQ( 1 ) )
236         CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
237     $                TSZO, B, LDB, WORKQ, -1, INFO2 )
238         LWO  = MAX( LWO, INT( WORKQ( 1 ) ) )
239         CALL DGEQR( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
240         TSZM = INT( TQ( 1 ) )
241         LWM  = INT( WORKQ( 1 ) )
242         CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
243     $                TSZM, B, LDB, WORKQ, -1, INFO2 )
244         LWM = MAX( LWM, INT( WORKQ( 1 ) ) )
245         WSIZEO = TSZO + LWO
246         WSIZEM = TSZM + LWM
247       ELSE
248         CALL DGELQ( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
249         TSZO = INT( TQ( 1 ) )
250         LWO  = INT( WORKQ( 1 ) )
251         CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
252     $                TSZO, B, LDB, WORKQ, -1, INFO2 )
253         LWO  = MAX( LWO, INT( WORKQ( 1 ) ) )
254         CALL DGELQ( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
255         TSZM = INT( TQ( 1 ) )
256         LWM  = INT( WORKQ( 1 ) )
257         CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
258     $                TSZM, B, LDB, WORKQ, -1, INFO2 )
259         LWM  = MAX( LWM, INT( WORKQ( 1 ) ) )
260         WSIZEO = TSZO + LWO
261         WSIZEM = TSZM + LWM
262       END IF
263*
264       IF( ( LWORK.LT.WSIZEM ).AND.( .NOT.LQUERY ) ) THEN
265          INFO = -10
266       END IF
267*
268      END IF
269*
270      IF( INFO.NE.0 ) THEN
271        CALL XERBLA( 'DGETSLS', -INFO )
272        WORK( 1 ) = DBLE( WSIZEO )
273        RETURN
274      END IF
275      IF( LQUERY ) THEN
276        IF( LWORK.EQ.-1 ) WORK( 1 ) = REAL( WSIZEO )
277        IF( LWORK.EQ.-2 ) WORK( 1 ) = REAL( WSIZEM )
278        RETURN
279      END IF
280      IF( LWORK.LT.WSIZEO ) THEN
281        LW1 = TSZM
282        LW2 = LWM
283      ELSE
284        LW1 = TSZO
285        LW2 = LWO
286      END IF
287*
288*     Quick return if possible
289*
290      IF( MIN( M, N, NRHS ).EQ.0 ) THEN
291           CALL DLASET( 'FULL', MAX( M, N ), NRHS, ZERO, ZERO,
292     $                  B, LDB )
293           RETURN
294      END IF
295*
296*     Get machine parameters
297*
298       SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
299       BIGNUM = ONE / SMLNUM
300       CALL DLABAD( SMLNUM, BIGNUM )
301*
302*     Scale A, B if max element outside range [SMLNUM,BIGNUM]
303*
304      ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
305      IASCL = 0
306      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
307*
308*        Scale matrix norm up to SMLNUM
309*
310         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
311         IASCL = 1
312      ELSE IF( ANRM.GT.BIGNUM ) THEN
313*
314*        Scale matrix norm down to BIGNUM
315*
316         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
317         IASCL = 2
318      ELSE IF( ANRM.EQ.ZERO ) THEN
319*
320*        Matrix all zero. Return zero solution.
321*
322         CALL DLASET( 'F', MAXMN, NRHS, ZERO, ZERO, B, LDB )
323         GO TO 50
324      END IF
325*
326      BROW = M
327      IF ( TRAN ) THEN
328        BROW = N
329      END IF
330      BNRM = DLANGE( 'M', BROW, NRHS, B, LDB, WORK )
331      IBSCL = 0
332      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
333*
334*        Scale matrix norm up to SMLNUM
335*
336         CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
337     $                INFO )
338         IBSCL = 1
339      ELSE IF( BNRM.GT.BIGNUM ) THEN
340*
341*        Scale matrix norm down to BIGNUM
342*
343         CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
344     $                INFO )
345         IBSCL = 2
346      END IF
347*
348      IF ( M.GE.N ) THEN
349*
350*        compute QR factorization of A
351*
352        CALL DGEQR( M, N, A, LDA, WORK( LW2+1 ), LW1,
353     $              WORK( 1 ), LW2, INFO )
354        IF ( .NOT.TRAN ) THEN
355*
356*           Least-Squares Problem min || A * X - B ||
357*
358*           B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
359*
360          CALL DGEMQR( 'L' , 'T', M, NRHS, N, A, LDA,
361     $                 WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
362     $                 INFO )
363*
364*           B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
365*
366          CALL DTRTRS( 'U', 'N', 'N', N, NRHS,
367     $                  A, LDA, B, LDB, INFO )
368          IF( INFO.GT.0 ) THEN
369            RETURN
370          END IF
371          SCLLEN = N
372        ELSE
373*
374*           Overdetermined system of equations A**T * X = B
375*
376*           B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
377*
378            CALL DTRTRS( 'U', 'T', 'N', N, NRHS,
379     $                   A, LDA, B, LDB, INFO )
380*
381            IF( INFO.GT.0 ) THEN
382               RETURN
383            END IF
384*
385*           B(N+1:M,1:NRHS) = ZERO
386*
387            DO 20 J = 1, NRHS
388               DO 10 I = N + 1, M
389                  B( I, J ) = ZERO
390   10          CONTINUE
391   20       CONTINUE
392*
393*           B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
394*
395            CALL DGEMQR( 'L', 'N', M, NRHS, N, A, LDA,
396     $                   WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
397     $                   INFO )
398*
399            SCLLEN = M
400*
401         END IF
402*
403      ELSE
404*
405*        Compute LQ factorization of A
406*
407         CALL DGELQ( M, N, A, LDA, WORK( LW2+1 ), LW1,
408     $               WORK( 1 ), LW2, INFO )
409*
410*        workspace at least M, optimally M*NB.
411*
412         IF( .NOT.TRAN ) THEN
413*
414*           underdetermined system of equations A * X = B
415*
416*           B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
417*
418            CALL DTRTRS( 'L', 'N', 'N', M, NRHS,
419     $                   A, LDA, B, LDB, INFO )
420*
421            IF( INFO.GT.0 ) THEN
422               RETURN
423            END IF
424*
425*           B(M+1:N,1:NRHS) = 0
426*
427            DO 40 J = 1, NRHS
428               DO 30 I = M + 1, N
429                  B( I, J ) = ZERO
430   30          CONTINUE
431   40       CONTINUE
432*
433*           B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
434*
435            CALL DGEMLQ( 'L', 'T', N, NRHS, M, A, LDA,
436     $                   WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
437     $                   INFO )
438*
439*           workspace at least NRHS, optimally NRHS*NB
440*
441            SCLLEN = N
442*
443         ELSE
444*
445*           overdetermined system min || A**T * X - B ||
446*
447*           B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
448*
449            CALL DGEMLQ( 'L', 'N', N, NRHS, M, A, LDA,
450     $                   WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
451     $                   INFO )
452*
453*           workspace at least NRHS, optimally NRHS*NB
454*
455*           B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
456*
457            CALL DTRTRS( 'Lower', 'Transpose', 'Non-unit', M, NRHS,
458     $                   A, LDA, B, LDB, INFO )
459*
460            IF( INFO.GT.0 ) THEN
461               RETURN
462            END IF
463*
464            SCLLEN = M
465*
466         END IF
467*
468      END IF
469*
470*     Undo scaling
471*
472      IF( IASCL.EQ.1 ) THEN
473        CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
474     $               INFO )
475      ELSE IF( IASCL.EQ.2 ) THEN
476        CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
477     $               INFO )
478      END IF
479      IF( IBSCL.EQ.1 ) THEN
480        CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
481     $               INFO )
482      ELSE IF( IBSCL.EQ.2 ) THEN
483        CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
484     $               INFO )
485      END IF
486*
487   50 CONTINUE
488      WORK( 1 ) = DBLE( TSZO + LWO )
489      RETURN
490*
491*     End of DGETSLS
492*
493      END
494