LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dgels.f
Go to the documentation of this file.
1 *> \brief <b> DGELS 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 DGELS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgels.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgels.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgels.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGELS( 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 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> DGELS 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 DOUBLE PRECISION 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 DGEQRF;
100 *> if M < N, A is overwritten by details of its LQ
101 *> factorization as returned by DGELQF.
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleGEsolve
181 *
182 * =====================================================================
183  SUBROUTINE dgels( 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  DOUBLE PRECISION a( lda, * ), b( ldb, * ), work( * )
197 * ..
198 *
199 * =====================================================================
200 *
201 * .. Parameters ..
202  DOUBLE PRECISION zero, one
203  parameter( zero = 0.0d0, one = 1.0d0 )
204 * ..
205 * .. Local Scalars ..
206  LOGICAL lquery, tpsd
207  INTEGER brow, i, iascl, ibscl, j, mn, nb, scllen, wsize
208  DOUBLE PRECISION anrm, bignum, bnrm, smlnum
209 * ..
210 * .. Local Arrays ..
211  DOUBLE PRECISION rwork( 1 )
212 * ..
213 * .. External Functions ..
214  LOGICAL lsame
215  INTEGER ilaenv
216  DOUBLE PRECISION dlamch, dlange
217  EXTERNAL lsame, ilaenv, dlabad, dlamch, dlange
218 * ..
219 * .. External Subroutines ..
220  EXTERNAL dgelqf, dgeqrf, dlascl, dlaset, dormlq, dormqr,
221  $ dtrtrs, xerbla
222 * ..
223 * .. Intrinsic Functions ..
224  INTRINSIC dble, max, min
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. .NOT.lquery )
246  $ 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, 'DGEQRF', ' ', m, n, -1, -1 )
260  IF( tpsd ) THEN
261  nb = max( nb, ilaenv( 1, 'DORMQR', 'LN', m, nrhs, n,
262  $ -1 ) )
263  ELSE
264  nb = max( nb, ilaenv( 1, 'DORMQR', 'LT', m, nrhs, n,
265  $ -1 ) )
266  END IF
267  ELSE
268  nb = ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
269  IF( tpsd ) THEN
270  nb = max( nb, ilaenv( 1, 'DORMLQ', 'LT', n, nrhs, m,
271  $ -1 ) )
272  ELSE
273  nb = max( nb, ilaenv( 1, 'DORMLQ', '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 ) = dble( wsize )
280 *
281  END IF
282 *
283  IF( info.NE.0 ) THEN
284  CALL xerbla( 'DGELS ', -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 dlaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
294  return
295  END IF
296 *
297 * Get machine parameters
298 *
299  smlnum = dlamch( 'S' ) / dlamch( 'P' )
300  bignum = one / smlnum
301  CALL dlabad( smlnum, bignum )
302 *
303 * Scale A, B if max element outside range [SMLNUM,BIGNUM]
304 *
305  anrm = dlange( '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 dlascl( '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 dlascl( '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 dlaset( '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 = dlange( '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 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 dgeqrf( 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 dormqr( '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 dtrtrs( '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 dtrtrs( '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 dormqr( '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 dgelqf( 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 dtrtrs( '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 dormlq( '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 dormlq( '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 dtrtrs( '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 dlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
484  $ info )
485  ELSE IF( iascl.EQ.2 ) THEN
486  CALL dlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
487  $ info )
488  END IF
489  IF( ibscl.EQ.1 ) THEN
490  CALL dlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
491  $ info )
492  ELSE IF( ibscl.EQ.2 ) THEN
493  CALL dlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
494  $ info )
495  END IF
496 *
497  50 continue
498  work( 1 ) = dble( wsize )
499 *
500  return
501 *
502 * End of DGELS
503 *
504  END