LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ zgetsls()

subroutine zgetsls ( character  TRANS,
integer  M,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGETSLS

Purpose:
 ZGETSLS solves overdetermined or underdetermined complex 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 = 'C' and m >= n:  find the minimum norm solution of
    an undetermined system A**T * X = B.

 4. If TRANS = 'C' 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;
          = 'C': the linear system involves A**H.
[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 COMPLEX*16 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 ZGEQR or ZGELQ.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[in,out]B
          B is COMPLEX*16 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 = 'C'.
          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 = 'C' and m >= n, rows 1 to M of B contain the
          minimum norm solution vectors;
          if TRANS = 'C' 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) COMPLEX*16 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 zgetsls.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  COMPLEX*16 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  COMPLEX*16 CZERO
185  parameter( czero = ( 0.0d+0, 0.0d+0 ) )
186 * ..
187 * .. Local Scalars ..
188  LOGICAL LQUERY, TRAN
189  INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
190  $ SCLLEN, MNK, TSZO, TSZM, LWO, LWM, LW1, LW2,
191  $ WSIZEO, WSIZEM, INFO2
192  DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM, DUM( 1 )
193  COMPLEX*16 TQ( 5 ), WORKQ( 1 )
194 * ..
195 * .. External Functions ..
196  LOGICAL LSAME
197  INTEGER ILAENV
198  DOUBLE PRECISION DLAMCH, ZLANGE
199  EXTERNAL lsame, ilaenv, dlabad, dlamch, zlange
200 * ..
201 * .. External Subroutines ..
202  EXTERNAL zgeqr, zgemqr, zlascl, zlaset,
203  $ ztrtrs, xerbla, zgelq, zgemlq
204 * ..
205 * .. Intrinsic Functions ..
206  INTRINSIC dble, max, min, int
207 * ..
208 * .. Executable Statements ..
209 *
210 * Test the input arguments.
211 *
212  info = 0
213  minmn = min( m, n )
214  maxmn = max( m, n )
215  mnk = max( minmn, nrhs )
216  tran = lsame( trans, 'C' )
217 *
218  lquery = ( lwork.EQ.-1 .OR. lwork.EQ.-2 )
219  IF( .NOT.( lsame( trans, 'N' ) .OR.
220  $ lsame( trans, 'C' ) ) ) THEN
221  info = -1
222  ELSE IF( m.LT.0 ) THEN
223  info = -2
224  ELSE IF( n.LT.0 ) THEN
225  info = -3
226  ELSE IF( nrhs.LT.0 ) THEN
227  info = -4
228  ELSE IF( lda.LT.max( 1, m ) ) THEN
229  info = -6
230  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
231  info = -8
232  END IF
233 *
234  IF( info.EQ.0 ) THEN
235 *
236 * Determine the block size and minimum LWORK
237 *
238  IF( m.GE.n ) THEN
239  CALL zgeqr( m, n, a, lda, tq, -1, workq, -1, info2 )
240  tszo = int( tq( 1 ) )
241  lwo = int( workq( 1 ) )
242  CALL zgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
243  $ tszo, b, ldb, workq, -1, info2 )
244  lwo = max( lwo, int( workq( 1 ) ) )
245  CALL zgeqr( m, n, a, lda, tq, -2, workq, -2, info2 )
246  tszm = int( tq( 1 ) )
247  lwm = int( workq( 1 ) )
248  CALL zgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
249  $ tszm, b, ldb, workq, -1, info2 )
250  lwm = max( lwm, int( workq( 1 ) ) )
251  wsizeo = tszo + lwo
252  wsizem = tszm + lwm
253  ELSE
254  CALL zgelq( m, n, a, lda, tq, -1, workq, -1, info2 )
255  tszo = int( tq( 1 ) )
256  lwo = int( workq( 1 ) )
257  CALL zgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
258  $ tszo, b, ldb, workq, -1, info2 )
259  lwo = max( lwo, int( workq( 1 ) ) )
260  CALL zgelq( m, n, a, lda, tq, -2, workq, -2, info2 )
261  tszm = int( tq( 1 ) )
262  lwm = int( workq( 1 ) )
263  CALL zgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
264  $ tszo, b, ldb, workq, -1, info2 )
265  lwm = max( lwm, int( workq( 1 ) ) )
266  wsizeo = tszo + lwo
267  wsizem = tszm + lwm
268  END IF
269 *
270  IF( ( lwork.LT.wsizem ).AND.( .NOT.lquery ) ) THEN
271  info = -10
272  END IF
273 *
274  END IF
275 *
276  IF( info.NE.0 ) THEN
277  CALL xerbla( 'ZGETSLS', -info )
278  work( 1 ) = dble( wsizeo )
279  RETURN
280  END IF
281  IF( lquery ) THEN
282  IF( lwork.EQ.-1 ) work( 1 ) = real( wsizeo )
283  IF( lwork.EQ.-2 ) work( 1 ) = real( wsizem )
284  RETURN
285  END IF
286  IF( lwork.LT.wsizeo ) THEN
287  lw1 = tszm
288  lw2 = lwm
289  ELSE
290  lw1 = tszo
291  lw2 = lwo
292  END IF
293 *
294 * Quick return if possible
295 *
296  IF( min( m, n, nrhs ).EQ.0 ) THEN
297  CALL zlaset( 'FULL', max( m, n ), nrhs, czero, czero,
298  $ b, ldb )
299  RETURN
300  END IF
301 *
302 * Get machine parameters
303 *
304  smlnum = dlamch( 'S' ) / dlamch( 'P' )
305  bignum = one / smlnum
306  CALL dlabad( smlnum, bignum )
307 *
308 * Scale A, B if max element outside range [SMLNUM,BIGNUM]
309 *
310  anrm = zlange( 'M', m, n, a, lda, dum )
311  iascl = 0
312  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
313 *
314 * Scale matrix norm up to SMLNUM
315 *
316  CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
317  iascl = 1
318  ELSE IF( anrm.GT.bignum ) THEN
319 *
320 * Scale matrix norm down to BIGNUM
321 *
322  CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
323  iascl = 2
324  ELSE IF( anrm.EQ.zero ) THEN
325 *
326 * Matrix all zero. Return zero solution.
327 *
328  CALL zlaset( 'F', maxmn, nrhs, czero, czero, b, ldb )
329  GO TO 50
330  END IF
331 *
332  brow = m
333  IF ( tran ) THEN
334  brow = n
335  END IF
336  bnrm = zlange( 'M', brow, nrhs, b, ldb, dum )
337  ibscl = 0
338  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
339 *
340 * Scale matrix norm up to SMLNUM
341 *
342  CALL zlascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
343  $ info )
344  ibscl = 1
345  ELSE IF( bnrm.GT.bignum ) THEN
346 *
347 * Scale matrix norm down to BIGNUM
348 *
349  CALL zlascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
350  $ info )
351  ibscl = 2
352  END IF
353 *
354  IF ( m.GE.n ) THEN
355 *
356 * compute QR factorization of A
357 *
358  CALL zgeqr( m, n, a, lda, work( lw2+1 ), lw1,
359  $ work( 1 ), lw2, info )
360  IF ( .NOT.tran ) THEN
361 *
362 * Least-Squares Problem min || A * X - B ||
363 *
364 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
365 *
366  CALL zgemqr( 'L' , 'C', m, nrhs, n, a, lda,
367  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
368  $ info )
369 *
370 * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
371 *
372  CALL ztrtrs( 'U', 'N', 'N', n, nrhs,
373  $ a, lda, b, ldb, info )
374  IF( info.GT.0 ) THEN
375  RETURN
376  END IF
377  scllen = n
378  ELSE
379 *
380 * Overdetermined system of equations A**T * X = B
381 *
382 * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
383 *
384  CALL ztrtrs( 'U', 'C', 'N', n, nrhs,
385  $ a, lda, b, ldb, info )
386 *
387  IF( info.GT.0 ) THEN
388  RETURN
389  END IF
390 *
391 * B(N+1:M,1:NRHS) = CZERO
392 *
393  DO 20 j = 1, nrhs
394  DO 10 i = n + 1, m
395  b( i, j ) = czero
396  10 CONTINUE
397  20 CONTINUE
398 *
399 * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
400 *
401  CALL zgemqr( 'L', 'N', m, nrhs, n, a, lda,
402  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
403  $ info )
404 *
405  scllen = m
406 *
407  END IF
408 *
409  ELSE
410 *
411 * Compute LQ factorization of A
412 *
413  CALL zgelq( m, n, a, lda, work( lw2+1 ), lw1,
414  $ work( 1 ), lw2, info )
415 *
416 * workspace at least M, optimally M*NB.
417 *
418  IF( .NOT.tran ) 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 ztrtrs( 'L', 'N', 'N', 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 ) = czero
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 zgemlq( 'L', 'C', n, nrhs, m, a, lda,
442  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
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 zgemlq( 'L', 'N', n, nrhs, m, a, lda,
456  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
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 ztrtrs( 'L', 'C', 'N', 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 zlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
480  $ info )
481  ELSE IF( iascl.EQ.2 ) THEN
482  CALL zlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
483  $ info )
484  END IF
485  IF( ibscl.EQ.1 ) THEN
486  CALL zlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
487  $ info )
488  ELSE IF( ibscl.EQ.2 ) THEN
489  CALL zlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
490  $ info )
491  END IF
492 *
493  50 CONTINUE
494  work( 1 ) = dble( tszo + lwo )
495  RETURN
496 *
497 * End of ZGETSLS
498 *
Here is the call graph for this function:
Here is the caller graph for this function:
zlange
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
zgelq
subroutine zgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
ZGELQ
Definition: zgelq.f:172
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:108
zgeqr
subroutine zgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
ZGEQR
Definition: zgeqr.f:174
zlascl
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:145
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
ztrtrs
subroutine ztrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
ZTRTRS
Definition: ztrtrs.f:142
zgemqr
subroutine zgemqr(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
ZGEMQR
Definition: zgemqr.f:172
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
zgemlq
subroutine zgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
ZGEMLQ
Definition: zgemlq.f:169
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:70