LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zgelst.f
Go to the documentation of this file.
1*> \brief <b> ZGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization with compact WY representation of Q.</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelst.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelst.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelst.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZGELST( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
22* INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER TRANS
26* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
27* ..
28* .. Array Arguments ..
29* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> ZGELST solves overdetermined or underdetermined real linear systems
39*> involving an M-by-N matrix A, or its conjugate-transpose, using a QR
40*> or LQ factorization of A with compact WY representation of Q.
41*> It is assumed that A has full rank.
42*>
43*> The following options are provided:
44*>
45*> 1. If TRANS = 'N' and m >= n: find the least squares solution of
46*> an overdetermined system, i.e., solve the least squares problem
47*> minimize || B - A*X ||.
48*>
49*> 2. If TRANS = 'N' and m < n: find the minimum norm solution of
50*> an underdetermined system A * X = B.
51*>
52*> 3. If TRANS = 'C' and m >= n: find the minimum norm solution of
53*> an underdetermined system A**T * X = B.
54*>
55*> 4. If TRANS = 'C' and m < n: find the least squares solution of
56*> an overdetermined system, i.e., solve the least squares problem
57*> minimize || B - A**T * X ||.
58*>
59*> Several right hand side vectors b and solution vectors x can be
60*> handled in a single call; they are stored as the columns of the
61*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
62*> matrix X.
63*> \endverbatim
64*
65* Arguments:
66* ==========
67*
68*> \param[in] TRANS
69*> \verbatim
70*> TRANS is CHARACTER*1
71*> = 'N': the linear system involves A;
72*> = 'C': the linear system involves A**H.
73*> \endverbatim
74*>
75*> \param[in] M
76*> \verbatim
77*> M is INTEGER
78*> The number of rows of the matrix A. M >= 0.
79*> \endverbatim
80*>
81*> \param[in] N
82*> \verbatim
83*> N is INTEGER
84*> The number of columns of the matrix A. N >= 0.
85*> \endverbatim
86*>
87*> \param[in] NRHS
88*> \verbatim
89*> NRHS is INTEGER
90*> The number of right hand sides, i.e., the number of
91*> columns of the matrices B and X. NRHS >=0.
92*> \endverbatim
93*>
94*> \param[in,out] A
95*> \verbatim
96*> A is COMPLEX*16 array, dimension (LDA,N)
97*> On entry, the M-by-N matrix A.
98*> On exit,
99*> if M >= N, A is overwritten by details of its QR
100*> factorization as returned by ZGEQRT;
101*> if M < N, A is overwritten by details of its LQ
102*> factorization as returned by ZGELQT.
103*> \endverbatim
104*>
105*> \param[in] LDA
106*> \verbatim
107*> LDA is INTEGER
108*> The leading dimension of the array A. LDA >= max(1,M).
109*> \endverbatim
110*>
111*> \param[in,out] B
112*> \verbatim
113*> B is COMPLEX*16 array, dimension (LDB,NRHS)
114*> On entry, the matrix B of right hand side vectors, stored
115*> columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
116*> if TRANS = 'C'.
117*> On exit, if INFO = 0, B is overwritten by the solution
118*> vectors, stored columnwise:
119*> if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
120*> squares solution vectors; the residual sum of squares for the
121*> solution in each column is given by the sum of squares of
122*> modulus of elements N+1 to M in that column;
123*> if TRANS = 'N' and m < n, rows 1 to N of B contain the
124*> minimum norm solution vectors;
125*> if TRANS = 'C' and m >= n, rows 1 to M of B contain the
126*> minimum norm solution vectors;
127*> if TRANS = 'C' and m < n, rows 1 to M of B contain the
128*> least squares solution vectors; the residual sum of squares
129*> for the solution in each column is given by the sum of
130*> squares of the modulus of elements M+1 to N in that column.
131*> \endverbatim
132*>
133*> \param[in] LDB
134*> \verbatim
135*> LDB is INTEGER
136*> The leading dimension of the array B. LDB >= MAX(1,M,N).
137*> \endverbatim
138*>
139*> \param[out] WORK
140*> \verbatim
141*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
142*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
143*> \endverbatim
144*>
145*> \param[in] LWORK
146*> \verbatim
147*> LWORK is INTEGER
148*> The dimension of the array WORK.
149*> LWORK >= max( 1, MN + max( MN, NRHS ) ).
150*> For optimal performance,
151*> LWORK >= max( 1, (MN + max( MN, NRHS ))*NB ).
152*> where MN = min(M,N) and NB is the optimum block size.
153*>
154*> If LWORK = -1, then a workspace query is assumed; the routine
155*> only calculates the optimal size of the WORK array, returns
156*> this value as the first entry of the WORK array, and no error
157*> message related to LWORK is issued by XERBLA.
158*> \endverbatim
159*>
160*> \param[out] INFO
161*> \verbatim
162*> INFO is INTEGER
163*> = 0: successful exit
164*> < 0: if INFO = -i, the i-th argument had an illegal value
165*> > 0: if INFO = i, the i-th diagonal element of the
166*> triangular factor of A is zero, so that A does not have
167*> full rank; the least squares solution could not be
168*> computed.
169*> \endverbatim
170*
171* Authors:
172* ========
173*
174*> \author Univ. of Tennessee
175*> \author Univ. of California Berkeley
176*> \author Univ. of Colorado Denver
177*> \author NAG Ltd.
178*
179*> \ingroup complex16GEsolve
180*
181*> \par Contributors:
182* ==================
183*>
184*> \verbatim
185*>
186*> November 2022, Igor Kozachenko,
187*> Computer Science Division,
188*> University of California, Berkeley
189*> \endverbatim
190*
191* =====================================================================
192 SUBROUTINE zgelst( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
193 \$ INFO )
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*
533 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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 zgelst(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
ZGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization ...
Definition: zgelst.f:194
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