LAPACK  3.8.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.
Date
December 2016
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 212 of file cgelsy.f.

212 *
213 * -- LAPACK driver routine (version 3.7.0) --
214 * -- LAPACK is a software package provided by Univ. of Tennessee, --
215 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216 * December 2016
217 *
218 * .. Scalar Arguments ..
219  INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
220  REAL rcond
221 * ..
222 * .. Array Arguments ..
223  INTEGER jpvt( * )
224  REAL rwork( * )
225  COMPLEX a( lda, * ), b( ldb, * ), work( * )
226 * ..
227 *
228 * =====================================================================
229 *
230 * .. Parameters ..
231  INTEGER imax, imin
232  parameter( imax = 1, imin = 2 )
233  REAL zero, one
234  parameter( zero = 0.0e+0, one = 1.0e+0 )
235  COMPLEX czero, cone
236  parameter( czero = ( 0.0e+0, 0.0e+0 ),
237  $ cone = ( 1.0e+0, 0.0e+0 ) )
238 * ..
239 * .. Local Scalars ..
240  LOGICAL lquery
241  INTEGER i, iascl, ibscl, ismax, ismin, j, lwkopt, mn,
242  $ nb, nb1, nb2, nb3, nb4
243  REAL anrm, bignum, bnrm, smax, smaxpr, smin, sminpr,
244  $ smlnum, wsize
245  COMPLEX c1, c2, s1, s2
246 * ..
247 * .. External Subroutines ..
248  EXTERNAL ccopy, cgeqp3, claic1, clascl, claset, ctrsm,
250 * ..
251 * .. External Functions ..
252  INTEGER ilaenv
253  REAL clange, slamch
254  EXTERNAL clange, ilaenv, slamch
255 * ..
256 * .. Intrinsic Functions ..
257  INTRINSIC abs, max, min, REAL, cmplx
258 * ..
259 * .. Executable Statements ..
260 *
261  mn = min( m, n )
262  ismin = mn + 1
263  ismax = 2*mn + 1
264 *
265 * Test the input arguments.
266 *
267  info = 0
268  nb1 = ilaenv( 1, 'CGEQRF', ' ', m, n, -1, -1 )
269  nb2 = ilaenv( 1, 'CGERQF', ' ', m, n, -1, -1 )
270  nb3 = ilaenv( 1, 'CUNMQR', ' ', m, n, nrhs, -1 )
271  nb4 = ilaenv( 1, 'CUNMRQ', ' ', m, n, nrhs, -1 )
272  nb = max( nb1, nb2, nb3, nb4 )
273  lwkopt = max( 1, mn+2*n+nb*(n+1), 2*mn+nb*nrhs )
274  work( 1 ) = cmplx( lwkopt )
275  lquery = ( lwork.EQ.-1 )
276  IF( m.LT.0 ) THEN
277  info = -1
278  ELSE IF( n.LT.0 ) THEN
279  info = -2
280  ELSE IF( nrhs.LT.0 ) THEN
281  info = -3
282  ELSE IF( lda.LT.max( 1, m ) ) THEN
283  info = -5
284  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
285  info = -7
286  ELSE IF( lwork.LT.( mn+max( 2*mn, n+1, mn+nrhs ) ) .AND.
287  $ .NOT.lquery ) THEN
288  info = -12
289  END IF
290 *
291  IF( info.NE.0 ) THEN
292  CALL xerbla( 'CGELSY', -info )
293  RETURN
294  ELSE IF( lquery ) THEN
295  RETURN
296  END IF
297 *
298 * Quick return if possible
299 *
300  IF( min( m, n, nrhs ).EQ.0 ) THEN
301  rank = 0
302  RETURN
303  END IF
304 *
305 * Get machine parameters
306 *
307  smlnum = slamch( 'S' ) / slamch( 'P' )
308  bignum = one / smlnum
309  CALL slabad( smlnum, bignum )
310 *
311 * Scale A, B if max entries outside range [SMLNUM,BIGNUM]
312 *
313  anrm = clange( 'M', m, n, a, lda, rwork )
314  iascl = 0
315  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
316 *
317 * Scale matrix norm up to SMLNUM
318 *
319  CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
320  iascl = 1
321  ELSE IF( anrm.GT.bignum ) THEN
322 *
323 * Scale matrix norm down to BIGNUM
324 *
325  CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
326  iascl = 2
327  ELSE IF( anrm.EQ.zero ) THEN
328 *
329 * Matrix all zero. Return zero solution.
330 *
331  CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
332  rank = 0
333  GO TO 70
334  END IF
335 *
336  bnrm = clange( 'M', m, nrhs, b, ldb, rwork )
337  ibscl = 0
338  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
339 *
340 * Scale matrix norm up to SMLNUM
341 *
342  CALL clascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
343  ibscl = 1
344  ELSE IF( bnrm.GT.bignum ) THEN
345 *
346 * Scale matrix norm down to BIGNUM
347 *
348  CALL clascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
349  ibscl = 2
350  END IF
351 *
352 * Compute QR factorization with column pivoting of A:
353 * A * P = Q * R
354 *
355  CALL cgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
356  $ lwork-mn, rwork, info )
357  wsize = mn + REAL( WORK( MN+1 ) )
358 *
359 * complex workspace: MN+NB*(N+1). real workspace 2*N.
360 * Details of Householder rotations stored in WORK(1:MN).
361 *
362 * Determine RANK using incremental condition estimation
363 *
364  work( ismin ) = cone
365  work( ismax ) = cone
366  smax = abs( a( 1, 1 ) )
367  smin = smax
368  IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
369  rank = 0
370  CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
371  GO TO 70
372  ELSE
373  rank = 1
374  END IF
375 *
376  10 CONTINUE
377  IF( rank.LT.mn ) THEN
378  i = rank + 1
379  CALL claic1( imin, rank, work( ismin ), smin, a( 1, i ),
380  $ a( i, i ), sminpr, s1, c1 )
381  CALL claic1( imax, rank, work( ismax ), smax, a( 1, i ),
382  $ a( i, i ), smaxpr, s2, c2 )
383 *
384  IF( smaxpr*rcond.LE.sminpr ) THEN
385  DO 20 i = 1, rank
386  work( ismin+i-1 ) = s1*work( ismin+i-1 )
387  work( ismax+i-1 ) = s2*work( ismax+i-1 )
388  20 CONTINUE
389  work( ismin+rank ) = c1
390  work( ismax+rank ) = c2
391  smin = sminpr
392  smax = smaxpr
393  rank = rank + 1
394  GO TO 10
395  END IF
396  END IF
397 *
398 * complex workspace: 3*MN.
399 *
400 * Logically partition R = [ R11 R12 ]
401 * [ 0 R22 ]
402 * where R11 = R(1:RANK,1:RANK)
403 *
404 * [R11,R12] = [ T11, 0 ] * Y
405 *
406  IF( rank.LT.n )
407  $ CALL ctzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
408  $ lwork-2*mn, info )
409 *
410 * complex workspace: 2*MN.
411 * Details of Householder rotations stored in WORK(MN+1:2*MN)
412 *
413 * B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
414 *
415  CALL cunmqr( 'Left', 'Conjugate transpose', m, nrhs, mn, a, lda,
416  $ work( 1 ), b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
417  wsize = max( wsize, 2*mn+REAL( WORK( 2*MN+1 ) ) )
418 *
419 * complex workspace: 2*MN+NB*NRHS.
420 *
421 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
422 *
423  CALL ctrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
424  $ nrhs, cone, a, lda, b, ldb )
425 *
426  DO 40 j = 1, nrhs
427  DO 30 i = rank + 1, n
428  b( i, j ) = czero
429  30 CONTINUE
430  40 CONTINUE
431 *
432 * B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
433 *
434  IF( rank.LT.n ) THEN
435  CALL cunmrz( 'Left', 'Conjugate transpose', n, nrhs, rank,
436  $ n-rank, a, lda, work( mn+1 ), b, ldb,
437  $ work( 2*mn+1 ), lwork-2*mn, info )
438  END IF
439 *
440 * complex workspace: 2*MN+NRHS.
441 *
442 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
443 *
444  DO 60 j = 1, nrhs
445  DO 50 i = 1, n
446  work( jpvt( i ) ) = b( i, j )
447  50 CONTINUE
448  CALL ccopy( n, work( 1 ), 1, b( 1, j ), 1 )
449  60 CONTINUE
450 *
451 * complex workspace: N.
452 *
453 * Undo scaling
454 *
455  IF( iascl.EQ.1 ) THEN
456  CALL clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
457  CALL clascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
458  $ info )
459  ELSE IF( iascl.EQ.2 ) THEN
460  CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
461  CALL clascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
462  $ info )
463  END IF
464  IF( ibscl.EQ.1 ) THEN
465  CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
466  ELSE IF( ibscl.EQ.2 ) THEN
467  CALL clascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
468  END IF
469 *
470  70 CONTINUE
471  work( 1 ) = cmplx( lwkopt )
472 *
473  RETURN
474 *
475 * End of CGELSY
476 *
subroutine cunmrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMRZ
Definition: cunmrz.f:189
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:108
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:145
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
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:117
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:182
subroutine claic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
CLAIC1 applies one step of incremental condition estimation.
Definition: claic1.f:137
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine cgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
CGEQP3
Definition: cgeqp3.f:161
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83
subroutine ctzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CTZRZF
Definition: ctzrzf.f:153
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:170
Here is the call graph for this function:
Here is the caller graph for this function: