LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zgelsx.f
Go to the documentation of this file.
1*> \brief <b> ZGELSX 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 ZGELSX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelsx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelsx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelsx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
22* WORK, RWORK, INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INFO, LDA, LDB, 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*> This routine is deprecated and has been replaced by routine ZGELSY.
41*>
42*> ZGELSX computes the minimum-norm solution to a complex linear least
43*> squares problem:
44*> minimize || A * X - B ||
45*> using a complete orthogonal factorization of A. A is an M-by-N
46*> matrix which may be rank-deficient.
47*>
48*> Several right hand side vectors b and solution vectors x can be
49*> handled in a single call; they are stored as the columns of the
50*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
51*> matrix X.
52*>
53*> The routine first computes a QR factorization with column pivoting:
54*> A * P = Q * [ R11 R12 ]
55*> [ 0 R22 ]
56*> with R11 defined as the largest leading submatrix whose estimated
57*> condition number is less than 1/RCOND. The order of R11, RANK,
58*> is the effective rank of A.
59*>
60*> Then, R22 is considered to be negligible, and R12 is annihilated
61*> by unitary transformations from the right, arriving at the
62*> complete orthogonal factorization:
63*> A * P = Q * [ T11 0 ] * Z
64*> [ 0 0 ]
65*> The minimum-norm solution is then
66*> X = P * Z**H [ inv(T11)*Q1**H*B ]
67*> [ 0 ]
68*> where Q1 consists of the first RANK columns of Q.
69*> \endverbatim
70*
71* Arguments:
72* ==========
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 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*> On exit, A has been overwritten by details of its
98*> complete orthogonal factorization.
99*> \endverbatim
100*>
101*> \param[in] LDA
102*> \verbatim
103*> LDA is INTEGER
104*> The leading dimension of the array A. LDA >= max(1,M).
105*> \endverbatim
106*>
107*> \param[in,out] B
108*> \verbatim
109*> B is COMPLEX*16 array, dimension (LDB,NRHS)
110*> On entry, the M-by-NRHS right hand side matrix B.
111*> On exit, the N-by-NRHS solution matrix X.
112*> If m >= n and RANK = n, the residual sum-of-squares for
113*> the solution in the i-th column is given by the sum of
114*> squares of elements N+1:M in that column.
115*> \endverbatim
116*>
117*> \param[in] LDB
118*> \verbatim
119*> LDB is INTEGER
120*> The leading dimension of the array B. LDB >= max(1,M,N).
121*> \endverbatim
122*>
123*> \param[in,out] JPVT
124*> \verbatim
125*> JPVT is INTEGER array, dimension (N)
126*> On entry, if JPVT(i) .ne. 0, the i-th column of A is an
127*> initial column, otherwise it is a free column. Before
128*> the QR factorization of A, all initial columns are
129*> permuted to the leading positions; only the remaining
130*> free columns are moved as a result of column pivoting
131*> during the factorization.
132*> On exit, if JPVT(i) = k, then the i-th column of A*P
133*> was the k-th column of A.
134*> \endverbatim
135*>
136*> \param[in] RCOND
137*> \verbatim
138*> RCOND is DOUBLE PRECISION
139*> RCOND is used to determine the effective rank of A, which
140*> is defined as the order of the largest leading triangular
141*> submatrix R11 in the QR factorization with pivoting of A,
142*> whose estimated condition number < 1/RCOND.
143*> \endverbatim
144*>
145*> \param[out] RANK
146*> \verbatim
147*> RANK is INTEGER
148*> The effective rank of A, i.e., the order of the submatrix
149*> R11. This is the same as the order of the submatrix T11
150*> in the complete orthogonal factorization of A.
151*> \endverbatim
152*>
153*> \param[out] WORK
154*> \verbatim
155*> WORK is COMPLEX*16 array, dimension
156*> (min(M,N) + max( N, 2*min(M,N)+NRHS )),
157*> \endverbatim
158*>
159*> \param[out] RWORK
160*> \verbatim
161*> RWORK is DOUBLE PRECISION array, dimension (2*N)
162*> \endverbatim
163*>
164*> \param[out] INFO
165*> \verbatim
166*> INFO is INTEGER
167*> = 0: successful exit
168*> < 0: if INFO = -i, the i-th argument had an illegal value
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* =====================================================================
182 SUBROUTINE zgelsx( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
183 $ WORK, RWORK, INFO )
184*
185* -- LAPACK driver routine --
186* -- LAPACK is a software package provided by Univ. of Tennessee, --
187* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188*
189* .. Scalar Arguments ..
190 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
191 DOUBLE PRECISION RCOND
192* ..
193* .. Array Arguments ..
194 INTEGER JPVT( * )
195 DOUBLE PRECISION RWORK( * )
196 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
197* ..
198*
199* =====================================================================
200*
201* .. Parameters ..
202 INTEGER IMAX, IMIN
203 parameter( imax = 1, imin = 2 )
204 DOUBLE PRECISION ZERO, ONE, DONE, NTDONE
205 parameter( zero = 0.0d+0, one = 1.0d+0, done = zero,
206 $ ntdone = one )
207 COMPLEX*16 CZERO, CONE
208 parameter( czero = ( 0.0d+0, 0.0d+0 ),
209 $ cone = ( 1.0d+0, 0.0d+0 ) )
210* ..
211* .. Local Scalars ..
212 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
213 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
214 $ smlnum
215 COMPLEX*16 C1, C2, S1, S2, T1, T2
216* ..
217* .. External Subroutines ..
218 EXTERNAL xerbla, zgeqpf, zlaic1, zlascl, zlaset, zlatzm,
220* ..
221* .. External Functions ..
222 DOUBLE PRECISION DLAMCH, ZLANGE
223 EXTERNAL dlamch, zlange
224* ..
225* .. Intrinsic Functions ..
226 INTRINSIC abs, dconjg, max, min
227* ..
228* .. Executable Statements ..
229*
230 mn = min( m, n )
231 ismin = mn + 1
232 ismax = 2*mn + 1
233*
234* Test the input arguments.
235*
236 info = 0
237 IF( m.LT.0 ) THEN
238 info = -1
239 ELSE IF( n.LT.0 ) THEN
240 info = -2
241 ELSE IF( nrhs.LT.0 ) THEN
242 info = -3
243 ELSE IF( lda.LT.max( 1, m ) ) THEN
244 info = -5
245 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
246 info = -7
247 END IF
248*
249 IF( info.NE.0 ) THEN
250 CALL xerbla( 'ZGELSX', -info )
251 RETURN
252 END IF
253*
254* Quick return if possible
255*
256 IF( min( m, n, nrhs ).EQ.0 ) THEN
257 rank = 0
258 RETURN
259 END IF
260*
261* Get machine parameters
262*
263 smlnum = dlamch( 'S' ) / dlamch( 'P' )
264 bignum = one / smlnum
265*
266* Scale A, B if max elements outside range [SMLNUM,BIGNUM]
267*
268 anrm = zlange( 'M', m, n, a, lda, rwork )
269 iascl = 0
270 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
271*
272* Scale matrix norm up to SMLNUM
273*
274 CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
275 iascl = 1
276 ELSE IF( anrm.GT.bignum ) THEN
277*
278* Scale matrix norm down to BIGNUM
279*
280 CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
281 iascl = 2
282 ELSE IF( anrm.EQ.zero ) THEN
283*
284* Matrix all zero. Return zero solution.
285*
286 CALL zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
287 rank = 0
288 GO TO 100
289 END IF
290*
291 bnrm = zlange( 'M', m, nrhs, b, ldb, rwork )
292 ibscl = 0
293 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
294*
295* Scale matrix norm up to SMLNUM
296*
297 CALL zlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
298 ibscl = 1
299 ELSE IF( bnrm.GT.bignum ) THEN
300*
301* Scale matrix norm down to BIGNUM
302*
303 CALL zlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
304 ibscl = 2
305 END IF
306*
307* Compute QR factorization with column pivoting of A:
308* A * P = Q * R
309*
310 CALL zgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ), rwork,
311 $ info )
312*
313* complex workspace MN+N. Real workspace 2*N. Details of Householder
314* rotations stored in WORK(1:MN).
315*
316* Determine RANK using incremental condition estimation
317*
318 work( ismin ) = cone
319 work( ismax ) = cone
320 smax = abs( a( 1, 1 ) )
321 smin = smax
322 IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
323 rank = 0
324 CALL zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
325 GO TO 100
326 ELSE
327 rank = 1
328 END IF
329*
330 10 CONTINUE
331 IF( rank.LT.mn ) THEN
332 i = rank + 1
333 CALL zlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
334 $ a( i, i ), sminpr, s1, c1 )
335 CALL zlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
336 $ a( i, i ), smaxpr, s2, c2 )
337*
338 IF( smaxpr*rcond.LE.sminpr ) THEN
339 DO 20 i = 1, rank
340 work( ismin+i-1 ) = s1*work( ismin+i-1 )
341 work( ismax+i-1 ) = s2*work( ismax+i-1 )
342 20 CONTINUE
343 work( ismin+rank ) = c1
344 work( ismax+rank ) = c2
345 smin = sminpr
346 smax = smaxpr
347 rank = rank + 1
348 GO TO 10
349 END IF
350 END IF
351*
352* Logically partition R = [ R11 R12 ]
353* [ 0 R22 ]
354* where R11 = R(1:RANK,1:RANK)
355*
356* [R11,R12] = [ T11, 0 ] * Y
357*
358 IF( rank.LT.n )
359 $ CALL ztzrqf( rank, n, a, lda, work( mn+1 ), info )
360*
361* Details of Householder rotations stored in WORK(MN+1:2*MN)
362*
363* B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
364*
365 CALL zunm2r( 'Left', 'Conjugate transpose', m, nrhs, mn, a, lda,
366 $ work( 1 ), b, ldb, work( 2*mn+1 ), info )
367*
368* workspace NRHS
369*
370* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
371*
372 CALL ztrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
373 $ nrhs, cone, a, lda, b, ldb )
374*
375 DO 40 i = rank + 1, n
376 DO 30 j = 1, nrhs
377 b( i, j ) = czero
378 30 CONTINUE
379 40 CONTINUE
380*
381* B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
382*
383 IF( rank.LT.n ) THEN
384 DO 50 i = 1, rank
385 CALL zlatzm( 'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
386 $ dconjg( work( mn+i ) ), b( i, 1 ),
387 $ b( rank+1, 1 ), ldb, work( 2*mn+1 ) )
388 50 CONTINUE
389 END IF
390*
391* workspace NRHS
392*
393* B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
394*
395 DO 90 j = 1, nrhs
396 DO 60 i = 1, n
397 work( 2*mn+i ) = ntdone
398 60 CONTINUE
399 DO 80 i = 1, n
400 IF( work( 2*mn+i ).EQ.ntdone ) THEN
401 IF( jpvt( i ).NE.i ) THEN
402 k = i
403 t1 = b( k, j )
404 t2 = b( jpvt( k ), j )
405 70 CONTINUE
406 b( jpvt( k ), j ) = t1
407 work( 2*mn+k ) = done
408 t1 = t2
409 k = jpvt( k )
410 t2 = b( jpvt( k ), j )
411 IF( jpvt( k ).NE.i )
412 $ GO TO 70
413 b( i, j ) = t1
414 work( 2*mn+k ) = done
415 END IF
416 END IF
417 80 CONTINUE
418 90 CONTINUE
419*
420* Undo scaling
421*
422 IF( iascl.EQ.1 ) THEN
423 CALL zlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
424 CALL zlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
425 $ info )
426 ELSE IF( iascl.EQ.2 ) THEN
427 CALL zlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
428 CALL zlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
429 $ info )
430 END IF
431 IF( ibscl.EQ.1 ) THEN
432 CALL zlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
433 ELSE IF( ibscl.EQ.2 ) THEN
434 CALL zlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
435 END IF
436*
437 100 CONTINUE
438*
439 RETURN
440*
441* End of ZGELSX
442*
443 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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 zunm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
ZUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition zunm2r.f:159
subroutine zgelsx(m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, rwork, info)
ZGELSX solves overdetermined or underdetermined systems for GE matrices
Definition zgelsx.f:184
subroutine zgeqpf(m, n, a, lda, jpvt, tau, work, rwork, info)
ZGEQPF
Definition zgeqpf.f:148
subroutine zlatzm(side, m, n, v, incv, tau, c1, c2, ldc, work)
ZLATZM
Definition zlatzm.f:152
subroutine ztzrqf(m, n, a, lda, tau, info)
ZTZRQF
Definition ztzrqf.f:138