LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
cgelsx.f
Go to the documentation of this file.
1 *> \brief <b> CGELSX 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 CGELSX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgelsx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgelsx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgelsx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGELSX( 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 * REAL RCOND
27 * ..
28 * .. Array Arguments ..
29 * INTEGER JPVT( * )
30 * REAL RWORK( * )
31 * COMPLEX 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 CGELSY.
41 *>
42 *> CGELSX 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 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 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 REAL
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 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 REAL 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 *> \date November 2011
180 *
181 *> \ingroup complexGEsolve
182 *
183 * =====================================================================
184  SUBROUTINE cgelsx( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
185  $ work, rwork, info )
186 *
187 * -- LAPACK driver routine (version 3.4.0) --
188 * -- LAPACK is a software package provided by Univ. of Tennessee, --
189 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
190 * November 2011
191 *
192 * .. Scalar Arguments ..
193  INTEGER info, lda, ldb, m, n, nrhs, rank
194  REAL rcond
195 * ..
196 * .. Array Arguments ..
197  INTEGER jpvt( * )
198  REAL rwork( * )
199  COMPLEX a( lda, * ), b( ldb, * ), work( * )
200 * ..
201 *
202 * =====================================================================
203 *
204 * .. Parameters ..
205  INTEGER imax, imin
206  parameter( imax = 1, imin = 2 )
207  REAL zero, one, done, ntdone
208  parameter( zero = 0.0e+0, one = 1.0e+0, done = zero,
209  $ ntdone = one )
210  COMPLEX czero, cone
211  parameter( czero = ( 0.0e+0, 0.0e+0 ),
212  $ cone = ( 1.0e+0, 0.0e+0 ) )
213 * ..
214 * .. Local Scalars ..
215  INTEGER i, iascl, ibscl, ismax, ismin, j, k, mn
216  REAL anrm, bignum, bnrm, smax, smaxpr, smin, sminpr,
217  $ smlnum
218  COMPLEX c1, c2, s1, s2, t1, t2
219 * ..
220 * .. External Subroutines ..
221  EXTERNAL cgeqpf, claic1, clascl, claset, clatzm, ctrsm,
223 * ..
224 * .. External Functions ..
225  REAL clange, slamch
226  EXTERNAL clange, slamch
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC abs, conjg, max, min
230 * ..
231 * .. Executable Statements ..
232 *
233  mn = min( m, n )
234  ismin = mn + 1
235  ismax = 2*mn + 1
236 *
237 * Test the input arguments.
238 *
239  info = 0
240  IF( m.LT.0 ) THEN
241  info = -1
242  ELSE IF( n.LT.0 ) THEN
243  info = -2
244  ELSE IF( nrhs.LT.0 ) THEN
245  info = -3
246  ELSE IF( lda.LT.max( 1, m ) ) THEN
247  info = -5
248  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
249  info = -7
250  END IF
251 *
252  IF( info.NE.0 ) THEN
253  CALL xerbla( 'CGELSX', -info )
254  return
255  END IF
256 *
257 * Quick return if possible
258 *
259  IF( min( m, n, nrhs ).EQ.0 ) THEN
260  rank = 0
261  return
262  END IF
263 *
264 * Get machine parameters
265 *
266  smlnum = slamch( 'S' ) / slamch( 'P' )
267  bignum = one / smlnum
268  CALL slabad( smlnum, bignum )
269 *
270 * Scale A, B if max elements outside range [SMLNUM,BIGNUM]
271 *
272  anrm = clange( 'M', m, n, a, lda, rwork )
273  iascl = 0
274  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
275 *
276 * Scale matrix norm up to SMLNUM
277 *
278  CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
279  iascl = 1
280  ELSE IF( anrm.GT.bignum ) THEN
281 *
282 * Scale matrix norm down to BIGNUM
283 *
284  CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
285  iascl = 2
286  ELSE IF( anrm.EQ.zero ) THEN
287 *
288 * Matrix all zero. Return zero solution.
289 *
290  CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
291  rank = 0
292  go to 100
293  END IF
294 *
295  bnrm = clange( 'M', m, nrhs, b, ldb, rwork )
296  ibscl = 0
297  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
298 *
299 * Scale matrix norm up to SMLNUM
300 *
301  CALL clascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
302  ibscl = 1
303  ELSE IF( bnrm.GT.bignum ) THEN
304 *
305 * Scale matrix norm down to BIGNUM
306 *
307  CALL clascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
308  ibscl = 2
309  END IF
310 *
311 * Compute QR factorization with column pivoting of A:
312 * A * P = Q * R
313 *
314  CALL cgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ), rwork,
315  $ info )
316 *
317 * complex workspace MN+N. Real workspace 2*N. Details of Householder
318 * rotations stored in WORK(1:MN).
319 *
320 * Determine RANK using incremental condition estimation
321 *
322  work( ismin ) = cone
323  work( ismax ) = cone
324  smax = abs( a( 1, 1 ) )
325  smin = smax
326  IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
327  rank = 0
328  CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
329  go to 100
330  ELSE
331  rank = 1
332  END IF
333 *
334  10 continue
335  IF( rank.LT.mn ) THEN
336  i = rank + 1
337  CALL claic1( imin, rank, work( ismin ), smin, a( 1, i ),
338  $ a( i, i ), sminpr, s1, c1 )
339  CALL claic1( imax, rank, work( ismax ), smax, a( 1, i ),
340  $ a( i, i ), smaxpr, s2, c2 )
341 *
342  IF( smaxpr*rcond.LE.sminpr ) THEN
343  DO 20 i = 1, rank
344  work( ismin+i-1 ) = s1*work( ismin+i-1 )
345  work( ismax+i-1 ) = s2*work( ismax+i-1 )
346  20 continue
347  work( ismin+rank ) = c1
348  work( ismax+rank ) = c2
349  smin = sminpr
350  smax = smaxpr
351  rank = rank + 1
352  go to 10
353  END IF
354  END IF
355 *
356 * Logically partition R = [ R11 R12 ]
357 * [ 0 R22 ]
358 * where R11 = R(1:RANK,1:RANK)
359 *
360 * [R11,R12] = [ T11, 0 ] * Y
361 *
362  IF( rank.LT.n )
363  $ CALL ctzrqf( rank, n, a, lda, work( mn+1 ), info )
364 *
365 * Details of Householder rotations stored in WORK(MN+1:2*MN)
366 *
367 * B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
368 *
369  CALL cunm2r( 'Left', 'Conjugate transpose', m, nrhs, mn, a, lda,
370  $ work( 1 ), b, ldb, work( 2*mn+1 ), info )
371 *
372 * workspace NRHS
373 *
374 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
375 *
376  CALL ctrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
377  $ nrhs, cone, a, lda, b, ldb )
378 *
379  DO 40 i = rank + 1, n
380  DO 30 j = 1, nrhs
381  b( i, j ) = czero
382  30 continue
383  40 continue
384 *
385 * B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
386 *
387  IF( rank.LT.n ) THEN
388  DO 50 i = 1, rank
389  CALL clatzm( 'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
390  $ conjg( work( mn+i ) ), b( i, 1 ),
391  $ b( rank+1, 1 ), ldb, work( 2*mn+1 ) )
392  50 continue
393  END IF
394 *
395 * workspace NRHS
396 *
397 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
398 *
399  DO 90 j = 1, nrhs
400  DO 60 i = 1, n
401  work( 2*mn+i ) = ntdone
402  60 continue
403  DO 80 i = 1, n
404  IF( work( 2*mn+i ).EQ.ntdone ) THEN
405  IF( jpvt( i ).NE.i ) THEN
406  k = i
407  t1 = b( k, j )
408  t2 = b( jpvt( k ), j )
409  70 continue
410  b( jpvt( k ), j ) = t1
411  work( 2*mn+k ) = done
412  t1 = t2
413  k = jpvt( k )
414  t2 = b( jpvt( k ), j )
415  IF( jpvt( k ).NE.i )
416  $ go to 70
417  b( i, j ) = t1
418  work( 2*mn+k ) = done
419  END IF
420  END IF
421  80 continue
422  90 continue
423 *
424 * Undo scaling
425 *
426  IF( iascl.EQ.1 ) THEN
427  CALL clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
428  CALL clascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
429  $ info )
430  ELSE IF( iascl.EQ.2 ) THEN
431  CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
432  CALL clascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
433  $ info )
434  END IF
435  IF( ibscl.EQ.1 ) THEN
436  CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
437  ELSE IF( ibscl.EQ.2 ) THEN
438  CALL clascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
439  END IF
440 *
441  100 continue
442 *
443  return
444 *
445 * End of CGELSX
446 *
447  END