LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
sgels.f
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 underdetermined 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 *> \ingroup gels
179 *
180 * =====================================================================
181  SUBROUTINE sgels( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
182  $ INFO )
183 *
184 * -- LAPACK driver routine --
185 * -- LAPACK is a software package provided by Univ. of Tennessee, --
186 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187 *
188 * .. Scalar Arguments ..
189  CHARACTER TRANS
190  INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
191 * ..
192 * .. Array Arguments ..
193  REAL A( lda, * ), B( ldb, * ), WORK( * )
194 * ..
195 *
196 * =====================================================================
197 *
198 * .. Parameters ..
199  REAL ZERO, ONE
200  parameter( zero = 0.0e0, one = 1.0e0 )
201 * ..
202 * .. Local Scalars ..
203  LOGICAL LQUERY, TPSD
204  INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
205  REAL ANRM, BIGNUM, BNRM, SMLNUM
206 * ..
207 * .. Local Arrays ..
208  REAL RWORK( 1 )
209 * ..
210 * .. External Functions ..
211  LOGICAL LSAME
212  INTEGER ILAENV
213  REAL SLAMCH, SLANGE
214  EXTERNAL lsame, ilaenv, slamch, slange
215 * ..
216 * .. External Subroutines ..
217  EXTERNAL sgelqf, sgeqrf, slabad, slascl, slaset, sormlq,
218  $ sormqr, strtrs, xerbla
219 * ..
220 * .. Intrinsic Functions ..
221  INTRINSIC max, min, real
222 * ..
223 * .. Executable Statements ..
224 *
225 * Test the input arguments.
226 *
227  info = 0
228  mn = min( m, n )
229  lquery = ( lwork.EQ.-1 )
230  IF( .NOT.( lsame( trans, 'N' ) .OR. lsame( trans, 'T' ) ) ) THEN
231  info = -1
232  ELSE IF( m.LT.0 ) THEN
233  info = -2
234  ELSE IF( n.LT.0 ) THEN
235  info = -3
236  ELSE IF( nrhs.LT.0 ) THEN
237  info = -4
238  ELSE IF( lda.LT.max( 1, m ) ) THEN
239  info = -6
240  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
241  info = -8
242  ELSE IF( lwork.LT.max( 1, mn + max( mn, nrhs ) ) .AND.
243  $ .NOT.lquery ) THEN
244  info = -10
245  END IF
246 *
247 * Figure out optimal block size
248 *
249  IF( info.EQ.0 .OR. info.EQ.-10 ) THEN
250 *
251  tpsd = .true.
252  IF( lsame( trans, 'N' ) )
253  $ tpsd = .false.
254 *
255  IF( m.GE.n ) THEN
256  nb = ilaenv( 1, 'SGEQRF', ' ', m, n, -1, -1 )
257  IF( tpsd ) THEN
258  nb = max( nb, ilaenv( 1, 'SORMQR', 'LN', m, nrhs, n,
259  $ -1 ) )
260  ELSE
261  nb = max( nb, ilaenv( 1, 'SORMQR', 'LT', m, nrhs, n,
262  $ -1 ) )
263  END IF
264  ELSE
265  nb = ilaenv( 1, 'SGELQF', ' ', m, n, -1, -1 )
266  IF( tpsd ) THEN
267  nb = max( nb, ilaenv( 1, 'SORMLQ', 'LT', n, nrhs, m,
268  $ -1 ) )
269  ELSE
270  nb = max( nb, ilaenv( 1, 'SORMLQ', 'LN', n, nrhs, m,
271  $ -1 ) )
272  END IF
273  END IF
274 *
275  wsize = max( 1, mn + max( mn, nrhs )*nb )
276  work( 1 ) = REAL( wsize )
277 *
278  END IF
279 *
280  IF( info.NE.0 ) THEN
281  CALL xerbla( 'SGELS ', -info )
282  RETURN
283  ELSE IF( lquery ) THEN
284  RETURN
285  END IF
286 *
287 * Quick return if possible
288 *
289  IF( min( m, n, nrhs ).EQ.0 ) THEN
290  CALL slaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
291  RETURN
292  END IF
293 *
294 * Get machine parameters
295 *
296  smlnum = slamch( 'S' ) / slamch( 'P' )
297  bignum = one / smlnum
298  CALL slabad( smlnum, bignum )
299 *
300 * Scale A, B if max element outside range [SMLNUM,BIGNUM]
301 *
302  anrm = slange( 'M', m, n, a, lda, rwork )
303  iascl = 0
304  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
305 *
306 * Scale matrix norm up to SMLNUM
307 *
308  CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
309  iascl = 1
310  ELSE IF( anrm.GT.bignum ) THEN
311 *
312 * Scale matrix norm down to BIGNUM
313 *
314  CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
315  iascl = 2
316  ELSE IF( anrm.EQ.zero ) THEN
317 *
318 * Matrix all zero. Return zero solution.
319 *
320  CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
321  GO TO 50
322  END IF
323 *
324  brow = m
325  IF( tpsd )
326  $ brow = n
327  bnrm = slange( 'M', brow, nrhs, b, ldb, rwork )
328  ibscl = 0
329  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
330 *
331 * Scale matrix norm up to SMLNUM
332 *
333  CALL slascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
334  $ info )
335  ibscl = 1
336  ELSE IF( bnrm.GT.bignum ) THEN
337 *
338 * Scale matrix norm down to BIGNUM
339 *
340  CALL slascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
341  $ info )
342  ibscl = 2
343  END IF
344 *
345  IF( m.GE.n ) THEN
346 *
347 * compute QR factorization of A
348 *
349  CALL sgeqrf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
350  $ info )
351 *
352 * workspace at least N, optimally N*NB
353 *
354  IF( .NOT.tpsd ) 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 sormqr( 'Left', 'Transpose', m, nrhs, n, a, lda,
361  $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
362  $ info )
363 *
364 * workspace at least NRHS, optimally NRHS*NB
365 *
366 * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
367 *
368  CALL strtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
369  $ a, lda, b, ldb, info )
370 *
371  IF( info.GT.0 ) THEN
372  RETURN
373  END IF
374 *
375  scllen = n
376 *
377  ELSE
378 *
379 * Underdetermined system of equations A**T * X = B
380 *
381 * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
382 *
383  CALL strtrs( 'Upper', 'Transpose', 'Non-unit', n, nrhs,
384  $ a, lda, b, ldb, info )
385 *
386  IF( info.GT.0 ) THEN
387  RETURN
388  END IF
389 *
390 * B(N+1:M,1:NRHS) = ZERO
391 *
392  DO 20 j = 1, nrhs
393  DO 10 i = n + 1, m
394  b( i, j ) = zero
395  10 CONTINUE
396  20 CONTINUE
397 *
398 * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
399 *
400  CALL sormqr( 'Left', 'No transpose', m, nrhs, n, a, lda,
401  $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
402  $ info )
403 *
404 * workspace at least NRHS, optimally NRHS*NB
405 *
406  scllen = m
407 *
408  END IF
409 *
410  ELSE
411 *
412 * Compute LQ factorization of A
413 *
414  CALL sgelqf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
415  $ info )
416 *
417 * workspace at least M, optimally M*NB.
418 *
419  IF( .NOT.tpsd ) THEN
420 *
421 * underdetermined system of equations A * X = B
422 *
423 * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
424 *
425  CALL strtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
426  $ a, lda, b, ldb, info )
427 *
428  IF( info.GT.0 ) THEN
429  RETURN
430  END IF
431 *
432 * B(M+1:N,1:NRHS) = 0
433 *
434  DO 40 j = 1, nrhs
435  DO 30 i = m + 1, n
436  b( i, j ) = zero
437  30 CONTINUE
438  40 CONTINUE
439 *
440 * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
441 *
442  CALL sormlq( 'Left', 'Transpose', n, nrhs, m, a, lda,
443  $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
444  $ info )
445 *
446 * workspace at least NRHS, optimally NRHS*NB
447 *
448  scllen = n
449 *
450  ELSE
451 *
452 * overdetermined system min || A**T * X - B ||
453 *
454 * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
455 *
456  CALL sormlq( 'Left', 'No transpose', n, nrhs, m, a, lda,
457  $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
458  $ info )
459 *
460 * workspace at least NRHS, optimally NRHS*NB
461 *
462 * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
463 *
464  CALL strtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
465  $ a, lda, b, ldb, info )
466 *
467  IF( info.GT.0 ) THEN
468  RETURN
469  END IF
470 *
471  scllen = m
472 *
473  END IF
474 *
475  END IF
476 *
477 * Undo scaling
478 *
479  IF( iascl.EQ.1 ) THEN
480  CALL slascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
481  $ info )
482  ELSE IF( iascl.EQ.2 ) THEN
483  CALL slascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
484  $ info )
485  END IF
486  IF( ibscl.EQ.1 ) THEN
487  CALL slascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
488  $ info )
489  ELSE IF( ibscl.EQ.2 ) THEN
490  CALL slascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
491  $ info )
492  END IF
493 *
494  50 CONTINUE
495  work( 1 ) = REAL( wsize )
496 *
497  RETURN
498 *
499 * End of SGELS
500 *
501  END
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
sgels
subroutine sgels(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
SGELS solves overdetermined or underdetermined systems for GE matrices
Definition: sgels.f:183
slascl
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
sormqr
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:168
slaset
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:110
sormlq
subroutine sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
Definition: sormlq.f:168
sgelqf
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
Definition: sgelqf.f:143
strtrs
subroutine strtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
STRTRS
Definition: strtrs.f:140
sgeqrf
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:146