LAPACK  3.10.0 LAPACK: Linear Algebra PACKage
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
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,
219  \$ ztrsm, ztzrqf, zunm2r
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  CALL dlabad( smlnum, bignum )
266 *
267 * Scale A, B if max elements outside range [SMLNUM,BIGNUM]
268 *
269  anrm = zlange( 'M', m, n, a, lda, rwork )
270  iascl = 0
271  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
272 *
273 * Scale matrix norm up to SMLNUM
274 *
275  CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
276  iascl = 1
277  ELSE IF( anrm.GT.bignum ) THEN
278 *
279 * Scale matrix norm down to BIGNUM
280 *
281  CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
282  iascl = 2
283  ELSE IF( anrm.EQ.zero ) THEN
284 *
285 * Matrix all zero. Return zero solution.
286 *
287  CALL zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
288  rank = 0
289  GO TO 100
290  END IF
291 *
292  bnrm = zlange( 'M', m, nrhs, b, ldb, rwork )
293  ibscl = 0
294  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
295 *
296 * Scale matrix norm up to SMLNUM
297 *
298  CALL zlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
299  ibscl = 1
300  ELSE IF( bnrm.GT.bignum ) THEN
301 *
302 * Scale matrix norm down to BIGNUM
303 *
304  CALL zlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
305  ibscl = 2
306  END IF
307 *
308 * Compute QR factorization with column pivoting of A:
309 * A * P = Q * R
310 *
311  CALL zgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ), rwork,
312  \$ info )
313 *
314 * complex workspace MN+N. Real workspace 2*N. Details of Householder
315 * rotations stored in WORK(1:MN).
316 *
317 * Determine RANK using incremental condition estimation
318 *
319  work( ismin ) = cone
320  work( ismax ) = cone
321  smax = abs( a( 1, 1 ) )
322  smin = smax
323  IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
324  rank = 0
325  CALL zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
326  GO TO 100
327  ELSE
328  rank = 1
329  END IF
330 *
331  10 CONTINUE
332  IF( rank.LT.mn ) THEN
333  i = rank + 1
334  CALL zlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
335  \$ a( i, i ), sminpr, s1, c1 )
336  CALL zlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
337  \$ a( i, i ), smaxpr, s2, c2 )
338 *
339  IF( smaxpr*rcond.LE.sminpr ) THEN
340  DO 20 i = 1, rank
341  work( ismin+i-1 ) = s1*work( ismin+i-1 )
342  work( ismax+i-1 ) = s2*work( ismax+i-1 )
343  20 CONTINUE
344  work( ismin+rank ) = c1
345  work( ismax+rank ) = c2
346  smin = sminpr
347  smax = smaxpr
348  rank = rank + 1
349  GO TO 10
350  END IF
351  END IF
352 *
353 * Logically partition R = [ R11 R12 ]
354 * [ 0 R22 ]
355 * where R11 = R(1:RANK,1:RANK)
356 *
357 * [R11,R12] = [ T11, 0 ] * Y
358 *
359  IF( rank.LT.n )
360  \$ CALL ztzrqf( rank, n, a, lda, work( mn+1 ), info )
361 *
362 * Details of Householder rotations stored in WORK(MN+1:2*MN)
363 *
364 * B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
365 *
366  CALL zunm2r( 'Left', 'Conjugate transpose', m, nrhs, mn, a, lda,
367  \$ work( 1 ), b, ldb, work( 2*mn+1 ), info )
368 *
369 * workspace NRHS
370 *
371 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
372 *
373  CALL ztrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
374  \$ nrhs, cone, a, lda, b, ldb )
375 *
376  DO 40 i = rank + 1, n
377  DO 30 j = 1, nrhs
378  b( i, j ) = czero
379  30 CONTINUE
380  40 CONTINUE
381 *
382 * B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
383 *
384  IF( rank.LT.n ) THEN
385  DO 50 i = 1, rank
386  CALL zlatzm( 'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
387  \$ dconjg( work( mn+i ) ), b( i, 1 ),
388  \$ b( rank+1, 1 ), ldb, work( 2*mn+1 ) )
389  50 CONTINUE
390  END IF
391 *
392 * workspace NRHS
393 *
394 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
395 *
396  DO 90 j = 1, nrhs
397  DO 60 i = 1, n
398  work( 2*mn+i ) = ntdone
399  60 CONTINUE
400  DO 80 i = 1, n
401  IF( work( 2*mn+i ).EQ.ntdone ) THEN
402  IF( jpvt( i ).NE.i ) THEN
403  k = i
404  t1 = b( k, j )
405  t2 = b( jpvt( k ), j )
406  70 CONTINUE
407  b( jpvt( k ), j ) = t1
408  work( 2*mn+k ) = done
409  t1 = t2
410  k = jpvt( k )
411  t2 = b( jpvt( k ), j )
412  IF( jpvt( k ).NE.i )
413  \$ GO TO 70
414  b( i, j ) = t1
415  work( 2*mn+k ) = done
416  END IF
417  END IF
418  80 CONTINUE
419  90 CONTINUE
420 *
421 * Undo scaling
422 *
423  IF( iascl.EQ.1 ) THEN
424  CALL zlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
425  CALL zlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
426  \$ info )
427  ELSE IF( iascl.EQ.2 ) THEN
428  CALL zlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
429  CALL zlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
430  \$ info )
431  END IF
432  IF( ibscl.EQ.1 ) THEN
433  CALL zlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
434  ELSE IF( ibscl.EQ.2 ) THEN
435  CALL zlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
436  END IF
437 *
438  100 CONTINUE
439 *
440  RETURN
441 *
442 * End of ZGELSX
443 *
444  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:180
subroutine zgeqpf(M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO)
ZGEQPF
Definition: zgeqpf.f:148
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 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 zlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
ZLAIC1 applies one step of incremental condition estimation.
Definition: zlaic1.f:135
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 ztzrqf(M, N, A, LDA, TAU, INFO)
ZTZRQF
Definition: ztzrqf.f:138
subroutine zlatzm(SIDE, M, N, V, INCV, TAU, C1, C2, LDC, WORK)
ZLATZM
Definition: zlatzm.f:152
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