LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ zgels()

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

ZGELS solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
 ZGELS solves overdetermined or underdetermined complex linear systems
 involving an M-by-N matrix A, or its conjugate-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 = 'C' and m >= n:  find the minimum norm solution of
    an underdetermined system A**H * 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**H * 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.
            if M >= N, A is overwritten by details of its QR
                       factorization as returned by ZGEQRF;
            if M <  N, A is overwritten by details of its LQ
                       factorization as returned by ZGELQF.
[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; the residual sum of squares for the
          solution in each column is given by the sum of squares of the
          modulus 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 = '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; the residual sum of squares
          for the solution in each column is given by the sum of
          squares of the modulus 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 COMPLEX*16 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.
Date
December 2016

Definition at line 184 of file zgels.f.

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