LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages

◆ zgelst()

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

ZGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization with compact WY representation of Q.

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

Purpose:
!>
!> ZGELST solves overdetermined or underdetermined real linear systems
!> involving an M-by-N matrix A, or its conjugate-transpose, using a QR
!> or LQ factorization of A with compact WY representation of Q.
!>
!> It is assumed that A has full rank, and only a rudimentary protection
!> against rank-deficient matrices is provided. This subroutine only detects
!> exact rank-deficiency, where a diagonal element of the triangular factor
!> of A is exactly zero.
!>
!> It is conceivable for one (or more) of the diagonal elements of the triangular
!> factor of A to be subnormally tiny numbers without this subroutine signalling
!> an error. The solutions computed for such almost-rank-deficient matrices may
!> be less accurate due to a loss of numerical precision.
!>
!> 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**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**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.
!>          On exit,
!>            if M >= N, A is overwritten by details of its QR
!>                       factorization as returned by ZGEQRT;
!>            if M <  N, A is overwritten by details of its LQ
!>                       factorization as returned by ZGELQT.
!>
[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
!>          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 exactly 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.
Contributors:
!>
!>  November 2022,  Igor Kozachenko,
!>                  Computer Science Division,
!>                  University of California, Berkeley
!>

Definition at line 199 of file zgelst.f.

202*
203* -- LAPACK driver routine --
204* -- LAPACK is a software package provided by Univ. of Tennessee, --
205* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
206*
207* .. Scalar Arguments ..
208 CHARACTER TRANS
209 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
210* ..
211* .. Array Arguments ..
212 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
213* ..
214*
215* =====================================================================
216*
217* .. Parameters ..
218 DOUBLE PRECISION ZERO, ONE
219 parameter( zero = 0.0d+0, one = 1.0d+0 )
220 COMPLEX*16 CZERO
221 parameter( czero = ( 0.0d+0, 0.0d+0 ) )
222* ..
223* .. Local Scalars ..
224 LOGICAL LQUERY, TPSD
225 INTEGER BROW, I, IASCL, IBSCL, J, LWOPT, MN, MNNRHS,
226 $ NB, NBMIN, SCLLEN
227 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
228* ..
229* .. Local Arrays ..
230 DOUBLE PRECISION RWORK( 1 )
231* ..
232* .. External Functions ..
233 LOGICAL LSAME
234 INTEGER ILAENV
235 DOUBLE PRECISION DLAMCH, ZLANGE
236 EXTERNAL lsame, ilaenv, dlamch, zlange
237* ..
238* .. External Subroutines ..
239 EXTERNAL zgelqt, zgeqrt, zgemlqt, zgemqrt,
240 $ zlascl,
241 $ zlaset, ztrtrs, xerbla
242* ..
243* .. Intrinsic Functions ..
244 INTRINSIC dble, max, min
245* ..
246* .. Executable Statements ..
247*
248* Test the input arguments.
249*
250 info = 0
251 mn = min( m, n )
252 lquery = ( lwork.EQ.-1 )
253 IF( .NOT.( lsame( trans, 'N' ) .OR.
254 $ lsame( trans, 'C' ) ) ) THEN
255 info = -1
256 ELSE IF( m.LT.0 ) THEN
257 info = -2
258 ELSE IF( n.LT.0 ) THEN
259 info = -3
260 ELSE IF( nrhs.LT.0 ) THEN
261 info = -4
262 ELSE IF( lda.LT.max( 1, m ) ) THEN
263 info = -6
264 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
265 info = -8
266 ELSE IF( lwork.LT.max( 1, mn+max( mn, nrhs ) ) .AND. .NOT.lquery )
267 $ THEN
268 info = -10
269 END IF
270*
271* Figure out optimal block size and optimal workspace size
272*
273 IF( info.EQ.0 .OR. info.EQ.-10 ) THEN
274*
275 tpsd = .true.
276 IF( lsame( trans, 'N' ) )
277 $ tpsd = .false.
278*
279 nb = ilaenv( 1, 'ZGELST', ' ', m, n, -1, -1 )
280*
281 mnnrhs = max( mn, nrhs )
282 lwopt = max( 1, (mn+mnnrhs)*nb )
283 work( 1 ) = dble( lwopt )
284*
285 END IF
286*
287 IF( info.NE.0 ) THEN
288 CALL xerbla( 'ZGELST ', -info )
289 RETURN
290 ELSE IF( lquery ) THEN
291 RETURN
292 END IF
293*
294* Quick return if possible
295*
296 IF( min( m, n, nrhs ).EQ.0 ) THEN
297 CALL zlaset( 'Full', max( m, n ), nrhs, czero, czero, b,
298 $ ldb )
299 work( 1 ) = dble( lwopt )
300 RETURN
301 END IF
302*
303* *GEQRT and *GELQT routines cannot accept NB larger than min(M,N)
304*
305 IF( nb.GT.mn ) nb = mn
306*
307* Determine the block size from the supplied LWORK
308* ( at this stage we know that LWORK >= (minimum required workspace,
309* but it may be less than optimal)
310*
311 nb = min( nb, lwork/( mn + mnnrhs ) )
312*
313* The minimum value of NB, when blocked code is used
314*
315 nbmin = max( 2, ilaenv( 2, 'ZGELST', ' ', m, n, -1, -1 ) )
316*
317 IF( nb.LT.nbmin ) THEN
318 nb = 1
319 END IF
320*
321* Get machine parameters
322*
323 smlnum = dlamch( 'S' ) / dlamch( 'P' )
324 bignum = one / smlnum
325*
326* Scale A, B if max element outside range [SMLNUM,BIGNUM]
327*
328 anrm = zlange( 'M', m, n, a, lda, rwork )
329 iascl = 0
330 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
331*
332* Scale matrix norm up to SMLNUM
333*
334 CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
335 iascl = 1
336 ELSE IF( anrm.GT.bignum ) THEN
337*
338* Scale matrix norm down to BIGNUM
339*
340 CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
341 iascl = 2
342 ELSE IF( anrm.EQ.zero ) THEN
343*
344* Matrix all zero. Return zero solution.
345*
346 CALL zlaset( 'Full', max( m, n ), nrhs, czero, czero, b,
347 $ ldb )
348 work( 1 ) = dble( lwopt )
349 RETURN
350 END IF
351*
352 brow = m
353 IF( tpsd )
354 $ brow = n
355 bnrm = zlange( 'M', brow, nrhs, b, ldb, rwork )
356 ibscl = 0
357 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
358*
359* Scale matrix norm up to SMLNUM
360*
361 CALL zlascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
362 $ info )
363 ibscl = 1
364 ELSE IF( bnrm.GT.bignum ) THEN
365*
366* Scale matrix norm down to BIGNUM
367*
368 CALL zlascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
369 $ info )
370 ibscl = 2
371 END IF
372*
373 IF( m.GE.n ) THEN
374*
375* M > N:
376* Compute the blocked QR factorization of A,
377* using the compact WY representation of Q,
378* workspace at least N, optimally N*NB.
379*
380 CALL zgeqrt( m, n, nb, a, lda, work( 1 ), nb,
381 $ work( mn*nb+1 ), info )
382*
383 IF( .NOT.tpsd ) THEN
384*
385* M > N, A is not transposed:
386* Overdetermined system of equations,
387* least-squares problem, min || A * X - B ||.
388*
389* Compute B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS),
390* using the compact WY representation of Q,
391* workspace at least NRHS, optimally NRHS*NB.
392*
393 CALL zgemqrt( 'Left', 'Conjugate transpose', m, nrhs, n,
394 $ nb,
395 $ a, lda, work( 1 ), nb, b, ldb,
396 $ work( mn*nb+1 ), info )
397*
398* Compute B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
399*
400 CALL ztrtrs( 'Upper', 'No transpose', 'Non-unit', n,
401 $ nrhs,
402 $ a, lda, b, ldb, info )
403*
404 IF( info.GT.0 ) THEN
405 RETURN
406 END IF
407*
408 scllen = n
409*
410 ELSE
411*
412* M > N, A is transposed:
413* Underdetermined system of equations,
414* minimum norm solution of A**T * X = B.
415*
416* Compute B := inv(R**T) * B in two row blocks of B.
417*
418* Block 1: B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
419*
420 CALL ztrtrs( 'Upper', 'Conjugate transpose', 'Non-unit',
421 $ n, nrhs, a, lda, b, ldb, info )
422*
423 IF( info.GT.0 ) THEN
424 RETURN
425 END IF
426*
427* Block 2: Zero out all rows below the N-th row in B:
428* B(N+1:M,1:NRHS) = ZERO
429*
430 DO j = 1, nrhs
431 DO i = n + 1, m
432 b( i, j ) = zero
433 END DO
434 END DO
435*
436* Compute B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS),
437* using the compact WY representation of Q,
438* workspace at least NRHS, optimally NRHS*NB.
439*
440 CALL zgemqrt( 'Left', 'No transpose', m, nrhs, n, nb,
441 $ a, lda, work( 1 ), nb, b, ldb,
442 $ work( mn*nb+1 ), info )
443*
444 scllen = m
445*
446 END IF
447*
448 ELSE
449*
450* M < N:
451* Compute the blocked LQ factorization of A,
452* using the compact WY representation of Q,
453* workspace at least M, optimally M*NB.
454*
455 CALL zgelqt( m, n, nb, a, lda, work( 1 ), nb,
456 $ work( mn*nb+1 ), info )
457*
458 IF( .NOT.tpsd ) THEN
459*
460* M < N, A is not transposed:
461* Underdetermined system of equations,
462* minimum norm solution of A * X = B.
463*
464* Compute B := inv(L) * B in two row blocks of B.
465*
466* Block 1: B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
467*
468 CALL ztrtrs( 'Lower', 'No transpose', 'Non-unit', m,
469 $ nrhs,
470 $ a, lda, b, ldb, info )
471*
472 IF( info.GT.0 ) THEN
473 RETURN
474 END IF
475*
476* Block 2: Zero out all rows below the M-th row in B:
477* B(M+1:N,1:NRHS) = ZERO
478*
479 DO j = 1, nrhs
480 DO i = m + 1, n
481 b( i, j ) = zero
482 END DO
483 END DO
484*
485* Compute B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS),
486* using the compact WY representation of Q,
487* workspace at least NRHS, optimally NRHS*NB.
488*
489 CALL zgemlqt( 'Left', 'Conjugate transpose', n, nrhs, m,
490 $ nb,
491 $ a, lda, work( 1 ), nb, b, ldb,
492 $ work( mn*nb+1 ), info )
493*
494 scllen = n
495*
496 ELSE
497*
498* M < N, A is transposed:
499* Overdetermined system of equations,
500* least-squares problem, min || A**T * X - B ||.
501*
502* Compute B(1:N,1:NRHS) := Q * B(1:N,1:NRHS),
503* using the compact WY representation of Q,
504* workspace at least NRHS, optimally NRHS*NB.
505*
506 CALL zgemlqt( 'Left', 'No transpose', n, nrhs, m, nb,
507 $ a, lda, work( 1 ), nb, b, ldb,
508 $ work( mn*nb+1), info )
509*
510* Compute B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
511*
512 CALL ztrtrs( 'Lower', 'Conjugate transpose', 'Non-unit',
513 $ m, nrhs, a, lda, b, ldb, info )
514*
515 IF( info.GT.0 ) THEN
516 RETURN
517 END IF
518*
519 scllen = m
520*
521 END IF
522*
523 END IF
524*
525* Undo scaling
526*
527 IF( iascl.EQ.1 ) THEN
528 CALL zlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
529 $ info )
530 ELSE IF( iascl.EQ.2 ) THEN
531 CALL zlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
532 $ info )
533 END IF
534 IF( ibscl.EQ.1 ) THEN
535 CALL zlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
536 $ info )
537 ELSE IF( ibscl.EQ.2 ) THEN
538 CALL zlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
539 $ info )
540 END IF
541*
542 work( 1 ) = dble( lwopt )
543*
544 RETURN
545*
546* End of ZGELST
547*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgelqt
subroutine zgelqt(m, n, mb, a, lda, t, ldt, work, info)
ZGELQT
Definition zgelqt.f:137
zgemlqt
subroutine zgemlqt(side, trans, m, n, k, mb, v, ldv, t, ldt, c, ldc, work, info)
ZGEMLQT
Definition zgemlqt.f:166
zgemqrt
subroutine zgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
ZGEMQRT
Definition zgemqrt.f:166
zgeqrt
subroutine zgeqrt(m, n, nb, a, lda, t, ldt, work, info)
ZGEQRT
Definition zgeqrt.f:139
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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:113
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:142
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:104
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
ztrtrs
subroutine ztrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
ZTRTRS
Definition ztrtrs.f:144
Here is the call graph for this function:
Here is the caller graph for this function: