LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dgels()

subroutine dgels ( character  TRANS,
integer  M,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGELS solves overdetermined or underdetermined systems for GE matrices

Download DGELS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DGELS solves overdetermined or underdetermined real linear systems
 involving an M-by-N matrix A, or its transpose, using a QR or LQ
 factorization of A.  It is assumed that A has full rank.

 The following options are provided:

 1. If TRANS = 'N' and m >= n:  find the least squares solution of
    an overdetermined system, i.e., solve the least squares problem
                 minimize || B - A*X ||.

 2. If TRANS = 'N' and m < n:  find the minimum norm solution of
    an underdetermined system A * X = B.

 3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
    an underdetermined system A**T * X = B.

 4. If TRANS = 'T' and m < n:  find the least squares solution of
    an overdetermined system, i.e., solve the least squares problem
                 minimize || B - A**T * X ||.

 Several right hand side vectors b and solution vectors x can be
 handled in a single call; they are stored as the columns of the
 M-by-NRHS right hand side matrix B and the N-by-NRHS solution
 matrix X.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
          = 'N': the linear system involves A;
          = 'T': the linear system involves A**T.
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of
          columns of the matrices B and X. NRHS >=0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit,
            if M >= N, A is overwritten by details of its QR
                       factorization as returned by DGEQRF;
            if M <  N, A is overwritten by details of its LQ
                       factorization as returned by DGELQF.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[in,out]B
          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
          On entry, the matrix B of right hand side vectors, stored
          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
          if TRANS = 'T'.
          On exit, if INFO = 0, B is overwritten by the solution
          vectors, stored columnwise:
          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
          squares solution vectors; the residual sum of squares for the
          solution in each column is given by the sum of squares of
          elements N+1 to M in that column;
          if TRANS = 'N' and m < n, rows 1 to N of B contain the
          minimum norm solution vectors;
          if TRANS = 'T' and m >= n, rows 1 to M of B contain the
          minimum norm solution vectors;
          if TRANS = 'T' and m < n, rows 1 to M of B contain the
          least squares solution vectors; the residual sum of squares
          for the solution in each column is given by the sum of
          squares of elements M+1 to N in that column.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= MAX(1,M,N).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          LWORK >= max( 1, MN + max( MN, NRHS ) ).
          For optimal performance,
          LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
          where MN = min(M,N) and NB is the optimum block size.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO =  i, the i-th diagonal element of the
                triangular factor of A is zero, so that A does not have
                full rank; the least squares solution could not be
                computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 185 of file dgels.f.

185 *
186 * -- LAPACK driver routine (version 3.7.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 * December 2016
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 * Underdetermined 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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:145
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
Definition: dgelqf.f:137
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:169
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:116
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
Definition: dormlq.f:169
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:138
subroutine dtrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
DTRTRS
Definition: dtrtrs.f:142
Here is the call graph for this function:
Here is the caller graph for this function: