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