LAPACK 3.11.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, 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:146
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: