LAPACK 3.11.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.
[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: