LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ dgetsls()

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

DGETSLS

Purpose:
 DGETSLS solves overdetermined or underdetermined real linear systems
 involving an M-by-N matrix A, using a tall skinny QR or short wide 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 undetermined 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,
          A is overwritten by details of its QR or LQ
          factorization as returned by DGEQR or DGELQ.
[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.
          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.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= MAX(1,M,N).
[out]WORK
          (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
          or optimal, if query was assumed) LWORK.
          See LWORK for details.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          If LWORK = -1 or -2, then a workspace query is assumed.
          If LWORK = -1, the routine calculates optimal size of WORK for the
          optimal performance and returns this value in WORK(1).
          If LWORK = -2, the routine calculates minimal size of WORK and 
          returns this value in WORK(1).
[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
June 2017

Definition at line 164 of file dgetsls.f.

164 *
165 * -- LAPACK driver routine (version 3.7.1) --
166 * -- LAPACK is a software package provided by Univ. of Tennessee, --
167 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168 * June 2017
169 *
170 * .. Scalar Arguments ..
171  CHARACTER TRANS
172  INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
173 * ..
174 * .. Array Arguments ..
175  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
176 *
177 * ..
178 *
179 * =====================================================================
180 *
181 * .. Parameters ..
182  DOUBLE PRECISION ZERO, ONE
183  parameter( zero = 0.0d0, one = 1.0d0 )
184 * ..
185 * .. Local Scalars ..
186  LOGICAL LQUERY, TRAN
187  INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
188  $ SCLLEN, MNK, TSZO, TSZM, LWO, LWM, LW1, LW2,
189  $ WSIZEO, WSIZEM, INFO2
190  DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM, TQ( 5 ), WORKQ( 1 )
191 * ..
192 * .. External Functions ..
193  LOGICAL LSAME
194  INTEGER ILAENV
195  DOUBLE PRECISION DLAMCH, DLANGE
196  EXTERNAL lsame, ilaenv, dlabad, dlamch, dlange
197 * ..
198 * .. External Subroutines ..
199  EXTERNAL dgeqr, dgemqr, dlascl, dlaset,
200  $ dtrtrs, xerbla, dgelq, dgemlq
201 * ..
202 * .. Intrinsic Functions ..
203  INTRINSIC dble, max, min, int
204 * ..
205 * .. Executable Statements ..
206 *
207 * Test the input arguments.
208 *
209  info = 0
210  minmn = min( m, n )
211  maxmn = max( m, n )
212  mnk = max( minmn, nrhs )
213  tran = lsame( trans, 'T' )
214 *
215  lquery = ( lwork.EQ.-1 .OR. lwork.EQ.-2 )
216  IF( .NOT.( lsame( trans, 'N' ) .OR.
217  $ lsame( trans, 'T' ) ) ) THEN
218  info = -1
219  ELSE IF( m.LT.0 ) THEN
220  info = -2
221  ELSE IF( n.LT.0 ) THEN
222  info = -3
223  ELSE IF( nrhs.LT.0 ) THEN
224  info = -4
225  ELSE IF( lda.LT.max( 1, m ) ) THEN
226  info = -6
227  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
228  info = -8
229  END IF
230 *
231  IF( info.EQ.0 ) THEN
232 *
233 * Determine the block size and minimum LWORK
234 *
235  IF( m.GE.n ) THEN
236  CALL dgeqr( m, n, a, lda, tq, -1, workq, -1, info2 )
237  tszo = int( tq( 1 ) )
238  lwo = int( workq( 1 ) )
239  CALL dgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
240  $ tszo, b, ldb, workq, -1, info2 )
241  lwo = max( lwo, int( workq( 1 ) ) )
242  CALL dgeqr( m, n, a, lda, tq, -2, workq, -2, info2 )
243  tszm = int( tq( 1 ) )
244  lwm = int( workq( 1 ) )
245  CALL dgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
246  $ tszm, b, ldb, workq, -1, info2 )
247  lwm = max( lwm, int( workq( 1 ) ) )
248  wsizeo = tszo + lwo
249  wsizem = tszm + lwm
250  ELSE
251  CALL dgelq( m, n, a, lda, tq, -1, workq, -1, info2 )
252  tszo = int( tq( 1 ) )
253  lwo = int( workq( 1 ) )
254  CALL dgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
255  $ tszo, b, ldb, workq, -1, info2 )
256  lwo = max( lwo, int( workq( 1 ) ) )
257  CALL dgelq( m, n, a, lda, tq, -2, workq, -2, info2 )
258  tszm = int( tq( 1 ) )
259  lwm = int( workq( 1 ) )
260  CALL dgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
261  $ tszo, b, ldb, workq, -1, info2 )
262  lwm = max( lwm, int( workq( 1 ) ) )
263  wsizeo = tszo + lwo
264  wsizem = tszm + lwm
265  END IF
266 *
267  IF( ( lwork.LT.wsizem ).AND.( .NOT.lquery ) ) THEN
268  info = -10
269  END IF
270 *
271  END IF
272 *
273  IF( info.NE.0 ) THEN
274  CALL xerbla( 'DGETSLS', -info )
275  work( 1 ) = dble( wsizeo )
276  RETURN
277  END IF
278  IF( lquery ) THEN
279  IF( lwork.EQ.-1 ) work( 1 ) = real( wsizeo )
280  IF( lwork.EQ.-2 ) work( 1 ) = real( wsizem )
281  RETURN
282  END IF
283  IF( lwork.LT.wsizeo ) THEN
284  lw1 = tszm
285  lw2 = lwm
286  ELSE
287  lw1 = tszo
288  lw2 = lwo
289  END IF
290 *
291 * Quick return if possible
292 *
293  IF( min( m, n, nrhs ).EQ.0 ) THEN
294  CALL dlaset( 'FULL', max( m, n ), nrhs, zero, zero,
295  $ b, ldb )
296  RETURN
297  END IF
298 *
299 * Get machine parameters
300 *
301  smlnum = dlamch( 'S' ) / dlamch( 'P' )
302  bignum = one / smlnum
303  CALL dlabad( smlnum, bignum )
304 *
305 * Scale A, B if max element outside range [SMLNUM,BIGNUM]
306 *
307  anrm = dlange( 'M', m, n, a, lda, work )
308  iascl = 0
309  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
310 *
311 * Scale matrix norm up to SMLNUM
312 *
313  CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
314  iascl = 1
315  ELSE IF( anrm.GT.bignum ) THEN
316 *
317 * Scale matrix norm down to BIGNUM
318 *
319  CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
320  iascl = 2
321  ELSE IF( anrm.EQ.zero ) THEN
322 *
323 * Matrix all zero. Return zero solution.
324 *
325  CALL dlaset( 'F', maxmn, nrhs, zero, zero, b, ldb )
326  GO TO 50
327  END IF
328 *
329  brow = m
330  IF ( tran ) THEN
331  brow = n
332  END IF
333  bnrm = dlange( 'M', brow, nrhs, b, ldb, work )
334  ibscl = 0
335  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
336 *
337 * Scale matrix norm up to SMLNUM
338 *
339  CALL dlascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
340  $ info )
341  ibscl = 1
342  ELSE IF( bnrm.GT.bignum ) THEN
343 *
344 * Scale matrix norm down to BIGNUM
345 *
346  CALL dlascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
347  $ info )
348  ibscl = 2
349  END IF
350 *
351  IF ( m.GE.n ) THEN
352 *
353 * compute QR factorization of A
354 *
355  CALL dgeqr( m, n, a, lda, work( lw2+1 ), lw1,
356  $ work( 1 ), lw2, info )
357  IF ( .NOT.tran ) 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 dgemqr( 'L' , 'T', m, nrhs, n, a, lda,
364  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
365  $ info )
366 *
367 * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
368 *
369  CALL dtrtrs( 'U', 'N', 'N', n, nrhs,
370  $ a, lda, b, ldb, info )
371  IF( info.GT.0 ) THEN
372  RETURN
373  END IF
374  scllen = n
375  ELSE
376 *
377 * Overdetermined system of equations A**T * X = B
378 *
379 * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
380 *
381  CALL dtrtrs( 'U', 'T', 'N', n, nrhs,
382  $ a, lda, b, ldb, info )
383 *
384  IF( info.GT.0 ) THEN
385  RETURN
386  END IF
387 *
388 * B(N+1:M,1:NRHS) = ZERO
389 *
390  DO 20 j = 1, nrhs
391  DO 10 i = n + 1, m
392  b( i, j ) = zero
393  10 CONTINUE
394  20 CONTINUE
395 *
396 * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
397 *
398  CALL dgemqr( 'L', 'N', m, nrhs, n, a, lda,
399  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
400  $ info )
401 *
402  scllen = m
403 *
404  END IF
405 *
406  ELSE
407 *
408 * Compute LQ factorization of A
409 *
410  CALL dgelq( m, n, a, lda, work( lw2+1 ), lw1,
411  $ work( 1 ), lw2, info )
412 *
413 * workspace at least M, optimally M*NB.
414 *
415  IF( .NOT.tran ) THEN
416 *
417 * underdetermined system of equations A * X = B
418 *
419 * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
420 *
421  CALL dtrtrs( 'L', 'N', 'N', m, nrhs,
422  $ a, lda, b, ldb, info )
423 *
424  IF( info.GT.0 ) THEN
425  RETURN
426  END IF
427 *
428 * B(M+1:N,1:NRHS) = 0
429 *
430  DO 40 j = 1, nrhs
431  DO 30 i = m + 1, n
432  b( i, j ) = zero
433  30 CONTINUE
434  40 CONTINUE
435 *
436 * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
437 *
438  CALL dgemlq( 'L', 'T', n, nrhs, m, a, lda,
439  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
440  $ info )
441 *
442 * workspace at least NRHS, optimally NRHS*NB
443 *
444  scllen = n
445 *
446  ELSE
447 *
448 * overdetermined system min || A**T * X - B ||
449 *
450 * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
451 *
452  CALL dgemlq( 'L', 'N', n, nrhs, m, a, lda,
453  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
454  $ info )
455 *
456 * workspace at least NRHS, optimally NRHS*NB
457 *
458 * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
459 *
460  CALL dtrtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
461  $ a, lda, b, ldb, info )
462 *
463  IF( info.GT.0 ) THEN
464  RETURN
465  END IF
466 *
467  scllen = m
468 *
469  END IF
470 *
471  END IF
472 *
473 * Undo scaling
474 *
475  IF( iascl.EQ.1 ) THEN
476  CALL dlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
477  $ info )
478  ELSE IF( iascl.EQ.2 ) THEN
479  CALL dlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
480  $ info )
481  END IF
482  IF( ibscl.EQ.1 ) THEN
483  CALL dlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
484  $ info )
485  ELSE IF( ibscl.EQ.2 ) THEN
486  CALL dlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
487  $ info )
488  END IF
489 *
490  50 CONTINUE
491  work( 1 ) = dble( tszo + lwo )
492  RETURN
493 *
494 * End of DGETSLS
495 *
Here is the call graph for this function:
Here is the caller graph for this function:
dlascl
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
dtrtrs
subroutine dtrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
DTRTRS
Definition: dtrtrs.f:142
dgemqr
subroutine dgemqr(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
DGEMQR
Definition: dgemqr.f:172
dlange
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
dgelq
subroutine dgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
DGELQ
Definition: dgelq.f:172
dgeqr
subroutine dgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
DGEQR
Definition: dgeqr.f:174
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
dlaset
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
dgemlq
subroutine dgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
DGEMLQ
Definition: dgemlq.f:171
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:70