LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ cgelsx()

subroutine cgelsx ( integer m,
integer n,
integer nrhs,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
integer, dimension( * ) jpvt,
real rcond,
integer rank,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer info )

CGELSX solves overdetermined or underdetermined systems for GE matrices

Download CGELSX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> This routine is deprecated and has been replaced by routine CGELSY.
!>
!> CGELSX computes the minimum-norm solution to a complex linear least
!> squares problem:
!>     minimize || A * X - B ||
!> using a complete orthogonal factorization of A.  A is an M-by-N
!> matrix which may be rank-deficient.
!>
!> Several right hand side vectors b and solution vectors x can be
!> handled in a single call; they are stored as the columns of the
!> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
!> matrix X.
!>
!> The routine first computes a QR factorization with column pivoting:
!>     A * P = Q * [ R11 R12 ]
!>                 [  0  R22 ]
!> with R11 defined as the largest leading submatrix whose estimated
!> condition number is less than 1/RCOND.  The order of R11, RANK,
!> is the effective rank of A.
!>
!> Then, R22 is considered to be negligible, and R12 is annihilated
!> by unitary transformations from the right, arriving at the
!> complete orthogonal factorization:
!>    A * P = Q * [ T11 0 ] * Z
!>                [  0  0 ]
!> The minimum-norm solution is then
!>    X = P * Z**H [ inv(T11)*Q1**H*B ]
!>                 [        0         ]
!> where Q1 consists of the first RANK columns of Q.
!>
Parameters
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of
!>          columns of matrices B and X. NRHS >= 0.
!>
[in,out]A
!>          A is COMPLEX array, dimension (LDA,N)
!>          On entry, the M-by-N matrix A.
!>          On exit, A has been overwritten by details of its
!>          complete orthogonal factorization.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[in,out]B
!>          B is COMPLEX array, dimension (LDB,NRHS)
!>          On entry, the M-by-NRHS right hand side matrix B.
!>          On exit, the N-by-NRHS solution matrix X.
!>          If m >= n and RANK = n, the residual sum-of-squares for
!>          the solution in the i-th column is given by the sum of
!>          squares of elements N+1:M in that column.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,M,N).
!>
[in,out]JPVT
!>          JPVT is INTEGER array, dimension (N)
!>          On entry, if JPVT(i) .ne. 0, the i-th column of A is an
!>          initial column, otherwise it is a free column.  Before
!>          the QR factorization of A, all initial columns are
!>          permuted to the leading positions; only the remaining
!>          free columns are moved as a result of column pivoting
!>          during the factorization.
!>          On exit, if JPVT(i) = k, then the i-th column of A*P
!>          was the k-th column of A.
!>
[in]RCOND
!>          RCOND is REAL
!>          RCOND is used to determine the effective rank of A, which
!>          is defined as the order of the largest leading triangular
!>          submatrix R11 in the QR factorization with pivoting of A,
!>          whose estimated condition number < 1/RCOND.
!>
[out]RANK
!>          RANK is INTEGER
!>          The effective rank of A, i.e., the order of the submatrix
!>          R11.  This is the same as the order of the submatrix T11
!>          in the complete orthogonal factorization of A.
!>
[out]WORK
!>          WORK is COMPLEX array, dimension
!>                      (min(M,N) + max( N, 2*min(M,N)+NRHS )),
!>
[out]RWORK
!>          RWORK is REAL array, dimension (2*N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 180 of file cgelsx.f.

182*
183* -- LAPACK driver routine --
184* -- LAPACK is a software package provided by Univ. of Tennessee, --
185* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186*
187* .. Scalar Arguments ..
188 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
189 REAL RCOND
190* ..
191* .. Array Arguments ..
192 INTEGER JPVT( * )
193 REAL RWORK( * )
194 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
195* ..
196*
197* =====================================================================
198*
199* .. Parameters ..
200 INTEGER IMAX, IMIN
201 parameter( imax = 1, imin = 2 )
202 REAL ZERO, ONE, DONE, NTDONE
203 parameter( zero = 0.0e+0, one = 1.0e+0, done = zero,
204 $ ntdone = one )
205 COMPLEX CZERO, CONE
206 parameter( czero = ( 0.0e+0, 0.0e+0 ),
207 $ cone = ( 1.0e+0, 0.0e+0 ) )
208* ..
209* .. Local Scalars ..
210 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
211 REAL ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
212 $ SMLNUM
213 COMPLEX C1, C2, S1, S2, T1, T2
214* ..
215* .. External Subroutines ..
216 EXTERNAL cgeqpf, claic1, clascl, claset, clatzm, ctrsm,
217 $ ctzrqf, cunm2r, xerbla
218* ..
219* .. External Functions ..
220 REAL CLANGE, SLAMCH
221 EXTERNAL clange, slamch
222* ..
223* .. Intrinsic Functions ..
224 INTRINSIC abs, conjg, max, min
225* ..
226* .. Executable Statements ..
227*
228 mn = min( m, n )
229 ismin = mn + 1
230 ismax = 2*mn + 1
231*
232* Test the input arguments.
233*
234 info = 0
235 IF( m.LT.0 ) THEN
236 info = -1
237 ELSE IF( n.LT.0 ) THEN
238 info = -2
239 ELSE IF( nrhs.LT.0 ) THEN
240 info = -3
241 ELSE IF( lda.LT.max( 1, m ) ) THEN
242 info = -5
243 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
244 info = -7
245 END IF
246*
247 IF( info.NE.0 ) THEN
248 CALL xerbla( 'CGELSX', -info )
249 RETURN
250 END IF
251*
252* Quick return if possible
253*
254 IF( min( m, n, nrhs ).EQ.0 ) THEN
255 rank = 0
256 RETURN
257 END IF
258*
259* Get machine parameters
260*
261 smlnum = slamch( 'S' ) / slamch( 'P' )
262 bignum = one / smlnum
263*
264* Scale A, B if max elements outside range [SMLNUM,BIGNUM]
265*
266 anrm = clange( 'M', m, n, a, lda, rwork )
267 iascl = 0
268 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
269*
270* Scale matrix norm up to SMLNUM
271*
272 CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
273 iascl = 1
274 ELSE IF( anrm.GT.bignum ) THEN
275*
276* Scale matrix norm down to BIGNUM
277*
278 CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
279 iascl = 2
280 ELSE IF( anrm.EQ.zero ) THEN
281*
282* Matrix all zero. Return zero solution.
283*
284 CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
285 rank = 0
286 GO TO 100
287 END IF
288*
289 bnrm = clange( 'M', m, nrhs, b, ldb, rwork )
290 ibscl = 0
291 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
292*
293* Scale matrix norm up to SMLNUM
294*
295 CALL clascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb,
296 $ info )
297 ibscl = 1
298 ELSE IF( bnrm.GT.bignum ) THEN
299*
300* Scale matrix norm down to BIGNUM
301*
302 CALL clascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb,
303 $ 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 cgeqpf( 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 claset( '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 claic1( imin, rank, work( ismin ), smin, a( 1, i ),
334 $ a( i, i ), sminpr, s1, c1 )
335 CALL claic1( 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 ctzrqf( 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 cunm2r( 'Left', 'Conjugate transpose', m, nrhs, mn, a,
366 $ lda, 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 ctrsm( '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 clatzm( 'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
386 $ conjg( 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 clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb,
424 $ info )
425 CALL clascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
426 $ info )
427 ELSE IF( iascl.EQ.2 ) THEN
428 CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb,
429 $ info )
430 CALL clascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
431 $ info )
432 END IF
433 IF( ibscl.EQ.1 ) THEN
434 CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb,
435 $ info )
436 ELSE IF( ibscl.EQ.2 ) THEN
437 CALL clascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb,
438 $ info )
439 END IF
440*
441 100 CONTINUE
442*
443 RETURN
444*
445* End of CGELSX
446*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cgeqpf
subroutine cgeqpf(m, n, a, lda, jpvt, tau, work, rwork, info)
CGEQPF
Definition cgeqpf.f:146
clatzm
subroutine clatzm(side, m, n, v, incv, tau, c1, c2, ldc, work)
CLATZM
Definition clatzm.f:151
ctzrqf
subroutine ctzrqf(m, n, a, lda, tau, info)
CTZRQF
Definition ctzrqf.f:136
claic1
subroutine claic1(job, j, x, sest, w, gamma, sestpr, s, c)
CLAIC1 applies one step of incremental condition estimation.
Definition claic1.f:133
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
clange
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:113
clascl
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:142
claset
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:104
ctrsm
subroutine ctrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRSM
Definition ctrsm.f:180
cunm2r
subroutine cunm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
CUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition cunm2r.f:157
Here is the call graph for this function: