001:       SUBROUTINE SGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
002:      $                  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:       CHARACTER          TRANS
011:       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               A( LDA, * ), B( LDB, * ), WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  SGELS solves overdetermined or underdetermined real linear systems
021: *  involving an M-by-N matrix A, or its transpose, using a QR or LQ
022: *  factorization of A.  It is assumed that A has full rank.
023: *
024: *  The following options are provided: 
025: *
026: *  1. If TRANS = 'N' and m >= n:  find the least squares solution of
027: *     an overdetermined system, i.e., solve the least squares problem
028: *                  minimize || B - A*X ||.
029: *
030: *  2. If TRANS = 'N' and m < n:  find the minimum norm solution of
031: *     an underdetermined system A * X = B.
032: *
033: *  3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
034: *     an undetermined system A**T * X = B.
035: *
036: *  4. If TRANS = 'T' and m < n:  find the least squares solution of
037: *     an overdetermined system, i.e., solve the least squares problem
038: *                  minimize || B - A**T * X ||.
039: *
040: *  Several right hand side vectors b and solution vectors x can be 
041: *  handled in a single call; they are stored as the columns of the
042: *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution 
043: *  matrix X.
044: *
045: *  Arguments
046: *  =========
047: *
048: *  TRANS   (input) CHARACTER*1
049: *          = 'N': the linear system involves A;
050: *          = 'T': the linear system involves A**T. 
051: *
052: *  M       (input) INTEGER
053: *          The number of rows of the matrix A.  M >= 0.
054: *
055: *  N       (input) INTEGER
056: *          The number of columns of the matrix A.  N >= 0.
057: *
058: *  NRHS    (input) INTEGER
059: *          The number of right hand sides, i.e., the number of
060: *          columns of the matrices B and X. NRHS >=0.
061: *
062: *  A       (input/output) REAL array, dimension (LDA,N)
063: *          On entry, the M-by-N matrix A.
064: *          On exit,
065: *            if M >= N, A is overwritten by details of its QR
066: *                       factorization as returned by SGEQRF;
067: *            if M <  N, A is overwritten by details of its LQ
068: *                       factorization as returned by SGELQF.
069: *
070: *  LDA     (input) INTEGER
071: *          The leading dimension of the array A.  LDA >= max(1,M).
072: *
073: *  B       (input/output) REAL array, dimension (LDB,NRHS)
074: *          On entry, the matrix B of right hand side vectors, stored
075: *          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
076: *          if TRANS = 'T'.  
077: *          On exit, if INFO = 0, B is overwritten by the solution
078: *          vectors, stored columnwise:
079: *          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
080: *          squares solution vectors; the residual sum of squares for the
081: *          solution in each column is given by the sum of squares of
082: *          elements N+1 to M in that column;
083: *          if TRANS = 'N' and m < n, rows 1 to N of B contain the
084: *          minimum norm solution vectors;
085: *          if TRANS = 'T' and m >= n, rows 1 to M of B contain the
086: *          minimum norm solution vectors;
087: *          if TRANS = 'T' and m < n, rows 1 to M of B contain the
088: *          least squares solution vectors; the residual sum of squares
089: *          for the solution in each column is given by the sum of
090: *          squares of elements M+1 to N in that column.
091: *
092: *  LDB     (input) INTEGER
093: *          The leading dimension of the array B. LDB >= MAX(1,M,N).
094: *
095: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
096: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
097: *
098: *  LWORK   (input) INTEGER
099: *          The dimension of the array WORK.
100: *          LWORK >= max( 1, MN + max( MN, NRHS ) ).
101: *          For optimal performance,
102: *          LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
103: *          where MN = min(M,N) and NB is the optimum block size.
104: *
105: *          If LWORK = -1, then a workspace query is assumed; the routine
106: *          only calculates the optimal size of the WORK array, returns
107: *          this value as the first entry of the WORK array, and no error
108: *          message related to LWORK is issued by XERBLA.
109: *
110: *  INFO    (output) INTEGER
111: *          = 0:  successful exit
112: *          < 0:  if INFO = -i, the i-th argument had an illegal value
113: *          > 0:  if INFO =  i, the i-th diagonal element of the
114: *                triangular factor of A is zero, so that A does not have
115: *                full rank; the least squares solution could not be
116: *                computed.
117: *
118: *  =====================================================================
119: *
120: *     .. Parameters ..
121:       REAL               ZERO, ONE
122:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
123: *     ..
124: *     .. Local Scalars ..
125:       LOGICAL            LQUERY, TPSD
126:       INTEGER            BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
127:       REAL               ANRM, BIGNUM, BNRM, SMLNUM
128: *     ..
129: *     .. Local Arrays ..
130:       REAL               RWORK( 1 )
131: *     ..
132: *     .. External Functions ..
133:       LOGICAL            LSAME
134:       INTEGER            ILAENV
135:       REAL               SLAMCH, SLANGE
136:       EXTERNAL           LSAME, ILAENV, SLAMCH, SLANGE
137: *     ..
138: *     .. External Subroutines ..
139:       EXTERNAL           SGELQF, SGEQRF, SLABAD, SLASCL, SLASET, SORMLQ,
140:      $                   SORMQR, STRTRS, XERBLA
141: *     ..
142: *     .. Intrinsic Functions ..
143:       INTRINSIC          MAX, MIN, REAL
144: *     ..
145: *     .. Executable Statements ..
146: *
147: *     Test the input arguments.
148: *
149:       INFO = 0
150:       MN = MIN( M, N )
151:       LQUERY = ( LWORK.EQ.-1 )
152:       IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'T' ) ) ) THEN
153:          INFO = -1
154:       ELSE IF( M.LT.0 ) THEN
155:          INFO = -2
156:       ELSE IF( N.LT.0 ) THEN
157:          INFO = -3
158:       ELSE IF( NRHS.LT.0 ) THEN
159:          INFO = -4
160:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
161:          INFO = -6
162:       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
163:          INFO = -8
164:       ELSE IF( LWORK.LT.MAX( 1, MN + MAX( MN, NRHS ) ) .AND.
165:      $   .NOT.LQUERY ) THEN
166:          INFO = -10
167:       END IF
168: *
169: *     Figure out optimal block size
170: *
171:       IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
172: *
173:          TPSD = .TRUE.
174:          IF( LSAME( TRANS, 'N' ) )
175:      $      TPSD = .FALSE.
176: *
177:          IF( M.GE.N ) THEN
178:             NB = ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
179:             IF( TPSD ) THEN
180:                NB = MAX( NB, ILAENV( 1, 'SORMQR', 'LN', M, NRHS, N,
181:      $              -1 ) )
182:             ELSE
183:                NB = MAX( NB, ILAENV( 1, 'SORMQR', 'LT', M, NRHS, N,
184:      $              -1 ) )
185:             END IF
186:          ELSE
187:             NB = ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
188:             IF( TPSD ) THEN
189:                NB = MAX( NB, ILAENV( 1, 'SORMLQ', 'LT', N, NRHS, M,
190:      $              -1 ) )
191:             ELSE
192:                NB = MAX( NB, ILAENV( 1, 'SORMLQ', 'LN', N, NRHS, M,
193:      $              -1 ) )
194:             END IF
195:          END IF
196: *
197:          WSIZE = MAX( 1, MN + MAX( MN, NRHS )*NB )
198:          WORK( 1 ) = REAL( WSIZE )
199: *
200:       END IF
201: *
202:       IF( INFO.NE.0 ) THEN
203:          CALL XERBLA( 'SGELS ', -INFO )
204:          RETURN
205:       ELSE IF( LQUERY ) THEN
206:          RETURN
207:       END IF
208: *
209: *     Quick return if possible
210: *
211:       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
212:          CALL SLASET( 'Full', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
213:          RETURN
214:       END IF
215: *
216: *     Get machine parameters
217: *
218:       SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
219:       BIGNUM = ONE / SMLNUM
220:       CALL SLABAD( SMLNUM, BIGNUM )
221: *
222: *     Scale A, B if max element outside range [SMLNUM,BIGNUM]
223: *
224:       ANRM = SLANGE( 'M', M, N, A, LDA, RWORK )
225:       IASCL = 0
226:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
227: *
228: *        Scale matrix norm up to SMLNUM
229: *
230:          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
231:          IASCL = 1
232:       ELSE IF( ANRM.GT.BIGNUM ) THEN
233: *
234: *        Scale matrix norm down to BIGNUM
235: *
236:          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
237:          IASCL = 2
238:       ELSE IF( ANRM.EQ.ZERO ) THEN
239: *
240: *        Matrix all zero. Return zero solution.
241: *
242:          CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
243:          GO TO 50
244:       END IF
245: *
246:       BROW = M
247:       IF( TPSD )
248:      $   BROW = N
249:       BNRM = SLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
250:       IBSCL = 0
251:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
252: *
253: *        Scale matrix norm up to SMLNUM
254: *
255:          CALL SLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
256:      $                INFO )
257:          IBSCL = 1
258:       ELSE IF( BNRM.GT.BIGNUM ) THEN
259: *
260: *        Scale matrix norm down to BIGNUM
261: *
262:          CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
263:      $                INFO )
264:          IBSCL = 2
265:       END IF
266: *
267:       IF( M.GE.N ) THEN
268: *
269: *        compute QR factorization of A
270: *
271:          CALL SGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
272:      $                INFO )
273: *
274: *        workspace at least N, optimally N*NB
275: *
276:          IF( .NOT.TPSD ) THEN
277: *
278: *           Least-Squares Problem min || A * X - B ||
279: *
280: *           B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
281: *
282:             CALL SORMQR( 'Left', 'Transpose', M, NRHS, N, A, LDA,
283:      $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
284:      $                   INFO )
285: *
286: *           workspace at least NRHS, optimally NRHS*NB
287: *
288: *           B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
289: *
290:             CALL STRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS,
291:      $                   A, LDA, B, LDB, INFO )
292: *
293:             IF( INFO.GT.0 ) THEN
294:                RETURN
295:             END IF
296: *
297:             SCLLEN = N
298: *
299:          ELSE
300: *
301: *           Overdetermined system of equations A' * X = B
302: *
303: *           B(1:N,1:NRHS) := inv(R') * B(1:N,1:NRHS)
304: *
305:             CALL STRTRS( 'Upper', 'Transpose', 'Non-unit', N, NRHS,
306:      $                   A, LDA, B, LDB, INFO )
307: *
308:             IF( INFO.GT.0 ) THEN
309:                RETURN
310:             END IF
311: *
312: *           B(N+1:M,1:NRHS) = ZERO
313: *
314:             DO 20 J = 1, NRHS
315:                DO 10 I = N + 1, M
316:                   B( I, J ) = ZERO
317:    10          CONTINUE
318:    20       CONTINUE
319: *
320: *           B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
321: *
322:             CALL SORMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA,
323:      $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
324:      $                   INFO )
325: *
326: *           workspace at least NRHS, optimally NRHS*NB
327: *
328:             SCLLEN = M
329: *
330:          END IF
331: *
332:       ELSE
333: *
334: *        Compute LQ factorization of A
335: *
336:          CALL SGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
337:      $                INFO )
338: *
339: *        workspace at least M, optimally M*NB.
340: *
341:          IF( .NOT.TPSD ) THEN
342: *
343: *           underdetermined system of equations A * X = B
344: *
345: *           B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
346: *
347:             CALL STRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS,
348:      $                   A, LDA, B, LDB, INFO )
349: *
350:             IF( INFO.GT.0 ) THEN
351:                RETURN
352:             END IF
353: *
354: *           B(M+1:N,1:NRHS) = 0
355: *
356:             DO 40 J = 1, NRHS
357:                DO 30 I = M + 1, N
358:                   B( I, J ) = ZERO
359:    30          CONTINUE
360:    40       CONTINUE
361: *
362: *           B(1:N,1:NRHS) := Q(1:N,:)' * B(1:M,1:NRHS)
363: *
364:             CALL SORMLQ( 'Left', 'Transpose', N, NRHS, M, A, LDA,
365:      $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
366:      $                   INFO )
367: *
368: *           workspace at least NRHS, optimally NRHS*NB
369: *
370:             SCLLEN = N
371: *
372:          ELSE
373: *
374: *           overdetermined system min || A' * X - B ||
375: *
376: *           B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
377: *
378:             CALL SORMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA,
379:      $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
380:      $                   INFO )
381: *
382: *           workspace at least NRHS, optimally NRHS*NB
383: *
384: *           B(1:M,1:NRHS) := inv(L') * B(1:M,1:NRHS)
385: *
386:             CALL STRTRS( 'Lower', 'Transpose', 'Non-unit', M, NRHS,
387:      $                   A, LDA, B, LDB, INFO )
388: *
389:             IF( INFO.GT.0 ) THEN
390:                RETURN
391:             END IF
392: *
393:             SCLLEN = M
394: *
395:          END IF
396: *
397:       END IF
398: *
399: *     Undo scaling
400: *
401:       IF( IASCL.EQ.1 ) THEN
402:          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
403:      $                INFO )
404:       ELSE IF( IASCL.EQ.2 ) THEN
405:          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
406:      $                INFO )
407:       END IF
408:       IF( IBSCL.EQ.1 ) THEN
409:          CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
410:      $                INFO )
411:       ELSE IF( IBSCL.EQ.2 ) THEN
412:          CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
413:      $                INFO )
414:       END IF
415: *
416:    50 CONTINUE
417:       WORK( 1 ) = REAL( WSIZE )
418: *
419:       RETURN
420: *
421: *     End of SGELS
422: *
423:       END
424: