LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zgelsy.f
Go to the documentation of this file.
1*> \brief <b> ZGELSY 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 ZGELSY + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelsy.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelsy.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelsy.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
22* WORK, LWORK, RWORK, INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
26* DOUBLE PRECISION RCOND
27* ..
28* .. Array Arguments ..
29* INTEGER JPVT( * )
30* DOUBLE PRECISION RWORK( * )
31* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> ZGELSY computes the minimum-norm solution to a complex linear least
41*> squares problem:
42*> minimize || A * X - B ||
43*> using a complete orthogonal factorization of A. A is an M-by-N
44*> matrix which may be rank-deficient.
45*>
46*> Several right hand side vectors b and solution vectors x can be
47*> handled in a single call; they are stored as the columns of the
48*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
49*> matrix X.
50*>
51*> The routine first computes a QR factorization with column pivoting:
52*> A * P = Q * [ R11 R12 ]
53*> [ 0 R22 ]
54*> with R11 defined as the largest leading submatrix whose estimated
55*> condition number is less than 1/RCOND. The order of R11, RANK,
56*> is the effective rank of A.
57*>
58*> Then, R22 is considered to be negligible, and R12 is annihilated
59*> by unitary transformations from the right, arriving at the
60*> complete orthogonal factorization:
61*> A * P = Q * [ T11 0 ] * Z
62*> [ 0 0 ]
63*> The minimum-norm solution is then
64*> X = P * Z**H [ inv(T11)*Q1**H*B ]
65*> [ 0 ]
66*> where Q1 consists of the first RANK columns of Q.
67*>
68*> This routine is basically identical to the original xGELSX except
69*> three differences:
70*> o The permutation of matrix B (the right hand side) is faster and
71*> more simple.
72*> o The call to the subroutine xGEQPF has been substituted by the
73*> the call to the subroutine xGEQP3. This subroutine is a Blas-3
74*> version of the QR factorization with column pivoting.
75*> o Matrix B (the right hand side) is updated with Blas-3.
76*> \endverbatim
77*
78* Arguments:
79* ==========
80*
81*> \param[in] M
82*> \verbatim
83*> M is INTEGER
84*> The number of rows of the matrix A. M >= 0.
85*> \endverbatim
86*>
87*> \param[in] N
88*> \verbatim
89*> N is INTEGER
90*> The number of columns of the matrix A. N >= 0.
91*> \endverbatim
92*>
93*> \param[in] NRHS
94*> \verbatim
95*> NRHS is INTEGER
96*> The number of right hand sides, i.e., the number of
97*> columns of matrices B and X. NRHS >= 0.
98*> \endverbatim
99*>
100*> \param[in,out] A
101*> \verbatim
102*> A is COMPLEX*16 array, dimension (LDA,N)
103*> On entry, the M-by-N matrix A.
104*> On exit, A has been overwritten by details of its
105*> complete orthogonal factorization.
106*> \endverbatim
107*>
108*> \param[in] LDA
109*> \verbatim
110*> LDA is INTEGER
111*> The leading dimension of the array A. LDA >= max(1,M).
112*> \endverbatim
113*>
114*> \param[in,out] B
115*> \verbatim
116*> B is COMPLEX*16 array, dimension (LDB,NRHS)
117*> On entry, the M-by-NRHS right hand side matrix B.
118*> On exit, the N-by-NRHS solution matrix X.
119*> If M = 0 or N = 0, B is not referenced.
120*> \endverbatim
121*>
122*> \param[in] LDB
123*> \verbatim
124*> LDB is INTEGER
125*> The leading dimension of the array B. LDB >= max(1,M,N).
126*> \endverbatim
127*>
128*> \param[in,out] JPVT
129*> \verbatim
130*> JPVT is INTEGER array, dimension (N)
131*> On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
132*> to the front of AP, otherwise column i is a free column.
133*> On exit, if JPVT(i) = k, then the i-th column of A*P
134*> was the k-th column of A.
135*> \endverbatim
136*>
137*> \param[in] RCOND
138*> \verbatim
139*> RCOND is DOUBLE PRECISION
140*> RCOND is used to determine the effective rank of A, which
141*> is defined as the order of the largest leading triangular
142*> submatrix R11 in the QR factorization with pivoting of A,
143*> whose estimated condition number < 1/RCOND.
144*> \endverbatim
145*>
146*> \param[out] RANK
147*> \verbatim
148*> RANK is INTEGER
149*> The effective rank of A, i.e., the order of the submatrix
150*> R11. This is the same as the order of the submatrix T11
151*> in the complete orthogonal factorization of A.
152*> If NRHS = 0, RANK = 0 on output.
153*> \endverbatim
154*>
155*> \param[out] WORK
156*> \verbatim
157*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
158*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
159*> \endverbatim
160*>
161*> \param[in] LWORK
162*> \verbatim
163*> LWORK is INTEGER
164*> The dimension of the array WORK.
165*> The unblocked strategy requires that:
166*> LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
167*> where MN = min(M,N).
168*> The block algorithm requires that:
169*> LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
170*> where NB is an upper bound on the blocksize returned
171*> by ILAENV for the routines ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR,
172*> and ZUNMRZ.
173*>
174*> If LWORK = -1, then a workspace query is assumed; the routine
175*> only calculates the optimal size of the WORK array, returns
176*> this value as the first entry of the WORK array, and no error
177*> message related to LWORK is issued by XERBLA.
178*> \endverbatim
179*>
180*> \param[out] RWORK
181*> \verbatim
182*> RWORK is DOUBLE PRECISION array, dimension (2*N)
183*> \endverbatim
184*>
185*> \param[out] INFO
186*> \verbatim
187*> INFO is INTEGER
188*> = 0: successful exit
189*> < 0: if INFO = -i, the i-th argument had an illegal value
190*> \endverbatim
191*
192* Authors:
193* ========
194*
195*> \author Univ. of Tennessee
196*> \author Univ. of California Berkeley
197*> \author Univ. of Colorado Denver
198*> \author NAG Ltd.
199*
200*> \ingroup gelsy
201*
202*> \par Contributors:
203* ==================
204*>
205*> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n
206*> E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
207*> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
208*>
209* =====================================================================
210 SUBROUTINE zgelsy( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
211 $ WORK, LWORK, RWORK, INFO )
212*
213* -- LAPACK driver routine --
214* -- LAPACK is a software package provided by Univ. of Tennessee, --
215* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216*
217* .. Scalar Arguments ..
218 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
219 DOUBLE PRECISION RCOND
220* ..
221* .. Array Arguments ..
222 INTEGER JPVT( * )
223 DOUBLE PRECISION RWORK( * )
224 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
225* ..
226*
227* =====================================================================
228*
229* .. Parameters ..
230 INTEGER IMAX, IMIN
231 parameter( imax = 1, imin = 2 )
232 DOUBLE PRECISION ZERO, ONE
233 parameter( zero = 0.0d+0, one = 1.0d+0 )
234 COMPLEX*16 CZERO, CONE
235 parameter( czero = ( 0.0d+0, 0.0d+0 ),
236 $ cone = ( 1.0d+0, 0.0d+0 ) )
237* ..
238* .. Local Scalars ..
239 LOGICAL LQUERY
240 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN,
241 $ nb, nb1, nb2, nb3, nb4
242 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
243 $ smlnum, wsize
244 COMPLEX*16 C1, C2, S1, S2
245* ..
246* .. External Subroutines ..
247 EXTERNAL xerbla, zcopy, zgeqp3, zlaic1, zlascl,
249* ..
250* .. External Functions ..
251 INTEGER ILAENV
252 DOUBLE PRECISION DLAMCH, ZLANGE
253 EXTERNAL ilaenv, dlamch, zlange
254* ..
255* .. Intrinsic Functions ..
256 INTRINSIC abs, dble, dcmplx, max, min
257* ..
258* .. Executable Statements ..
259*
260 mn = min( m, n )
261 ismin = mn + 1
262 ismax = 2*mn + 1
263*
264* Test the input arguments.
265*
266 info = 0
267 nb1 = ilaenv( 1, 'ZGEQRF', ' ', m, n, -1, -1 )
268 nb2 = ilaenv( 1, 'ZGERQF', ' ', m, n, -1, -1 )
269 nb3 = ilaenv( 1, 'ZUNMQR', ' ', m, n, nrhs, -1 )
270 nb4 = ilaenv( 1, 'ZUNMRQ', ' ', m, n, nrhs, -1 )
271 nb = max( nb1, nb2, nb3, nb4 )
272 lwkopt = max( 1, mn+2*n+nb*( n+1 ), 2*mn+nb*nrhs )
273 work( 1 ) = dcmplx( lwkopt )
274 lquery = ( lwork.EQ.-1 )
275 IF( m.LT.0 ) THEN
276 info = -1
277 ELSE IF( n.LT.0 ) THEN
278 info = -2
279 ELSE IF( nrhs.LT.0 ) THEN
280 info = -3
281 ELSE IF( lda.LT.max( 1, m ) ) THEN
282 info = -5
283 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
284 info = -7
285 ELSE IF( lwork.LT.( mn+max( 2*mn, n+1, mn+nrhs ) ) .AND. .NOT.
286 $ lquery ) THEN
287 info = -12
288 END IF
289*
290 IF( info.NE.0 ) THEN
291 CALL xerbla( 'ZGELSY', -info )
292 RETURN
293 ELSE IF( lquery ) THEN
294 RETURN
295 END IF
296*
297* Quick return if possible
298*
299 IF( min( m, n, nrhs ).EQ.0 ) THEN
300 rank = 0
301 RETURN
302 END IF
303*
304* Get machine parameters
305*
306 smlnum = dlamch( 'S' ) / dlamch( 'P' )
307 bignum = one / smlnum
308*
309* Scale A, B if max entries outside range [SMLNUM,BIGNUM]
310*
311 anrm = zlange( 'M', m, n, a, lda, rwork )
312 iascl = 0
313 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
314*
315* Scale matrix norm up to SMLNUM
316*
317 CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
318 iascl = 1
319 ELSE IF( anrm.GT.bignum ) THEN
320*
321* Scale matrix norm down to BIGNUM
322*
323 CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
324 iascl = 2
325 ELSE IF( anrm.EQ.zero ) THEN
326*
327* Matrix all zero. Return zero solution.
328*
329 CALL zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
330 rank = 0
331 GO TO 70
332 END IF
333*
334 bnrm = zlange( 'M', m, nrhs, b, ldb, rwork )
335 ibscl = 0
336 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
337*
338* Scale matrix norm up to SMLNUM
339*
340 CALL zlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
341 ibscl = 1
342 ELSE IF( bnrm.GT.bignum ) THEN
343*
344* Scale matrix norm down to BIGNUM
345*
346 CALL zlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
347 ibscl = 2
348 END IF
349*
350* Compute QR factorization with column pivoting of A:
351* A * P = Q * R
352*
353 CALL zgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
354 $ lwork-mn, rwork, info )
355 wsize = mn + dble( work( mn+1 ) )
356*
357* complex workspace: MN+NB*(N+1). real workspace 2*N.
358* Details of Householder rotations stored in WORK(1:MN).
359*
360* Determine RANK using incremental condition estimation
361*
362 work( ismin ) = cone
363 work( ismax ) = cone
364 smax = abs( a( 1, 1 ) )
365 smin = smax
366 IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
367 rank = 0
368 CALL zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
369 GO TO 70
370 ELSE
371 rank = 1
372 END IF
373*
374 10 CONTINUE
375 IF( rank.LT.mn ) THEN
376 i = rank + 1
377 CALL zlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
378 $ a( i, i ), sminpr, s1, c1 )
379 CALL zlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
380 $ a( i, i ), smaxpr, s2, c2 )
381*
382 IF( smaxpr*rcond.LE.sminpr ) THEN
383 DO 20 i = 1, rank
384 work( ismin+i-1 ) = s1*work( ismin+i-1 )
385 work( ismax+i-1 ) = s2*work( ismax+i-1 )
386 20 CONTINUE
387 work( ismin+rank ) = c1
388 work( ismax+rank ) = c2
389 smin = sminpr
390 smax = smaxpr
391 rank = rank + 1
392 GO TO 10
393 END IF
394 END IF
395*
396* complex workspace: 3*MN.
397*
398* Logically partition R = [ R11 R12 ]
399* [ 0 R22 ]
400* where R11 = R(1:RANK,1:RANK)
401*
402* [R11,R12] = [ T11, 0 ] * Y
403*
404 IF( rank.LT.n )
405 $ CALL ztzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
406 $ lwork-2*mn, info )
407*
408* complex workspace: 2*MN.
409* Details of Householder rotations stored in WORK(MN+1:2*MN)
410*
411* B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
412*
413 CALL zunmqr( 'Left', 'Conjugate transpose', m, nrhs, mn, a, lda,
414 $ work( 1 ), b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
415 wsize = max( wsize, 2*mn+dble( work( 2*mn+1 ) ) )
416*
417* complex workspace: 2*MN+NB*NRHS.
418*
419* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
420*
421 CALL ztrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
422 $ nrhs, cone, a, lda, b, ldb )
423*
424 DO 40 j = 1, nrhs
425 DO 30 i = rank + 1, n
426 b( i, j ) = czero
427 30 CONTINUE
428 40 CONTINUE
429*
430* B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
431*
432 IF( rank.LT.n ) THEN
433 CALL zunmrz( 'Left', 'Conjugate transpose', n, nrhs, rank,
434 $ n-rank, a, lda, work( mn+1 ), b, ldb,
435 $ work( 2*mn+1 ), lwork-2*mn, info )
436 END IF
437*
438* complex workspace: 2*MN+NRHS.
439*
440* B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
441*
442 DO 60 j = 1, nrhs
443 DO 50 i = 1, n
444 work( jpvt( i ) ) = b( i, j )
445 50 CONTINUE
446 CALL zcopy( n, work( 1 ), 1, b( 1, j ), 1 )
447 60 CONTINUE
448*
449* complex workspace: N.
450*
451* Undo scaling
452*
453 IF( iascl.EQ.1 ) THEN
454 CALL zlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
455 CALL zlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
456 $ info )
457 ELSE IF( iascl.EQ.2 ) THEN
458 CALL zlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
459 CALL zlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
460 $ info )
461 END IF
462 IF( ibscl.EQ.1 ) THEN
463 CALL zlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
464 ELSE IF( ibscl.EQ.2 ) THEN
465 CALL zlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
466 END IF
467*
468 70 CONTINUE
469 work( 1 ) = dcmplx( lwkopt )
470*
471 RETURN
472*
473* End of ZGELSY
474*
475 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
subroutine zgelsy(m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, rwork, info)
ZGELSY solves overdetermined or underdetermined systems for GE matrices
Definition zgelsy.f:212
subroutine zgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info)
ZGEQP3
Definition zgeqp3.f:159
subroutine zlaic1(job, j, x, sest, w, gamma, sestpr, s, c)
ZLAIC1 applies one step of incremental condition estimation.
Definition zlaic1.f:135
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 ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM
Definition ztrsm.f:180
subroutine ztzrzf(m, n, a, lda, tau, work, lwork, info)
ZTZRZF
Definition ztzrzf.f:151
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR
Definition zunmqr.f:167
subroutine zunmrz(side, trans, m, n, k, l, a, lda, tau, c, ldc, work, lwork, info)
ZUNMRZ
Definition zunmrz.f:187