 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ 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

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```

Definition at line 182 of file cgelsx.f.

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  REAL RCOND
192 * ..
193 * .. Array Arguments ..
194  INTEGER JPVT( * )
195  REAL RWORK( * )
196  COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
197 * ..
198 *
199 * =====================================================================
200 *
201 * .. Parameters ..
202  INTEGER IMAX, IMIN
203  parameter( imax = 1, imin = 2 )
204  REAL ZERO, ONE, DONE, NTDONE
205  parameter( zero = 0.0e+0, one = 1.0e+0, done = zero,
206  \$ ntdone = one )
207  COMPLEX CZERO, CONE
208  parameter( czero = ( 0.0e+0, 0.0e+0 ),
209  \$ cone = ( 1.0e+0, 0.0e+0 ) )
210 * ..
211 * .. Local Scalars ..
212  INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
213  REAL ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
214  \$ SMLNUM
215  COMPLEX C1, C2, S1, S2, T1, T2
216 * ..
217 * .. External Subroutines ..
218  EXTERNAL cgeqpf, claic1, clascl, claset, clatzm, ctrsm,
220 * ..
221 * .. External Functions ..
222  REAL CLANGE, SLAMCH
223  EXTERNAL clange, slamch
224 * ..
225 * .. Intrinsic Functions ..
226  INTRINSIC abs, conjg, 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( 'CGELSX', -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 = slamch( 'S' ) / slamch( 'P' )
264  bignum = one / smlnum
265  CALL slabad( smlnum, bignum )
266 *
267 * Scale A, B if max elements outside range [SMLNUM,BIGNUM]
268 *
269  anrm = clange( '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 clascl( '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 clascl( '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 claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
288  rank = 0
289  GO TO 100
290  END IF
291 *
292  bnrm = clange( '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 clascl( '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 clascl( '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 cgeqpf( 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 claset( '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 claic1( imin, rank, work( ismin ), smin, a( 1, i ),
335  \$ a( i, i ), sminpr, s1, c1 )
336  CALL claic1( 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 ctzrqf( 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 cunm2r( '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 ctrsm( '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 clatzm( 'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
387  \$ conjg( 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 clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, 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, info )
429  CALL clascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
430  \$ info )
431  END IF
432  IF( ibscl.EQ.1 ) THEN
433  CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
434  ELSE IF( ibscl.EQ.2 ) THEN
435  CALL clascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
436  END IF
437 *
438  100 CONTINUE
439 *
440  RETURN
441 *
442 * End of CGELSX
443 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:180
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:115
subroutine cgeqpf(M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO)
CGEQPF
Definition: cgeqpf.f:148
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 ctzrqf(M, N, A, LDA, TAU, INFO)
CTZRQF
Definition: ctzrqf.f:138
subroutine clatzm(SIDE, M, N, V, INCV, TAU, C1, C2, LDC, WORK)
CLATZM
Definition: clatzm.f:152
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:159
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function: