SRC/sgels.f

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