SRC/cgels.f

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