LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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: