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