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