LAPACK  3.10.0 LAPACK: Linear Algebra PACKage
cgelsy.f
Go to the documentation of this file.
1 *> \brief <b> CGELSY 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/cgelsy.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgelsy.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgelsy.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGELSY( 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 * 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 *> CGELSY 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 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 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 *> \endverbatim
120 *>
121 *> \param[in] LDB
122 *> \verbatim
123 *> LDB is INTEGER
124 *> The leading dimension of the array B. LDB >= max(1,M,N).
125 *> \endverbatim
126 *>
127 *> \param[in,out] JPVT
128 *> \verbatim
129 *> JPVT is INTEGER array, dimension (N)
130 *> On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
131 *> to the front of AP, otherwise column i is a free column.
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 (MAX(1,LWORK))
156 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
157 *> \endverbatim
158 *>
159 *> \param[in] LWORK
160 *> \verbatim
161 *> LWORK is INTEGER
162 *> The dimension of the array WORK.
163 *> The unblocked strategy requires that:
164 *> LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
165 *> where MN = min(M,N).
166 *> The block algorithm requires that:
167 *> LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
168 *> where NB is an upper bound on the blocksize returned
169 *> by ILAENV for the routines CGEQP3, CTZRZF, CTZRQF, CUNMQR,
170 *> and CUNMRZ.
171 *>
172 *> If LWORK = -1, then a workspace query is assumed; the routine
173 *> only calculates the optimal size of the WORK array, returns
174 *> this value as the first entry of the WORK array, and no error
175 *> message related to LWORK is issued by XERBLA.
176 *> \endverbatim
177 *>
178 *> \param[out] RWORK
179 *> \verbatim
180 *> RWORK is REAL array, dimension (2*N)
181 *> \endverbatim
182 *>
183 *> \param[out] INFO
184 *> \verbatim
185 *> INFO is INTEGER
186 *> = 0: successful exit
187 *> < 0: if INFO = -i, the i-th argument had an illegal value
188 *> \endverbatim
189 *
190 * Authors:
191 * ========
192 *
193 *> \author Univ. of Tennessee
194 *> \author Univ. of California Berkeley
195 *> \author Univ. of Colorado Denver
196 *> \author NAG Ltd.
197 *
198 *> \ingroup complexGEsolve
199 *
200 *> \par Contributors:
201 * ==================
202 *>
203 *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n
204 *> E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
205 *> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
206 *>
207 * =====================================================================
208  SUBROUTINE cgelsy( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
209  \$ WORK, LWORK, RWORK, INFO )
210 *
211 * -- LAPACK driver routine --
212 * -- LAPACK is a software package provided by Univ. of Tennessee, --
213 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
214 *
215 * .. Scalar Arguments ..
216  INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
217  REAL RCOND
218 * ..
219 * .. Array Arguments ..
220  INTEGER JPVT( * )
221  REAL RWORK( * )
222  COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
223 * ..
224 *
225 * =====================================================================
226 *
227 * .. Parameters ..
228  INTEGER IMAX, IMIN
229  parameter( imax = 1, imin = 2 )
230  REAL ZERO, ONE
231  parameter( zero = 0.0e+0, one = 1.0e+0 )
232  COMPLEX CZERO, CONE
233  parameter( czero = ( 0.0e+0, 0.0e+0 ),
234  \$ cone = ( 1.0e+0, 0.0e+0 ) )
235 * ..
236 * .. Local Scalars ..
237  LOGICAL LQUERY
238  INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN,
239  \$ nb, nb1, nb2, nb3, nb4
240  REAL ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
241  \$ smlnum, wsize
242  COMPLEX C1, C2, S1, S2
243 * ..
244 * .. External Subroutines ..
245  EXTERNAL ccopy, cgeqp3, claic1, clascl, claset, ctrsm,
247 * ..
248 * .. External Functions ..
249  INTEGER ILAENV
250  REAL CLANGE, SLAMCH
251  EXTERNAL clange, ilaenv, slamch
252 * ..
253 * .. Intrinsic Functions ..
254  INTRINSIC abs, max, min, real, cmplx
255 * ..
256 * .. Executable Statements ..
257 *
258  mn = min( m, n )
259  ismin = mn + 1
260  ismax = 2*mn + 1
261 *
262 * Test the input arguments.
263 *
264  info = 0
265  nb1 = ilaenv( 1, 'CGEQRF', ' ', m, n, -1, -1 )
266  nb2 = ilaenv( 1, 'CGERQF', ' ', m, n, -1, -1 )
267  nb3 = ilaenv( 1, 'CUNMQR', ' ', m, n, nrhs, -1 )
268  nb4 = ilaenv( 1, 'CUNMRQ', ' ', m, n, nrhs, -1 )
269  nb = max( nb1, nb2, nb3, nb4 )
270  lwkopt = max( 1, mn+2*n+nb*(n+1), 2*mn+nb*nrhs )
271  work( 1 ) = cmplx( lwkopt )
272  lquery = ( lwork.EQ.-1 )
273  IF( m.LT.0 ) THEN
274  info = -1
275  ELSE IF( n.LT.0 ) THEN
276  info = -2
277  ELSE IF( nrhs.LT.0 ) THEN
278  info = -3
279  ELSE IF( lda.LT.max( 1, m ) ) THEN
280  info = -5
281  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
282  info = -7
283  ELSE IF( lwork.LT.( mn+max( 2*mn, n+1, mn+nrhs ) ) .AND.
284  \$ .NOT.lquery ) THEN
285  info = -12
286  END IF
287 *
288  IF( info.NE.0 ) THEN
289  CALL xerbla( 'CGELSY', -info )
290  RETURN
291  ELSE IF( lquery ) THEN
292  RETURN
293  END IF
294 *
295 * Quick return if possible
296 *
297  IF( min( m, n, nrhs ).EQ.0 ) THEN
298  rank = 0
299  RETURN
300  END IF
301 *
302 * Get machine parameters
303 *
304  smlnum = slamch( 'S' ) / slamch( 'P' )
305  bignum = one / smlnum
306  CALL slabad( smlnum, bignum )
307 *
308 * Scale A, B if max entries outside range [SMLNUM,BIGNUM]
309 *
310  anrm = clange( 'M', m, n, a, lda, rwork )
311  iascl = 0
312  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
313 *
314 * Scale matrix norm up to SMLNUM
315 *
316  CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
317  iascl = 1
318  ELSE IF( anrm.GT.bignum ) THEN
319 *
320 * Scale matrix norm down to BIGNUM
321 *
322  CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
323  iascl = 2
324  ELSE IF( anrm.EQ.zero ) THEN
325 *
326 * Matrix all zero. Return zero solution.
327 *
328  CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
329  rank = 0
330  GO TO 70
331  END IF
332 *
333  bnrm = clange( 'M', m, nrhs, b, ldb, rwork )
334  ibscl = 0
335  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
336 *
337 * Scale matrix norm up to SMLNUM
338 *
339  CALL clascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
340  ibscl = 1
341  ELSE IF( bnrm.GT.bignum ) THEN
342 *
343 * Scale matrix norm down to BIGNUM
344 *
345  CALL clascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
346  ibscl = 2
347  END IF
348 *
349 * Compute QR factorization with column pivoting of A:
350 * A * P = Q * R
351 *
352  CALL cgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
353  \$ lwork-mn, rwork, info )
354  wsize = mn + real( work( mn+1 ) )
355 *
356 * complex workspace: MN+NB*(N+1). real workspace 2*N.
357 * Details of Householder rotations stored in WORK(1:MN).
358 *
359 * Determine RANK using incremental condition estimation
360 *
361  work( ismin ) = cone
362  work( ismax ) = cone
363  smax = abs( a( 1, 1 ) )
364  smin = smax
365  IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
366  rank = 0
367  CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
368  GO TO 70
369  ELSE
370  rank = 1
371  END IF
372 *
373  10 CONTINUE
374  IF( rank.LT.mn ) THEN
375  i = rank + 1
376  CALL claic1( imin, rank, work( ismin ), smin, a( 1, i ),
377  \$ a( i, i ), sminpr, s1, c1 )
378  CALL claic1( imax, rank, work( ismax ), smax, a( 1, i ),
379  \$ a( i, i ), smaxpr, s2, c2 )
380 *
381  IF( smaxpr*rcond.LE.sminpr ) THEN
382  DO 20 i = 1, rank
383  work( ismin+i-1 ) = s1*work( ismin+i-1 )
384  work( ismax+i-1 ) = s2*work( ismax+i-1 )
385  20 CONTINUE
386  work( ismin+rank ) = c1
387  work( ismax+rank ) = c2
388  smin = sminpr
389  smax = smaxpr
390  rank = rank + 1
391  GO TO 10
392  END IF
393  END IF
394 *
395 * complex workspace: 3*MN.
396 *
397 * Logically partition R = [ R11 R12 ]
398 * [ 0 R22 ]
399 * where R11 = R(1:RANK,1:RANK)
400 *
401 * [R11,R12] = [ T11, 0 ] * Y
402 *
403  IF( rank.LT.n )
404  \$ CALL ctzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
405  \$ lwork-2*mn, info )
406 *
407 * complex workspace: 2*MN.
408 * Details of Householder rotations stored in WORK(MN+1:2*MN)
409 *
410 * B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
411 *
412  CALL cunmqr( 'Left', 'Conjugate transpose', m, nrhs, mn, a, lda,
413  \$ work( 1 ), b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
414  wsize = max( wsize, 2*mn+real( work( 2*mn+1 ) ) )
415 *
416 * complex workspace: 2*MN+NB*NRHS.
417 *
418 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
419 *
420  CALL ctrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
421  \$ nrhs, cone, a, lda, b, ldb )
422 *
423  DO 40 j = 1, nrhs
424  DO 30 i = rank + 1, n
425  b( i, j ) = czero
426  30 CONTINUE
427  40 CONTINUE
428 *
429 * B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
430 *
431  IF( rank.LT.n ) THEN
432  CALL cunmrz( 'Left', 'Conjugate transpose', n, nrhs, rank,
433  \$ n-rank, a, lda, work( mn+1 ), b, ldb,
434  \$ work( 2*mn+1 ), lwork-2*mn, info )
435  END IF
436 *
437 * complex workspace: 2*MN+NRHS.
438 *
439 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
440 *
441  DO 60 j = 1, nrhs
442  DO 50 i = 1, n
443  work( jpvt( i ) ) = b( i, j )
444  50 CONTINUE
445  CALL ccopy( n, work( 1 ), 1, b( 1, j ), 1 )
446  60 CONTINUE
447 *
448 * complex workspace: N.
449 *
450 * Undo scaling
451 *
452  IF( iascl.EQ.1 ) THEN
453  CALL clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
454  CALL clascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
455  \$ info )
456  ELSE IF( iascl.EQ.2 ) THEN
457  CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
458  CALL clascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
459  \$ info )
460  END IF
461  IF( ibscl.EQ.1 ) THEN
462  CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
463  ELSE IF( ibscl.EQ.2 ) THEN
464  CALL clascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
465  END IF
466 *
467  70 CONTINUE
468  work( 1 ) = cmplx( lwkopt )
469 *
470  RETURN
471 *
472 * End of CGELSY
473 *
474  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:180
subroutine cgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
CGEQP3
Definition: cgeqp3.f:159
subroutine cgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, RWORK, INFO)
CGELSY solves overdetermined or underdetermined systems for GE matrices
Definition: cgelsy.f:210
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:106
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:143
subroutine claic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
CLAIC1 applies one step of incremental condition estimation.
Definition: claic1.f:135
subroutine ctzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CTZRZF
Definition: ctzrzf.f:151
subroutine cunmrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMRZ
Definition: cunmrz.f:187
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:168