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

◆ 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.
          If M = 0 or N = 0, B is not referenced.
[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.
          If NRHS = 0, RANK = 0 on output.
[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 210 of file cgelsy.f.

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