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

◆ dgelst()

subroutine dgelst ( 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 
)

DGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization with compact WY representation of Q.

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

Purpose:
 DGELST 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 with compact WY representation of Q.
 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 DGEQRT;
            if M <  N, A is overwritten by details of its LQ
                       factorization as returned by DGELQT.
[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.
Contributors:
  November 2022,  Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

Definition at line 192 of file dgelst.f.

194*
195* -- LAPACK driver routine --
196* -- LAPACK is a software package provided by Univ. of Tennessee, --
197* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
198*
199* .. Scalar Arguments ..
200 CHARACTER TRANS
201 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
202* ..
203* .. Array Arguments ..
204 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
205* ..
206*
207* =====================================================================
208*
209* .. Parameters ..
210 DOUBLE PRECISION ZERO, ONE
211 parameter( zero = 0.0d+0, one = 1.0d+0 )
212* ..
213* .. Local Scalars ..
214 LOGICAL LQUERY, TPSD
215 INTEGER BROW, I, IASCL, IBSCL, J, LWOPT, MN, MNNRHS,
216 $ NB, NBMIN, SCLLEN
217 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
218* ..
219* .. Local Arrays ..
220 DOUBLE PRECISION RWORK( 1 )
221* ..
222* .. External Functions ..
223 LOGICAL LSAME
224 INTEGER ILAENV
225 DOUBLE PRECISION DLAMCH, DLANGE
226 EXTERNAL lsame, ilaenv, dlamch, dlange
227* ..
228* .. External Subroutines ..
229 EXTERNAL dgelqt, dgeqrt, dgemlqt, dgemqrt, dlascl,
231* ..
232* .. Intrinsic Functions ..
233 INTRINSIC dble, max, min
234* ..
235* .. Executable Statements ..
236*
237* Test the input arguments.
238*
239 info = 0
240 mn = min( m, n )
241 lquery = ( lwork.EQ.-1 )
242 IF( .NOT.( lsame( trans, 'N' ) .OR. lsame( trans, 'T' ) ) ) THEN
243 info = -1
244 ELSE IF( m.LT.0 ) THEN
245 info = -2
246 ELSE IF( n.LT.0 ) THEN
247 info = -3
248 ELSE IF( nrhs.LT.0 ) THEN
249 info = -4
250 ELSE IF( lda.LT.max( 1, m ) ) THEN
251 info = -6
252 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
253 info = -8
254 ELSE IF( lwork.LT.max( 1, mn+max( mn, nrhs ) ) .AND. .NOT.lquery )
255 $ THEN
256 info = -10
257 END IF
258*
259* Figure out optimal block size and optimal workspace size
260*
261 IF( info.EQ.0 .OR. info.EQ.-10 ) THEN
262*
263 tpsd = .true.
264 IF( lsame( trans, 'N' ) )
265 $ tpsd = .false.
266*
267 nb = ilaenv( 1, 'DGELST', ' ', m, n, -1, -1 )
268*
269 mnnrhs = max( mn, nrhs )
270 lwopt = max( 1, (mn+mnnrhs)*nb )
271 work( 1 ) = dble( lwopt )
272*
273 END IF
274*
275 IF( info.NE.0 ) THEN
276 CALL xerbla( 'DGELST ', -info )
277 RETURN
278 ELSE IF( lquery ) THEN
279 RETURN
280 END IF
281*
282* Quick return if possible
283*
284 IF( min( m, n, nrhs ).EQ.0 ) THEN
285 CALL dlaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
286 work( 1 ) = dble( lwopt )
287 RETURN
288 END IF
289*
290* *GEQRT and *GELQT routines cannot accept NB larger than min(M,N)
291*
292 IF( nb.GT.mn ) nb = mn
293*
294* Determine the block size from the supplied LWORK
295* ( at this stage we know that LWORK >= (minimum required workspace,
296* but it may be less than optimal)
297*
298 nb = min( nb, lwork/( mn + mnnrhs ) )
299*
300* The minimum value of NB, when blocked code is used
301*
302 nbmin = max( 2, ilaenv( 2, 'DGELST', ' ', m, n, -1, -1 ) )
303*
304 IF( nb.LT.nbmin ) THEN
305 nb = 1
306 END IF
307*
308* Get machine parameters
309*
310 smlnum = dlamch( 'S' ) / dlamch( 'P' )
311 bignum = one / smlnum
312*
313* Scale A, B if max element outside range [SMLNUM,BIGNUM]
314*
315 anrm = dlange( 'M', m, n, a, lda, rwork )
316 iascl = 0
317 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
318*
319* Scale matrix norm up to SMLNUM
320*
321 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
322 iascl = 1
323 ELSE IF( anrm.GT.bignum ) THEN
324*
325* Scale matrix norm down to BIGNUM
326*
327 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
328 iascl = 2
329 ELSE IF( anrm.EQ.zero ) THEN
330*
331* Matrix all zero. Return zero solution.
332*
333 CALL dlaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
334 work( 1 ) = dble( lwopt )
335 RETURN
336 END IF
337*
338 brow = m
339 IF( tpsd )
340 $ brow = n
341 bnrm = dlange( 'M', brow, nrhs, b, ldb, rwork )
342 ibscl = 0
343 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
344*
345* Scale matrix norm up to SMLNUM
346*
347 CALL dlascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
348 $ info )
349 ibscl = 1
350 ELSE IF( bnrm.GT.bignum ) THEN
351*
352* Scale matrix norm down to BIGNUM
353*
354 CALL dlascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
355 $ info )
356 ibscl = 2
357 END IF
358*
359 IF( m.GE.n ) THEN
360*
361* M > N:
362* Compute the blocked QR factorization of A,
363* using the compact WY representation of Q,
364* workspace at least N, optimally N*NB.
365*
366 CALL dgeqrt( m, n, nb, a, lda, work( 1 ), nb,
367 $ work( mn*nb+1 ), info )
368*
369 IF( .NOT.tpsd ) THEN
370*
371* M > N, A is not transposed:
372* Overdetermined system of equations,
373* least-squares problem, min || A * X - B ||.
374*
375* Compute B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS),
376* using the compact WY representation of Q,
377* workspace at least NRHS, optimally NRHS*NB.
378*
379 CALL dgemqrt( 'Left', 'Transpose', m, nrhs, n, nb, a, lda,
380 $ work( 1 ), nb, b, ldb, work( mn*nb+1 ),
381 $ info )
382*
383* Compute B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
384*
385 CALL dtrtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
386 $ a, lda, b, ldb, info )
387*
388 IF( info.GT.0 ) THEN
389 RETURN
390 END IF
391*
392 scllen = n
393*
394 ELSE
395*
396* M > N, A is transposed:
397* Underdetermined system of equations,
398* minimum norm solution of A**T * X = B.
399*
400* Compute B := inv(R**T) * B in two row blocks of B.
401*
402* Block 1: B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
403*
404 CALL dtrtrs( 'Upper', 'Transpose', 'Non-unit', n, nrhs,
405 $ a, lda, b, ldb, info )
406*
407 IF( info.GT.0 ) THEN
408 RETURN
409 END IF
410*
411* Block 2: Zero out all rows below the N-th row in B:
412* B(N+1:M,1:NRHS) = ZERO
413*
414 DO j = 1, nrhs
415 DO i = n + 1, m
416 b( i, j ) = zero
417 END DO
418 END DO
419*
420* Compute B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS),
421* using the compact WY representation of Q,
422* workspace at least NRHS, optimally NRHS*NB.
423*
424 CALL dgemqrt( 'Left', 'No transpose', m, nrhs, n, nb,
425 $ a, lda, work( 1 ), nb, b, ldb,
426 $ work( mn*nb+1 ), info )
427*
428 scllen = m
429*
430 END IF
431*
432 ELSE
433*
434* M < N:
435* Compute the blocked LQ factorization of A,
436* using the compact WY representation of Q,
437* workspace at least M, optimally M*NB.
438*
439 CALL dgelqt( m, n, nb, a, lda, work( 1 ), nb,
440 $ work( mn*nb+1 ), info )
441*
442 IF( .NOT.tpsd ) THEN
443*
444* M < N, A is not transposed:
445* Underdetermined system of equations,
446* minimum norm solution of A * X = B.
447*
448* Compute B := inv(L) * B in two row blocks of B.
449*
450* Block 1: B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
451*
452 CALL dtrtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
453 $ a, lda, b, ldb, info )
454*
455 IF( info.GT.0 ) THEN
456 RETURN
457 END IF
458*
459* Block 2: Zero out all rows below the M-th row in B:
460* B(M+1:N,1:NRHS) = ZERO
461*
462 DO j = 1, nrhs
463 DO i = m + 1, n
464 b( i, j ) = zero
465 END DO
466 END DO
467*
468* Compute B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS),
469* using the compact WY representation of Q,
470* workspace at least NRHS, optimally NRHS*NB.
471*
472 CALL dgemlqt( 'Left', 'Transpose', n, nrhs, m, nb, a, lda,
473 $ work( 1 ), nb, b, ldb,
474 $ work( mn*nb+1 ), info )
475*
476 scllen = n
477*
478 ELSE
479*
480* M < N, A is transposed:
481* Overdetermined system of equations,
482* least-squares problem, min || A**T * X - B ||.
483*
484* Compute B(1:N,1:NRHS) := Q * B(1:N,1:NRHS),
485* using the compact WY representation of Q,
486* workspace at least NRHS, optimally NRHS*NB.
487*
488 CALL dgemlqt( 'Left', 'No transpose', n, nrhs, m, nb,
489 $ a, lda, work( 1 ), nb, b, ldb,
490 $ work( mn*nb+1), info )
491*
492* Compute B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
493*
494 CALL dtrtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
495 $ a, lda, b, ldb, info )
496*
497 IF( info.GT.0 ) THEN
498 RETURN
499 END IF
500*
501 scllen = m
502*
503 END IF
504*
505 END IF
506*
507* Undo scaling
508*
509 IF( iascl.EQ.1 ) THEN
510 CALL dlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
511 $ info )
512 ELSE IF( iascl.EQ.2 ) THEN
513 CALL dlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
514 $ info )
515 END IF
516 IF( ibscl.EQ.1 ) THEN
517 CALL dlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
518 $ info )
519 ELSE IF( ibscl.EQ.2 ) THEN
520 CALL dlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
521 $ info )
522 END IF
523*
524 work( 1 ) = dble( lwopt )
525*
526 RETURN
527*
528* End of DGELST
529*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgelqt(m, n, mb, a, lda, t, ldt, work, info)
DGELQT
Definition dgelqt.f:139
subroutine dgemlqt(side, trans, m, n, k, mb, v, ldv, t, ldt, c, ldc, work, info)
DGEMLQT
Definition dgemlqt.f:168
subroutine dgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
DGEMQRT
Definition dgemqrt.f:168
subroutine dgeqrt(m, n, nb, a, lda, t, ldt, work, info)
DGEQRT
Definition dgeqrt.f:141
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
Here is the call graph for this function:
Here is the caller graph for this function: