LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
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 
)
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**C.
[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
December 2016

Definition at line 162 of file zgetsls.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: