LAPACK 3.12.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 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*
300* Scale A, B if max element outside range [SMLNUM,BIGNUM]
301*
302 anrm = zlange( 'M', m, n, a, lda, rwork )
303 iascl = 0
304 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
305*
306* Scale matrix norm up to SMLNUM
307*
308 CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
309 iascl = 1
310 ELSE IF( anrm.GT.bignum ) THEN
311*
312* Scale matrix norm down to BIGNUM
313*
314 CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
315 iascl = 2
316 ELSE IF( anrm.EQ.zero ) THEN
317*
318* Matrix all zero. Return zero solution.
319*
320 CALL zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
321 GO TO 50
322 END IF
323*
324 brow = m
325 IF( tpsd )
326 $ brow = n
327 bnrm = zlange( 'M', brow, nrhs, b, ldb, rwork )
328 ibscl = 0
329 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
330*
331* Scale matrix norm up to SMLNUM
332*
333 CALL zlascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
334 $ info )
335 ibscl = 1
336 ELSE IF( bnrm.GT.bignum ) THEN
337*
338* Scale matrix norm down to BIGNUM
339*
340 CALL zlascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
341 $ info )
342 ibscl = 2
343 END IF
344*
345 IF( m.GE.n ) THEN
346*
347* compute QR factorization of A
348*
349 CALL zgeqrf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
350 $ info )
351*
352* workspace at least N, optimally N*NB
353*
354 IF( .NOT.tpsd ) THEN
355*
356* Least-Squares Problem min || A * X - B ||
357*
358* B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
359*
360 CALL zunmqr( 'Left', 'Conjugate transpose', m, nrhs, n, a,
361 $ lda, work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
362 $ info )
363*
364* workspace at least NRHS, optimally NRHS*NB
365*
366* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
367*
368 CALL ztrtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
369 $ a, lda, b, ldb, info )
370*
371 IF( info.GT.0 ) THEN
372 RETURN
373 END IF
374*
375 scllen = n
376*
377 ELSE
378*
379* Underdetermined system of equations A**T * X = B
380*
381* B(1:N,1:NRHS) := inv(R**H) * B(1:N,1:NRHS)
382*
383 CALL ztrtrs( 'Upper', 'Conjugate transpose','Non-unit',
384 $ n, nrhs, a, lda, b, ldb, info )
385*
386 IF( info.GT.0 ) THEN
387 RETURN
388 END IF
389*
390* B(N+1:M,1:NRHS) = ZERO
391*
392 DO 20 j = 1, nrhs
393 DO 10 i = n + 1, m
394 b( i, j ) = czero
395 10 CONTINUE
396 20 CONTINUE
397*
398* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
399*
400 CALL zunmqr( 'Left', 'No transpose', m, nrhs, n, a, lda,
401 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
402 $ info )
403*
404* workspace at least NRHS, optimally NRHS*NB
405*
406 scllen = m
407*
408 END IF
409*
410 ELSE
411*
412* Compute LQ factorization of A
413*
414 CALL zgelqf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
415 $ info )
416*
417* workspace at least M, optimally M*NB.
418*
419 IF( .NOT.tpsd ) THEN
420*
421* underdetermined system of equations A * X = B
422*
423* B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
424*
425 CALL ztrtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
426 $ a, lda, b, ldb, info )
427*
428 IF( info.GT.0 ) THEN
429 RETURN
430 END IF
431*
432* B(M+1:N,1:NRHS) = 0
433*
434 DO 40 j = 1, nrhs
435 DO 30 i = m + 1, n
436 b( i, j ) = czero
437 30 CONTINUE
438 40 CONTINUE
439*
440* B(1:N,1:NRHS) := Q(1:N,:)**H * B(1:M,1:NRHS)
441*
442 CALL zunmlq( 'Left', 'Conjugate transpose', n, nrhs, m, a,
443 $ lda, work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
444 $ info )
445*
446* workspace at least NRHS, optimally NRHS*NB
447*
448 scllen = n
449*
450 ELSE
451*
452* overdetermined system min || A**H * X - B ||
453*
454* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
455*
456 CALL zunmlq( 'Left', 'No transpose', n, nrhs, m, a, lda,
457 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
458 $ info )
459*
460* workspace at least NRHS, optimally NRHS*NB
461*
462* B(1:M,1:NRHS) := inv(L**H) * B(1:M,1:NRHS)
463*
464 CALL ztrtrs( 'Lower', 'Conjugate transpose', 'Non-unit',
465 $ m, nrhs, a, lda, b, ldb, info )
466*
467 IF( info.GT.0 ) THEN
468 RETURN
469 END IF
470*
471 scllen = m
472*
473 END IF
474*
475 END IF
476*
477* Undo scaling
478*
479 IF( iascl.EQ.1 ) THEN
480 CALL zlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
481 $ info )
482 ELSE IF( iascl.EQ.2 ) THEN
483 CALL zlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
484 $ info )
485 END IF
486 IF( ibscl.EQ.1 ) THEN
487 CALL zlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
488 $ info )
489 ELSE IF( ibscl.EQ.2 ) THEN
490 CALL zlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
491 $ info )
492 END IF
493*
494 50 CONTINUE
495 work( 1 ) = dble( wsize )
496*
497 RETURN
498*
499* End of ZGELS
500*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zgelqf(m, n, a, lda, tau, work, lwork, info)
ZGELQF
Definition zgelqf.f:143
subroutine zgeqrf(m, n, a, lda, tau, work, lwork, info)
ZGEQRF
Definition zgeqrf.f:146
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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 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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine ztrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
ZTRTRS
Definition ztrtrs.f:140
subroutine zunmlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMLQ
Definition zunmlq.f:167
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR
Definition zunmqr.f:167
Here is the call graph for this function:
Here is the caller graph for this function: