SRC/dgelsy.f

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