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

## ◆ 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.

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.

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 zero, so that A does not have full rank; the least squares solution could not be computed.```
Contributors:
```  November 2022,  Igor Kozachenko,
Computer Science Division,
University of California, Berkeley```

Definition at line 192 of file zgelst.f.

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