LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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, 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*
299* Scale A, B if max element outside range [SMLNUM,BIGNUM]
300*
301 anrm = dlange( 'M', m, n, a, lda, rwork )
302 iascl = 0
303 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
304*
305* Scale matrix norm up to SMLNUM
306*
307 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
308 iascl = 1
309 ELSE IF( anrm.GT.bignum ) THEN
310*
311* Scale matrix norm down to BIGNUM
312*
313 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
314 iascl = 2
315 ELSE IF( anrm.EQ.zero ) THEN
316*
317* Matrix all zero. Return zero solution.
318*
319 CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
320 GO TO 50
321 END IF
322*
323 brow = m
324 IF( tpsd )
325 $ brow = n
326 bnrm = dlange( 'M', brow, nrhs, b, ldb, rwork )
327 ibscl = 0
328 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
329*
330* Scale matrix norm up to SMLNUM
331*
332 CALL dlascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
333 $ info )
334 ibscl = 1
335 ELSE IF( bnrm.GT.bignum ) THEN
336*
337* Scale matrix norm down to BIGNUM
338*
339 CALL dlascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
340 $ info )
341 ibscl = 2
342 END IF
343*
344 IF( m.GE.n ) THEN
345*
346* compute QR factorization of A
347*
348 CALL dgeqrf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
349 $ info )
350*
351* workspace at least N, optimally N*NB
352*
353 IF( .NOT.tpsd ) THEN
354*
355* Least-Squares Problem min || A * X - B ||
356*
357* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
358*
359 CALL dormqr( 'Left', 'Transpose', m, nrhs, n, a, lda,
360 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
361 $ info )
362*
363* workspace at least NRHS, optimally NRHS*NB
364*
365* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
366*
367 CALL dtrtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
368 $ a, lda, b, ldb, info )
369*
370 IF( info.GT.0 ) THEN
371 RETURN
372 END IF
373*
374 scllen = n
375*
376 ELSE
377*
378* Underdetermined system of equations A**T * X = B
379*
380* B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
381*
382 CALL dtrtrs( 'Upper', 'Transpose', 'Non-unit', n, nrhs,
383 $ a, lda, b, ldb, info )
384*
385 IF( info.GT.0 ) THEN
386 RETURN
387 END IF
388*
389* B(N+1:M,1:NRHS) = ZERO
390*
391 DO 20 j = 1, nrhs
392 DO 10 i = n + 1, m
393 b( i, j ) = zero
394 10 CONTINUE
395 20 CONTINUE
396*
397* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
398*
399 CALL dormqr( 'Left', 'No transpose', m, nrhs, n, a, lda,
400 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
401 $ info )
402*
403* workspace at least NRHS, optimally NRHS*NB
404*
405 scllen = m
406*
407 END IF
408*
409 ELSE
410*
411* Compute LQ factorization of A
412*
413 CALL dgelqf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
414 $ info )
415*
416* workspace at least M, optimally M*NB.
417*
418 IF( .NOT.tpsd ) THEN
419*
420* underdetermined system of equations A * X = B
421*
422* B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
423*
424 CALL dtrtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
425 $ a, lda, b, ldb, info )
426*
427 IF( info.GT.0 ) THEN
428 RETURN
429 END IF
430*
431* B(M+1:N,1:NRHS) = 0
432*
433 DO 40 j = 1, nrhs
434 DO 30 i = m + 1, n
435 b( i, j ) = zero
436 30 CONTINUE
437 40 CONTINUE
438*
439* B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
440*
441 CALL dormlq( 'Left', 'Transpose', n, nrhs, m, a, lda,
442 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
443 $ info )
444*
445* workspace at least NRHS, optimally NRHS*NB
446*
447 scllen = n
448*
449 ELSE
450*
451* overdetermined system min || A**T * X - B ||
452*
453* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
454*
455 CALL dormlq( 'Left', 'No transpose', n, nrhs, m, a, lda,
456 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
457 $ info )
458*
459* workspace at least NRHS, optimally NRHS*NB
460*
461* B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
462*
463 CALL dtrtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
464 $ a, lda, b, ldb, info )
465*
466 IF( info.GT.0 ) THEN
467 RETURN
468 END IF
469*
470 scllen = m
471*
472 END IF
473*
474 END IF
475*
476* Undo scaling
477*
478 IF( iascl.EQ.1 ) THEN
479 CALL dlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
480 $ info )
481 ELSE IF( iascl.EQ.2 ) THEN
482 CALL dlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
483 $ info )
484 END IF
485 IF( ibscl.EQ.1 ) THEN
486 CALL dlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
487 $ info )
488 ELSE IF( ibscl.EQ.2 ) THEN
489 CALL dlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
490 $ info )
491 END IF
492*
493 50 CONTINUE
494 work( 1 ) = dble( wsize )
495*
496 RETURN
497*
498* End of DGELS
499*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgelqf(m, n, a, lda, tau, work, lwork, info)
DGELQF
Definition dgelqf.f:143
subroutine dgeqrf(m, n, a, lda, tau, work, lwork, info)
DGEQRF
Definition dgeqrf.f:146
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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 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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:167
Here is the call graph for this function:
Here is the caller graph for this function: