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

◆ zgelsy()

subroutine zgelsy ( integer m,
integer n,
integer nrhs,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldb, * ) b,
integer ldb,
integer, dimension( * ) jpvt,
double precision rcond,
integer rank,
complex*16, dimension( * ) work,
integer lwork,
double precision, dimension( * ) rwork,
integer info )

ZGELSY solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
!>
!> ZGELSY 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*16 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*16 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 DOUBLE PRECISION
!>          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*16 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 ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR,
!>          and ZUNMRZ.
!>
!>          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 DOUBLE PRECISION 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 zgelsy.f.

211*
212* -- LAPACK driver routine --
213* -- LAPACK is a software package provided by Univ. of Tennessee, --
214* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215*
216* .. Scalar Arguments ..
217 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
218 DOUBLE PRECISION RCOND
219* ..
220* .. Array Arguments ..
221 INTEGER JPVT( * )
222 DOUBLE PRECISION RWORK( * )
223 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
224* ..
225*
226* =====================================================================
227*
228* .. Parameters ..
229 INTEGER IMAX, IMIN
230 parameter( imax = 1, imin = 2 )
231 DOUBLE PRECISION ZERO, ONE
232 parameter( zero = 0.0d+0, one = 1.0d+0 )
233 COMPLEX*16 CZERO, CONE
234 parameter( czero = ( 0.0d+0, 0.0d+0 ),
235 $ cone = ( 1.0d+0, 0.0d+0 ) )
236* ..
237* .. Local Scalars ..
238 LOGICAL LQUERY
239 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN,
240 $ NB, NB1, NB2, NB3, NB4
241 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
242 $ SMLNUM, WSIZE
243 COMPLEX*16 C1, C2, S1, S2
244* ..
245* .. External Subroutines ..
246 EXTERNAL xerbla, zcopy, zgeqp3, zlaic1,
247 $ zlascl,
248 $ zlaset, ztrsm, ztzrzf, zunmqr, zunmrz
249* ..
250* .. External Functions ..
251 INTEGER ILAENV
252 DOUBLE PRECISION DLAMCH, ZLANGE
253 EXTERNAL ilaenv, dlamch, zlange
254* ..
255* .. Intrinsic Functions ..
256 INTRINSIC abs, dble, dcmplx, max, min
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, 'ZGEQRF', ' ', m, n, -1, -1 )
268 nb2 = ilaenv( 1, 'ZGERQF', ' ', m, n, -1, -1 )
269 nb3 = ilaenv( 1, 'ZUNMQR', ' ', m, n, nrhs, -1 )
270 nb4 = ilaenv( 1, 'ZUNMRQ', ' ', 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 ) = dcmplx( 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. .NOT.
286 $ lquery ) THEN
287 info = -12
288 END IF
289*
290 IF( info.NE.0 ) THEN
291 CALL xerbla( 'ZGELSY', -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 = dlamch( 'S' ) / dlamch( 'P' )
307 bignum = one / smlnum
308*
309* Scale A, B if max entries outside range [SMLNUM,BIGNUM]
310*
311 anrm = zlange( '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 zlascl( '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 zlascl( '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 zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
330 rank = 0
331 GO TO 70
332 END IF
333*
334 bnrm = zlange( '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 zlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb,
341 $ info )
342 ibscl = 1
343 ELSE IF( bnrm.GT.bignum ) THEN
344*
345* Scale matrix norm down to BIGNUM
346*
347 CALL zlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb,
348 $ 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 zgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
356 $ lwork-mn, rwork, info )
357 wsize = mn + dble( 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 zlaset( '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 zlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
380 $ a( i, i ), sminpr, s1, c1 )
381 CALL zlaic1( 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 ztzrzf( 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 zunmqr( 'Left', 'Conjugate transpose', m, nrhs, mn, a,
416 $ lda,
417 $ work( 1 ), b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
418 wsize = max( wsize, 2*mn+dble( work( 2*mn+1 ) ) )
419*
420* complex workspace: 2*MN+NB*NRHS.
421*
422* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
423*
424 CALL ztrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
425 $ nrhs, cone, a, lda, b, ldb )
426*
427 DO 40 j = 1, nrhs
428 DO 30 i = rank + 1, n
429 b( i, j ) = czero
430 30 CONTINUE
431 40 CONTINUE
432*
433* B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
434*
435 IF( rank.LT.n ) THEN
436 CALL zunmrz( 'Left', 'Conjugate transpose', n, nrhs, rank,
437 $ n-rank, a, lda, work( mn+1 ), b, ldb,
438 $ work( 2*mn+1 ), lwork-2*mn, info )
439 END IF
440*
441* complex workspace: 2*MN+NRHS.
442*
443* B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
444*
445 DO 60 j = 1, nrhs
446 DO 50 i = 1, n
447 work( jpvt( i ) ) = b( i, j )
448 50 CONTINUE
449 CALL zcopy( n, work( 1 ), 1, b( 1, j ), 1 )
450 60 CONTINUE
451*
452* complex workspace: N.
453*
454* Undo scaling
455*
456 IF( iascl.EQ.1 ) THEN
457 CALL zlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb,
458 $ info )
459 CALL zlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
460 $ info )
461 ELSE IF( iascl.EQ.2 ) THEN
462 CALL zlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb,
463 $ info )
464 CALL zlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
465 $ info )
466 END IF
467 IF( ibscl.EQ.1 ) THEN
468 CALL zlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb,
469 $ info )
470 ELSE IF( ibscl.EQ.2 ) THEN
471 CALL zlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb,
472 $ info )
473 END IF
474*
475 70 CONTINUE
476 work( 1 ) = dcmplx( lwkopt )
477*
478 RETURN
479*
480* End of ZGELSY
481*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zgeqp3
subroutine zgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info)
ZGEQP3
Definition zgeqp3.f:157
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
zlaic1
subroutine zlaic1(job, j, x, sest, w, gamma, sestpr, s, c)
ZLAIC1 applies one step of incremental condition estimation.
Definition zlaic1.f:133
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlange
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition zlange.f:113
zlascl
subroutine zlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition zlascl.f:142
zlaset
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:104
ztrsm
subroutine ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM
Definition ztrsm.f:180
ztzrzf
subroutine ztzrzf(m, n, a, lda, tau, work, lwork, info)
ZTZRZF
Definition ztzrzf.f:149
zunmqr
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR
Definition zunmqr.f:165
zunmrz
subroutine zunmrz(side, trans, m, n, k, l, a, lda, tau, c, ldc, work, lwork, info)
ZUNMRZ
Definition zunmrz.f:186
Here is the call graph for this function:
Here is the caller graph for this function: