LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cgelsy()

subroutine cgelsy ( 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,
integer  LWORK,
real, dimension( * )  RWORK,
integer  INFO 
)

CGELSY solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
 CGELSY 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.

 This routine is basically identical to the original xGELSX except
 three differences:
   o The permutation of matrix B (the right hand side) is faster and
     more simple.
   o The call to the subroutine xGEQPF has been substituted by the
     the call to the subroutine xGEQP3. This subroutine is a Blas-3
     version of the QR factorization with column pivoting.
   o Matrix B (the right hand side) is updated with Blas-3.
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.
[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 permuted
          to the front of AP, otherwise column i is a free column.
          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 (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          The unblocked strategy requires that:
            LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
          where MN = min(M,N).
          The block algorithm requires that:
            LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
          where NB is an upper bound on the blocksize returned
          by ILAENV for the routines CGEQP3, CTZRZF, CTZRQF, CUNMQR,
          and CUNMRZ.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[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.
Contributors:
A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain

Definition at line 208 of file cgelsy.f.

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 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
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
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 cgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
CGEQP3
Definition: cgeqp3.f:159
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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: