SRC/sgetsls.f

   1 *  Definition:
   2 *  ===========
   3 *
   4 *       SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
   5 *     $                   , WORK, LWORK, INFO )
   6
   7 *
   8 *       .. Scalar Arguments ..
   9 *       CHARACTER          TRANS
  10 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
  11 *       ..
  12 *       .. Array Arguments ..
  13 *       REAL   A( LDA, * ), B( LDB, * ), WORK( * )
  14 *       ..
  15 *
  16 *
  17 *> \par Purpose:
  18 *  =============
  19 *>
  20 *> \verbatim
  21 *>
  22 *> SGETSLS solves overdetermined or underdetermined real linear systems
  23 *> involving an M-by-N matrix A, using a tall skinny QR or short wide LQ
  24 *> factorization of A.  It is assumed that A has full rank.
  25 *>
  26 *>
  27 *>
  28 *> The following options are provided:
  29 *>
  30 *> 1. If TRANS = 'N' and m >= n:  find the least squares solution of
  31 *>    an overdetermined system, i.e., solve the least squares problem
  32 *>                 minimize || B - A*X ||.
  33
  34 *> 2. If TRANS = 'N' and m < n:  find the minimum norm solution of
  35 *>    an underdetermined system A * X = B.
  36
  37 *> 3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
  38 *>    an undetermined system A**T * X = B.
  39
  40 *> 4. If TRANS = 'T' and m < n:  find the least squares solution of
  41 *>    an overdetermined system, i.e., solve the least squares problem
  42 *>                 minimize || B - A**T * X ||.
  43 *> \endverbatim
  44 *
  45 *  Arguments:
  46 *  ==========
  47 *
  48 *> \param[in] TRANS
  49 *> \verbatim
  50 *>          TRANS is CHARACTER*1
  51 *>          = 'N': the linear system involves A;
  52 *>          = 'T': the linear system involves A**T.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] M
  56 *> \verbatim
  57 *>          M is INTEGER
  58 *>          The number of rows of the matrix A.  M >= 0.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] N
  62 *> \verbatim
  63 *>          N is INTEGER
  64 *>          The number of columns of the matrix A.  N >= 0.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] NRHS
  68 *> \verbatim
  69 *>          NRHS is INTEGER
  70 *>          The number of right hand sides, i.e., the number of
  71 *>          columns of the matrices B and X. NRHS >=0.
  72 *> \endverbatim
  73 *>
  74 *> \param[in,out] A
  75 *> \verbatim
  76 *>          A is REAL array, dimension (LDA,N)
  77 *>          On entry, the M-by-N matrix A.
  78 *>          On exit,
  79 *>          A is overwritten by details of its QR or LQ
  80 *>          factorization as returned by DGETSQR.
  81 *> \endverbatim
  82 *>
  83 *> \param[in] LDA
  84 *> \verbatim
  85 *>          LDA is INTEGER
  86 *>          The leading dimension of the array A.  LDA >= max(1,M).
  87 *> \endverbatim
  88 *>
  89 *> \param[in,out] B
  90 *> \verbatim
  91 *>          B is REAL array, dimension (LDB,NRHS)
  92 *>          On entry, the matrix B of right hand side vectors, stored
  93 *>          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
  94 *>          if TRANS = 'T'.
  95 *>          On exit, if INFO = 0, B is overwritten by the solution
  96 *>          vectors, stored columnwise:
  97 *>          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
  98 *>          squares solution vectors.
  99 *>          if TRANS = 'N' and m < n, rows 1 to N of B contain the
 100 *>          minimum norm solution vectors;
 101 *>          if TRANS = 'T' and m >= n, rows 1 to M of B contain the
 102 *>          minimum norm solution vectors;
 103 *>          if TRANS = 'T' and m < n, rows 1 to M of B contain the
 104 *>          least squares solution vectors.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] LDB
 108 *> \verbatim
 109 *>          LDB is INTEGER
 110 *>          The leading dimension of the array B. LDB >= MAX(1,M,N).
 111 *> \endverbatim
 112 *>
 113 *> \param[out] WORK
 114 *> \verbatim
 115 *>          WORK is REAL array, dimension (MAX(1,LWORK))
 116 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK,
 117 *>          and WORK(2) returns the minimum LWORK.
 118 *> \endverbatim
 119 *>
 120 *> \param[in] LWORK
 121 *> \verbatim
 122 *>          LWORK is INTEGER
 123 *>          The dimension of the array WORK.
 124 *>          IF LWORK=-1, workspace query is assumed, and
 125 *>          WORK(1) returns the optimal LWORK,
 126 *>          and WORK(2) returns the minimum LWORK.
 127 *> \endverbatim
 128 *>
 129 *> \param[out] INFO
 130 *> \verbatim
 131 *>          INFO is INTEGER
 132 *>          = 0:  successful exit
 133 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 134 *>          > 0:  if INFO =  i, the i-th diagonal element of the
 135 *>                triangular factor of A is zero, so that A does not have
 136 *>                full rank; the least squares solution could not be
 137 *>                computed.
 138 *> \endverbatim
 139 *
 140 *  Authors:
 141 *  ========
 142 *
 143 *> \author Univ. of Tennessee
 144 *> \author Univ. of California Berkeley
 145 *> \author Univ. of Colorado Denver
 146 *> \author NAG Ltd.
 147 *
 148 *> \date November 2011
 149 *
 150 *> \ingroup doubleGEsolve
 151 *
 152 *  =====================================================================
 153       SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
 154      $                   , WORK, LWORK, INFO )
 155 *
 156 *  -- LAPACK driver routine (version 3.4.0) --
 157 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 158 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 159 *     November 2011
 160 *
 161 *     .. Scalar Arguments ..
 162       CHARACTER          TRANS
 163       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, MB
 164 *     ..
 165 *     .. Array Arguments ..
 166       REAL               A( LDA, * ), B( LDB, * ), WORK( * )
 167 *
 168 *     ..
 169 *
 170 *  =====================================================================
 171 *
 172 *     .. Parameters ..
 173       REAL               ZERO, ONE
 174       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
 175 *     ..
 176 *     .. Local Scalars ..
 177       LOGICAL            LQUERY, TRAN
 178       INTEGER            I, IASCL, IBSCL, J, MINMN, MAXMN, BROW, LW,
 179      $                   SCLLEN, MNK, WSIZEO, WSIZEM, LW1, LW2, INFO2,
 180      $                   NB
 181       REAL               ANRM, BIGNUM, BNRM, SMLNUM
 182 *     ..
 183 *     .. External Functions ..
 184       LOGICAL            LSAME
 185       INTEGER            ILAENV
 186       REAL               SLAMCH, SLANGE
 187       EXTERNAL           LSAME, ILAENV, SLABAD, SLAMCH, SLANGE
 188 *     ..
 189 *     .. External Subroutines ..
 190       EXTERNAL           SGEQR, SGEMQR, SLASCL, SLASET,
 191      $                   STRTRS, XERBLA, SGELQ, SGEMLQ
 192 *     ..
 193 *     .. Intrinsic Functions ..
 194       INTRINSIC          REAL, MAX, MIN
 195 *     ..
 196 *     .. Executable Statements ..
 197 *
 198 *     Test the input arguments.
 199 *
 200       INFO=0
 201       MINMN = MIN( M, N )
 202       MAXMN = MAX( M, N )
 203       MNK   = MAX(MINMN,NRHS)
 204       TRAN  = LSAME( TRANS, 'T' )
 205 *
 206       LQUERY = ( LWORK.EQ.-1 )
 207       IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
 208      $    LSAME( TRANS, 'T' ) ) ) THEN
 209          INFO = -1
 210       ELSE IF( M.LT.0 ) THEN
 211          INFO = -2
 212       ELSE IF( N.LT.0 ) THEN
 213          INFO = -3
 214       ELSE IF( NRHS.LT.0 ) THEN
 215          INFO = -4
 216       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 217          INFO = -6
 218       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
 219          INFO = -8
 220       END IF
 221 *
 222       IF( INFO.EQ.0)  THEN
 223 *
 224 *     Determine the block size and minimum LWORK
 225 *
 226        IF ( M.GE.N ) THEN
 227         CALL SGEQR( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
 228      $   INFO2)
 229         MB = INT(WORK(4))
 230         NB = INT(WORK(5))
 231         LW = INT(WORK(6))
 232         CALL SGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WORK(1),
 233      $        INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
 234         WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
 235         WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
 236        ELSE
 237         CALL SGELQ( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
 238      $   INFO2)
 239         MB = INT(WORK(4))
 240         NB = INT(WORK(5))
 241         LW = INT(WORK(6))
 242         CALL SGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WORK(1),
 243      $        INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
 244         WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
 245         WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
 246        END IF
 247 *
 248        IF((LWORK.LT.WSIZEO).AND.(.NOT.LQUERY)) THEN
 249           INFO=-10
 250        END IF
 251       END IF
 252 *
 253       IF( INFO.NE.0 ) THEN
 254         CALL XERBLA( 'SGETSLS', -INFO )
 255         WORK( 1 ) = REAL( WSIZEO )
 256         WORK( 2 ) = REAL( WSIZEM )
 257         RETURN
 258       ELSE IF (LQUERY) THEN
 259         WORK( 1 ) = REAL( WSIZEO )
 260         WORK( 2 ) = REAL( WSIZEM )
 261         RETURN
 262       END IF
 263       IF(LWORK.LT.WSIZEO) THEN
 264         LW1=INT(WORK(3))
 265         LW2=MAX(LW,INT(WORK(6)))
 266       ELSE
 267         LW1=INT(WORK(2))
 268         LW2=MAX(LW,INT(WORK(6)))
 269       END IF
 270 *
 271 *     Quick return if possible
 272 *
 273       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
 274            CALL SLASET( 'FULL', MAX( M, N ), NRHS, ZERO, ZERO,
 275      $       B, LDB )
 276            RETURN
 277       END IF
 278 *
 279 *     Get machine parameters
 280 *
 281        SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
 282        BIGNUM = ONE / SMLNUM
 283        CALL SLABAD( SMLNUM, BIGNUM )
 284 *
 285 *     Scale A, B if max element outside range [SMLNUM,BIGNUM]
 286 *
 287       ANRM = SLANGE( 'M', M, N, A, LDA, WORK )
 288       IASCL = 0
 289       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 290 *
 291 *        Scale matrix norm up to SMLNUM
 292 *
 293          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
 294          IASCL = 1
 295       ELSE IF( ANRM.GT.BIGNUM ) THEN
 296 *
 297 *        Scale matrix norm down to BIGNUM
 298 *
 299          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
 300          IASCL = 2
 301       ELSE IF( ANRM.EQ.ZERO ) THEN
 302 *
 303 *        Matrix all zero. Return zero solution.
 304 *
 305          CALL SLASET( 'F', MAXMN, NRHS, ZERO, ZERO, B, LDB )
 306          GO TO 50
 307       END IF
 308 *
 309       BROW = M
 310       IF ( TRAN ) THEN
 311         BROW = N
 312       END IF
 313       BNRM = SLANGE( 'M', BROW, NRHS, B, LDB, WORK )
 314       IBSCL = 0
 315       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 316 *
 317 *        Scale matrix norm up to SMLNUM
 318 *
 319          CALL SLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
 320      $                INFO )
 321          IBSCL = 1
 322       ELSE IF( BNRM.GT.BIGNUM ) THEN
 323 *
 324 *        Scale matrix norm down to BIGNUM
 325 *
 326          CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
 327      $                INFO )
 328          IBSCL = 2
 329       END IF
 330 *
 331       IF ( M.GE.N) THEN
 332 *
 333 *        compute QR factorization of A
 334 *
 335         CALL SGEQR( M, N, A, LDA, WORK(LW2+1), LW1
 336      $    , WORK(1), LW2, INFO )
 337         IF (.NOT.TRAN) THEN
 338 *
 339 *           Least-Squares Problem min || A * X - B ||
 340 *
 341 *           B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
 342 *
 343           CALL SGEMQR( 'L' , 'T', M, NRHS, N, A, LDA,
 344      $         WORK(LW2+1), LW1, B, LDB, WORK(1), LW2, INFO )
 345 *
 346 *           B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
 347 *
 348           CALL STRTRS( 'U', 'N', 'N', N, NRHS,
 349      $                   A, LDA, B, LDB, INFO )
 350           IF(INFO.GT.0) THEN
 351             RETURN
 352           END IF
 353           SCLLEN = N
 354         ELSE
 355 *
 356 *           Overdetermined system of equations A**T * X = B
 357 *
 358 *           B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
 359 *
 360             CALL STRTRS( 'U', 'T', 'N', N, NRHS,
 361      $                   A, LDA, B, LDB, INFO )
 362 *
 363             IF( INFO.GT.0 ) THEN
 364                RETURN
 365             END IF
 366 *
 367 *           B(N+1:M,1:NRHS) = ZERO
 368 *
 369             DO 20 J = 1, NRHS
 370                DO 10 I = N + 1, M
 371                   B( I, J ) = ZERO
 372    10          CONTINUE
 373    20       CONTINUE
 374 *
 375 *           B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
 376 *
 377             CALL SGEMQR( 'L', 'N', M, NRHS, N, A, LDA,
 378      $                   WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
 379      $                   INFO )
 380 *
 381             SCLLEN = M
 382 *
 383          END IF
 384 *
 385       ELSE
 386 *
 387 *        Compute LQ factorization of A
 388 *
 389         CALL SGELQ( M, N, A, LDA, WORK(LW2+1), LW1
 390      $    , WORK(1), LW2, INFO )
 391 *
 392 *        workspace at least M, optimally M*NB.
 393 *
 394          IF( .NOT.TRAN ) THEN
 395 *
 396 *           underdetermined system of equations A * X = B
 397 *
 398 *           B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
 399 *
 400             CALL STRTRS( 'L', 'N', 'N', M, NRHS,
 401      $                   A, LDA, B, LDB, INFO )
 402 *
 403             IF( INFO.GT.0 ) THEN
 404                RETURN
 405             END IF
 406 *
 407 *           B(M+1:N,1:NRHS) = 0
 408 *
 409             DO 40 J = 1, NRHS
 410                DO 30 I = M + 1, N
 411                   B( I, J ) = ZERO
 412    30          CONTINUE
 413    40       CONTINUE
 414 *
 415 *           B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
 416 *
 417             CALL SGEMLQ( 'L', 'T', N, NRHS, M, A, LDA,
 418      $                   WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
 419      $                   INFO )
 420 *
 421 *           workspace at least NRHS, optimally NRHS*NB
 422 *
 423             SCLLEN = N
 424 *
 425          ELSE
 426 *
 427 *           overdetermined system min || A**T * X - B ||
 428 *
 429 *           B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
 430 *
 431             CALL SGEMLQ( 'L', 'N', N, NRHS, M, A, LDA,
 432      $                   WORK( LW2+1), LW1, B, LDB, WORK( 1 ), LW2,
 433      $                   INFO )
 434 *
 435 *           workspace at least NRHS, optimally NRHS*NB
 436 *
 437 *           B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
 438 *
 439             CALL STRTRS( 'Lower', 'Transpose', 'Non-unit', M, NRHS,
 440      $                   A, LDA, B, LDB, INFO )
 441 *
 442             IF( INFO.GT.0 ) THEN
 443                RETURN
 444             END IF
 445 *
 446             SCLLEN = M
 447 *
 448          END IF
 449 *
 450       END IF
 451 *
 452 *     Undo scaling
 453 *
 454       IF( IASCL.EQ.1 ) THEN
 455         CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
 456      $                INFO )
 457       ELSE IF( IASCL.EQ.2 ) THEN
 458         CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
 459      $                INFO )
 460       END IF
 461       IF( IBSCL.EQ.1 ) THEN
 462          CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
 463      $                INFO )
 464       ELSE IF( IBSCL.EQ.2 ) THEN
 465          CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
 466      $                INFO )
 467       END IF
 468 *
 469    50 CONTINUE
 470        WORK( 1 ) = REAL( WSIZEO )
 471        WORK( 2 ) = REAL( WSIZEM )
 472       RETURN
 473 *
 474 *     End of SGETSLS
 475 *
 476       END