LAPACK  3.9.1
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.

Definition at line 181 of file dgels.f.

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  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
194 * ..
195 *
196 * =====================================================================
197 *
198 * .. Parameters ..
199  DOUBLE PRECISION ZERO, ONE
200  parameter( zero = 0.0d0, one = 1.0d0 )
201 * ..
202 * .. Local Scalars ..
203  LOGICAL LQUERY, TPSD
204  INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
205  DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
206 * ..
207 * .. Local Arrays ..
208  DOUBLE PRECISION RWORK( 1 )
209 * ..
210 * .. External Functions ..
211  LOGICAL LSAME
212  INTEGER ILAENV
213  DOUBLE PRECISION DLAMCH, DLANGE
214  EXTERNAL lsame, ilaenv, dlabad, dlamch, dlange
215 * ..
216 * .. External Subroutines ..
217  EXTERNAL dgelqf, dgeqrf, dlascl, dlaset, dormlq, dormqr,
218  $ dtrtrs, xerbla
219 * ..
220 * .. Intrinsic Functions ..
221  INTRINSIC dble, max, min
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. .NOT.lquery )
243  $ 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, 'DGEQRF', ' ', m, n, -1, -1 )
257  IF( tpsd ) THEN
258  nb = max( nb, ilaenv( 1, 'DORMQR', 'LN', m, nrhs, n,
259  $ -1 ) )
260  ELSE
261  nb = max( nb, ilaenv( 1, 'DORMQR', 'LT', m, nrhs, n,
262  $ -1 ) )
263  END IF
264  ELSE
265  nb = ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
266  IF( tpsd ) THEN
267  nb = max( nb, ilaenv( 1, 'DORMLQ', 'LT', n, nrhs, m,
268  $ -1 ) )
269  ELSE
270  nb = max( nb, ilaenv( 1, 'DORMLQ', '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 ) = dble( wsize )
277 *
278  END IF
279 *
280  IF( info.NE.0 ) THEN
281  CALL xerbla( 'DGELS ', -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 dlaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
291  RETURN
292  END IF
293 *
294 * Get machine parameters
295 *
296  smlnum = dlamch( 'S' ) / dlamch( 'P' )
297  bignum = one / smlnum
298  CALL dlabad( smlnum, bignum )
299 *
300 * Scale A, B if max element outside range [SMLNUM,BIGNUM]
301 *
302  anrm = dlange( '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 dlascl( '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 dlascl( '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 dlaset( '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 = dlange( '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 dlascl( '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 dlascl( '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 dgeqrf( 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 dormqr( '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 dtrtrs( '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 dtrtrs( '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 dormqr( '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 dgelqf( 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 dtrtrs( '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 dormlq( '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 dormlq( '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 dtrtrs( '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 dlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
481  $ info )
482  ELSE IF( iascl.EQ.2 ) THEN
483  CALL dlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
484  $ info )
485  END IF
486  IF( ibscl.EQ.1 ) THEN
487  CALL dlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
488  $ info )
489  ELSE IF( ibscl.EQ.2 ) THEN
490  CALL dlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
491  $ info )
492  END IF
493 *
494  50 CONTINUE
495  work( 1 ) = dble( wsize )
496 *
497  RETURN
498 *
499 * End of DGELS
500 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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:143
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:110
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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:114
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:145
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
Definition: dgelqf.f:143
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:167
subroutine dtrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
DTRTRS
Definition: dtrtrs.f:140
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
Definition: dormlq.f:167
Here is the call graph for this function:
Here is the caller graph for this function: