LAPACK  3.10.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.

Definition at line 180 of file zgels.f.

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