LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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,
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:152
Here is the call graph for this function:
Here is the caller graph for this function: