SRC/zgelsy.f

   1 *> \brief <b> ZGELSY solves overdetermined or underdetermined systems for GE matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZGELSY + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelsy.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelsy.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelsy.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
  22 *                          WORK, LWORK, RWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
  26 *       DOUBLE PRECISION   RCOND
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            JPVT( * )
  30 *       DOUBLE PRECISION   RWORK( * )
  31 *       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> ZGELSY computes the minimum-norm solution to a complex linear least
  41 *> squares problem:
  42 *>     minimize || A * X - B ||
  43 *> using a complete orthogonal factorization of A.  A is an M-by-N
  44 *> matrix which may be rank-deficient.
  45 *>
  46 *> Several right hand side vectors b and solution vectors x can be
  47 *> handled in a single call; they are stored as the columns of the
  48 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
  49 *> matrix X.
  50 *>
  51 *> The routine first computes a QR factorization with column pivoting:
  52 *>     A * P = Q * [ R11 R12 ]
  53 *>                 [  0  R22 ]
  54 *> with R11 defined as the largest leading submatrix whose estimated
  55 *> condition number is less than 1/RCOND.  The order of R11, RANK,
  56 *> is the effective rank of A.
  57 *>
  58 *> Then, R22 is considered to be negligible, and R12 is annihilated
  59 *> by unitary transformations from the right, arriving at the
  60 *> complete orthogonal factorization:
  61 *>    A * P = Q * [ T11 0 ] * Z
  62 *>                [  0  0 ]
  63 *> The minimum-norm solution is then
  64 *>    X = P * Z**H [ inv(T11)*Q1**H*B ]
  65 *>                 [        0         ]
  66 *> where Q1 consists of the first RANK columns of Q.
  67 *>
  68 *> This routine is basically identical to the original xGELSX except
  69 *> three differences:
  70 *>   o The permutation of matrix B (the right hand side) is faster and
  71 *>     more simple.
  72 *>   o The call to the subroutine xGEQPF has been substituted by the
  73 *>     the call to the subroutine xGEQP3. This subroutine is a Blas-3
  74 *>     version of the QR factorization with column pivoting.
  75 *>   o Matrix B (the right hand side) is updated with Blas-3.
  76 *> \endverbatim
  77 *
  78 *  Arguments:
  79 *  ==========
  80 *
  81 *> \param[in] M
  82 *> \verbatim
  83 *>          M is INTEGER
  84 *>          The number of rows of the matrix A.  M >= 0.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] N
  88 *> \verbatim
  89 *>          N is INTEGER
  90 *>          The number of columns of the matrix A.  N >= 0.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] NRHS
  94 *> \verbatim
  95 *>          NRHS is INTEGER
  96 *>          The number of right hand sides, i.e., the number of
  97 *>          columns of matrices B and X. NRHS >= 0.
  98 *> \endverbatim
  99 *>
 100 *> \param[in,out] A
 101 *> \verbatim
 102 *>          A is COMPLEX*16 array, dimension (LDA,N)
 103 *>          On entry, the M-by-N matrix A.
 104 *>          On exit, A has been overwritten by details of its
 105 *>          complete orthogonal factorization.
 106 *> \endverbatim
 107 *>
 108 *> \param[in] LDA
 109 *> \verbatim
 110 *>          LDA is INTEGER
 111 *>          The leading dimension of the array A.  LDA >= max(1,M).
 112 *> \endverbatim
 113 *>
 114 *> \param[in,out] B
 115 *> \verbatim
 116 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
 117 *>          On entry, the M-by-NRHS right hand side matrix B.
 118 *>          On exit, the N-by-NRHS solution matrix X.
 119 *> \endverbatim
 120 *>
 121 *> \param[in] LDB
 122 *> \verbatim
 123 *>          LDB is INTEGER
 124 *>          The leading dimension of the array B. LDB >= max(1,M,N).
 125 *> \endverbatim
 126 *>
 127 *> \param[in,out] JPVT
 128 *> \verbatim
 129 *>          JPVT is INTEGER array, dimension (N)
 130 *>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
 131 *>          to the front of AP, otherwise column i is a free column.
 132 *>          On exit, if JPVT(i) = k, then the i-th column of A*P
 133 *>          was the k-th column of A.
 134 *> \endverbatim
 135 *>
 136 *> \param[in] RCOND
 137 *> \verbatim
 138 *>          RCOND is DOUBLE PRECISION
 139 *>          RCOND is used to determine the effective rank of A, which
 140 *>          is defined as the order of the largest leading triangular
 141 *>          submatrix R11 in the QR factorization with pivoting of A,
 142 *>          whose estimated condition number < 1/RCOND.
 143 *> \endverbatim
 144 *>
 145 *> \param[out] RANK
 146 *> \verbatim
 147 *>          RANK is INTEGER
 148 *>          The effective rank of A, i.e., the order of the submatrix
 149 *>          R11.  This is the same as the order of the submatrix T11
 150 *>          in the complete orthogonal factorization of A.
 151 *> \endverbatim
 152 *>
 153 *> \param[out] WORK
 154 *> \verbatim
 155 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 156 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 157 *> \endverbatim
 158 *>
 159 *> \param[in] LWORK
 160 *> \verbatim
 161 *>          LWORK is INTEGER
 162 *>          The dimension of the array WORK.
 163 *>          The unblocked strategy requires that:
 164 *>            LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
 165 *>          where MN = min(M,N).
 166 *>          The block algorithm requires that:
 167 *>            LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
 168 *>          where NB is an upper bound on the blocksize returned
 169 *>          by ILAENV for the routines ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR,
 170 *>          and ZUNMRZ.
 171 *>
 172 *>          If LWORK = -1, then a workspace query is assumed; the routine
 173 *>          only calculates the optimal size of the WORK array, returns
 174 *>          this value as the first entry of the WORK array, and no error
 175 *>          message related to LWORK is issued by XERBLA.
 176 *> \endverbatim
 177 *>
 178 *> \param[out] RWORK
 179 *> \verbatim
 180 *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
 181 *> \endverbatim
 182 *>
 183 *> \param[out] INFO
 184 *> \verbatim
 185 *>          INFO is INTEGER
 186 *>          = 0: successful exit
 187 *>          < 0: if INFO = -i, the i-th argument had an illegal value
 188 *> \endverbatim
 189 *
 190 *  Authors:
 191 *  ========
 192 *
 193 *> \author Univ. of Tennessee
 194 *> \author Univ. of California Berkeley
 195 *> \author Univ. of Colorado Denver
 196 *> \author NAG Ltd.
 197 *
 198 *> \date November 2011
 199 *
 200 *> \ingroup complex16GEsolve
 201 *
 202 *> \par Contributors:
 203 *  ==================
 204 *>
 205 *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n
 206 *>    E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
 207 *>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
 208 *>
 209 *  =====================================================================
 210       SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
 211      $                   WORK, LWORK, RWORK, INFO )
 212 *
 213 *  -- LAPACK driver routine (version 3.4.0) --
 214 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 215 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 216 *     November 2011
 217 *
 218 *     .. Scalar Arguments ..
 219       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
 220       DOUBLE PRECISION   RCOND
 221 *     ..
 222 *     .. Array Arguments ..
 223       INTEGER            JPVT( * )
 224       DOUBLE PRECISION   RWORK( * )
 225       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
 226 *     ..
 227 *
 228 *  =====================================================================
 229 *
 230 *     .. Parameters ..
 231       INTEGER            IMAX, IMIN
 232       PARAMETER          ( IMAX = 1, IMIN = 2 )
 233       DOUBLE PRECISION   ZERO, ONE
 234       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 235       COMPLEX*16         CZERO, CONE
 236       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
 237      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
 238 *     ..
 239 *     .. Local Scalars ..
 240       LOGICAL            LQUERY
 241       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN,
 242      $                   NB, NB1, NB2, NB3, NB4
 243       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
 244      $                   SMLNUM, WSIZE
 245       COMPLEX*16         C1, C2, S1, S2
 246 *     ..
 247 *     .. External Subroutines ..
 248       EXTERNAL           DLABAD, XERBLA, ZCOPY, ZGEQP3, ZLAIC1, ZLASCL,
 249      $                   ZLASET, ZTRSM, ZTZRZF, ZUNMQR, ZUNMRZ
 250 *     ..
 251 *     .. External Functions ..
 252       INTEGER            ILAENV
 253       DOUBLE PRECISION   DLAMCH, ZLANGE
 254       EXTERNAL           ILAENV, DLAMCH, ZLANGE
 255 *     ..
 256 *     .. Intrinsic Functions ..
 257       INTRINSIC          ABS, DBLE, DCMPLX, MAX, MIN
 258 *     ..
 259 *     .. Executable Statements ..
 260 *
 261       MN = MIN( M, N )
 262       ISMIN = MN + 1
 263       ISMAX = 2*MN + 1
 264 *
 265 *     Test the input arguments.
 266 *
 267       INFO = 0
 268       NB1 = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
 269       NB2 = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
 270       NB3 = ILAENV( 1, 'ZUNMQR', ' ', M, N, NRHS, -1 )
 271       NB4 = ILAENV( 1, 'ZUNMRQ', ' ', M, N, NRHS, -1 )
 272       NB = MAX( NB1, NB2, NB3, NB4 )
 273       LWKOPT = MAX( 1, MN+2*N+NB*( N+1 ), 2*MN+NB*NRHS )
 274       WORK( 1 ) = DCMPLX( LWKOPT )
 275       LQUERY = ( LWORK.EQ.-1 )
 276       IF( M.LT.0 ) THEN
 277          INFO = -1
 278       ELSE IF( N.LT.0 ) THEN
 279          INFO = -2
 280       ELSE IF( NRHS.LT.0 ) THEN
 281          INFO = -3
 282       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 283          INFO = -5
 284       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
 285          INFO = -7
 286       ELSE IF( LWORK.LT.( MN+MAX( 2*MN, N+1, MN+NRHS ) ) .AND. .NOT.
 287      $         LQUERY ) THEN
 288          INFO = -12
 289       END IF
 290 *
 291       IF( INFO.NE.0 ) THEN
 292          CALL XERBLA( 'ZGELSY', -INFO )
 293          RETURN
 294       ELSE IF( LQUERY ) THEN
 295          RETURN
 296       END IF
 297 *
 298 *     Quick return if possible
 299 *
 300       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
 301          RANK = 0
 302          RETURN
 303       END IF
 304 *
 305 *     Get machine parameters
 306 *
 307       SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
 308       BIGNUM = ONE / SMLNUM
 309       CALL DLABAD( SMLNUM, BIGNUM )
 310 *
 311 *     Scale A, B if max entries outside range [SMLNUM,BIGNUM]
 312 *
 313       ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
 314       IASCL = 0
 315       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 316 *
 317 *        Scale matrix norm up to SMLNUM
 318 *
 319          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
 320          IASCL = 1
 321       ELSE IF( ANRM.GT.BIGNUM ) THEN
 322 *
 323 *        Scale matrix norm down to BIGNUM
 324 *
 325          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
 326          IASCL = 2
 327       ELSE IF( ANRM.EQ.ZERO ) THEN
 328 *
 329 *        Matrix all zero. Return zero solution.
 330 *
 331          CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
 332          RANK = 0
 333          GO TO 70
 334       END IF
 335 *
 336       BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
 337       IBSCL = 0
 338       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 339 *
 340 *        Scale matrix norm up to SMLNUM
 341 *
 342          CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
 343          IBSCL = 1
 344       ELSE IF( BNRM.GT.BIGNUM ) THEN
 345 *
 346 *        Scale matrix norm down to BIGNUM
 347 *
 348          CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, 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, LDA,
 416      $             WORK( 1 ), B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO )
 417       WSIZE = MAX( WSIZE, 2*MN+DBLE( WORK( 2*MN+1 ) ) )
 418 *
 419 *     complex workspace: 2*MN+NB*NRHS.
 420 *
 421 *     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
 422 *
 423       CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
 424      $            NRHS, CONE, A, LDA, B, LDB )
 425 *
 426       DO 40 J = 1, NRHS
 427          DO 30 I = RANK + 1, N
 428             B( I, J ) = CZERO
 429    30    CONTINUE
 430    40 CONTINUE
 431 *
 432 *     B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
 433 *
 434       IF( RANK.LT.N ) THEN
 435          CALL ZUNMRZ( 'Left', 'Conjugate transpose', N, NRHS, RANK,
 436      $                N-RANK, A, LDA, WORK( MN+1 ), B, LDB,
 437      $                WORK( 2*MN+1 ), LWORK-2*MN, INFO )
 438       END IF
 439 *
 440 *     complex workspace: 2*MN+NRHS.
 441 *
 442 *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
 443 *
 444       DO 60 J = 1, NRHS
 445          DO 50 I = 1, N
 446             WORK( JPVT( I ) ) = B( I, J )
 447    50    CONTINUE
 448          CALL ZCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 )
 449    60 CONTINUE
 450 *
 451 *     complex workspace: N.
 452 *
 453 *     Undo scaling
 454 *
 455       IF( IASCL.EQ.1 ) THEN
 456          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
 457          CALL ZLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
 458      $                INFO )
 459       ELSE IF( IASCL.EQ.2 ) THEN
 460          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
 461          CALL ZLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
 462      $                INFO )
 463       END IF
 464       IF( IBSCL.EQ.1 ) THEN
 465          CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
 466       ELSE IF( IBSCL.EQ.2 ) THEN
 467          CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
 468       END IF
 469 *
 470    70 CONTINUE
 471       WORK( 1 ) = DCMPLX( LWKOPT )
 472 *
 473       RETURN
 474 *
 475 *     End of ZGELSY
 476 *
 477       END