LAPACK 3.11.0
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.
[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.
[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.

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 DOUBLE PRECISION RCOND
218* ..
219* .. Array Arguments ..
220 INTEGER JPVT( * )
221 DOUBLE PRECISION RWORK( * )
222 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
223* ..
224*
225* =====================================================================
226*
227* .. Parameters ..
228 INTEGER IMAX, IMIN
229 parameter( imax = 1, imin = 2 )
230 DOUBLE PRECISION ZERO, ONE
231 parameter( zero = 0.0d+0, one = 1.0d+0 )
232 COMPLEX*16 CZERO, CONE
233 parameter( czero = ( 0.0d+0, 0.0d+0 ),
234 $ cone = ( 1.0d+0, 0.0d+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 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
241 $ SMLNUM, WSIZE
242 COMPLEX*16 C1, C2, S1, S2
243* ..
244* .. External Subroutines ..
245 EXTERNAL dlabad, xerbla, zcopy, zgeqp3, zlaic1, zlascl,
247* ..
248* .. External Functions ..
249 INTEGER ILAENV
250 DOUBLE PRECISION DLAMCH, ZLANGE
251 EXTERNAL ilaenv, dlamch, zlange
252* ..
253* .. Intrinsic Functions ..
254 INTRINSIC abs, dble, dcmplx, max, min
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, 'ZGEQRF', ' ', m, n, -1, -1 )
266 nb2 = ilaenv( 1, 'ZGERQF', ' ', m, n, -1, -1 )
267 nb3 = ilaenv( 1, 'ZUNMQR', ' ', m, n, nrhs, -1 )
268 nb4 = ilaenv( 1, 'ZUNMRQ', ' ', 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 ) = dcmplx( 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. .NOT.
284 $ lquery ) THEN
285 info = -12
286 END IF
287*
288 IF( info.NE.0 ) THEN
289 CALL xerbla( 'ZGELSY', -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 = dlamch( 'S' ) / dlamch( 'P' )
305 bignum = one / smlnum
306 CALL dlabad( smlnum, bignum )
307*
308* Scale A, B if max entries outside range [SMLNUM,BIGNUM]
309*
310 anrm = zlange( '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 zlascl( '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 zlascl( '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 zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
329 rank = 0
330 GO TO 70
331 END IF
332*
333 bnrm = zlange( '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 zlascl( '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 zlascl( '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 zgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
353 $ lwork-mn, rwork, info )
354 wsize = mn + dble( 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 zlaset( '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 zlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
377 $ a( i, i ), sminpr, s1, c1 )
378 CALL zlaic1( 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 ztzrzf( 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 zunmqr( '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+dble( 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 ztrsm( '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 zunmrz( '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 zcopy( 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 zlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
454 CALL zlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
455 $ info )
456 ELSE IF( iascl.EQ.2 ) THEN
457 CALL zlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
458 CALL zlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
459 $ info )
460 END IF
461 IF( ibscl.EQ.1 ) THEN
462 CALL zlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
463 ELSE IF( ibscl.EQ.2 ) THEN
464 CALL zlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
465 END IF
466*
467 70 CONTINUE
468 work( 1 ) = dcmplx( lwkopt )
469*
470 RETURN
471*
472* End of ZGELSY
473*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.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 zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:180
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:115
subroutine zgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
ZGEQP3
Definition: zgeqp3.f:159
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:143
subroutine zlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
ZLAIC1 applies one step of incremental condition estimation.
Definition: zlaic1.f:135
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:106
subroutine zunmrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMRZ
Definition: zunmrz.f:187
subroutine ztzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZTZRZF
Definition: ztzrzf.f:151
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:167
Here is the call graph for this function:
Here is the caller graph for this function: