Annotation of rpl/lapack/lapack/zgelst.f, revision 1.1
1.1 ! bertrand 1: *> \brief <b> ZGELST 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 ZGELST + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelst.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelst.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelst.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZGELST( 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*16 A( LDA, * ), B( LDB, * ), WORK( * )
! 30: * ..
! 31: *
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *>
! 38: *> ZGELST solves overdetermined or underdetermined real linear systems
! 39: *> involving an M-by-N matrix A, or its conjugate-transpose, using a QR
! 40: *> or LQ 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 = 'C' and m >= n: find the minimum norm solution of
! 53: *> an underdetermined system A**T * X = B.
! 54: *>
! 55: *> 4. If TRANS = 'C' 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: *> = 'C': the linear system involves A**H.
! 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 COMPLEX*16 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 ZGEQRT;
! 101: *> if M < N, A is overwritten by details of its LQ
! 102: *> factorization as returned by ZGELQT.
! 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 COMPLEX*16 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 = 'C'.
! 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: *> modulus of 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 = 'C' and m >= n, rows 1 to M of B contain the
! 126: *> minimum norm solution vectors;
! 127: *> if TRANS = 'C' 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 the modulus 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 COMPLEX*16 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 complex16GEsolve
! 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 ZGELST( 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: COMPLEX*16 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: COMPLEX*16 CZERO
! 213: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
! 214: * ..
! 215: * .. Local Scalars ..
! 216: LOGICAL LQUERY, TPSD
! 217: INTEGER BROW, I, IASCL, IBSCL, J, LWOPT, MN, MNNRHS,
! 218: $ NB, NBMIN, SCLLEN
! 219: DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
! 220: * ..
! 221: * .. Local Arrays ..
! 222: DOUBLE PRECISION RWORK( 1 )
! 223: * ..
! 224: * .. External Functions ..
! 225: LOGICAL LSAME
! 226: INTEGER ILAENV
! 227: DOUBLE PRECISION DLAMCH, ZLANGE
! 228: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
! 229: * ..
! 230: * .. External Subroutines ..
! 231: EXTERNAL ZGELQT, ZGEQRT, ZGEMLQT, ZGEMQRT, DLABAD,
! 232: $ ZLASCL, ZLASET, ZTRTRS, XERBLA
! 233: * ..
! 234: * .. Intrinsic Functions ..
! 235: INTRINSIC DBLE, MAX, MIN
! 236: * ..
! 237: * .. Executable Statements ..
! 238: *
! 239: * Test the input arguments.
! 240: *
! 241: INFO = 0
! 242: MN = MIN( M, N )
! 243: LQUERY = ( LWORK.EQ.-1 )
! 244: IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'C' ) ) ) THEN
! 245: INFO = -1
! 246: ELSE IF( M.LT.0 ) THEN
! 247: INFO = -2
! 248: ELSE IF( N.LT.0 ) THEN
! 249: INFO = -3
! 250: ELSE IF( NRHS.LT.0 ) THEN
! 251: INFO = -4
! 252: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 253: INFO = -6
! 254: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
! 255: INFO = -8
! 256: ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY )
! 257: $ THEN
! 258: INFO = -10
! 259: END IF
! 260: *
! 261: * Figure out optimal block size and optimal workspace size
! 262: *
! 263: IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
! 264: *
! 265: TPSD = .TRUE.
! 266: IF( LSAME( TRANS, 'N' ) )
! 267: $ TPSD = .FALSE.
! 268: *
! 269: NB = ILAENV( 1, 'ZGELST', ' ', M, N, -1, -1 )
! 270: *
! 271: MNNRHS = MAX( MN, NRHS )
! 272: LWOPT = MAX( 1, (MN+MNNRHS)*NB )
! 273: WORK( 1 ) = DBLE( LWOPT )
! 274: *
! 275: END IF
! 276: *
! 277: IF( INFO.NE.0 ) THEN
! 278: CALL XERBLA( 'ZGELST ', -INFO )
! 279: RETURN
! 280: ELSE IF( LQUERY ) THEN
! 281: RETURN
! 282: END IF
! 283: *
! 284: * Quick return if possible
! 285: *
! 286: IF( MIN( M, N, NRHS ).EQ.0 ) THEN
! 287: CALL ZLASET( 'Full', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
! 288: WORK( 1 ) = DBLE( LWOPT )
! 289: RETURN
! 290: END IF
! 291: *
! 292: * *GEQRT and *GELQT routines cannot accept NB larger than min(M,N)
! 293: *
! 294: IF( NB.GT.MN ) NB = MN
! 295: *
! 296: * Determine the block size from the supplied LWORK
! 297: * ( at this stage we know that LWORK >= (minimum required workspace,
! 298: * but it may be less than optimal)
! 299: *
! 300: NB = MIN( NB, LWORK/( MN + MNNRHS ) )
! 301: *
! 302: * The minimum value of NB, when blocked code is used
! 303: *
! 304: NBMIN = MAX( 2, ILAENV( 2, 'ZGELST', ' ', M, N, -1, -1 ) )
! 305: *
! 306: IF( NB.LT.NBMIN ) THEN
! 307: NB = 1
! 308: END IF
! 309: *
! 310: * Get machine parameters
! 311: *
! 312: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
! 313: BIGNUM = ONE / SMLNUM
! 314: CALL DLABAD( SMLNUM, BIGNUM )
! 315: *
! 316: * Scale A, B if max element outside range [SMLNUM,BIGNUM]
! 317: *
! 318: ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
! 319: IASCL = 0
! 320: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 321: *
! 322: * Scale matrix norm up to SMLNUM
! 323: *
! 324: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
! 325: IASCL = 1
! 326: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 327: *
! 328: * Scale matrix norm down to BIGNUM
! 329: *
! 330: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
! 331: IASCL = 2
! 332: ELSE IF( ANRM.EQ.ZERO ) THEN
! 333: *
! 334: * Matrix all zero. Return zero solution.
! 335: *
! 336: CALL ZLASET( 'Full', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
! 337: WORK( 1 ) = DBLE( LWOPT )
! 338: RETURN
! 339: END IF
! 340: *
! 341: BROW = M
! 342: IF( TPSD )
! 343: $ BROW = N
! 344: BNRM = ZLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
! 345: IBSCL = 0
! 346: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
! 347: *
! 348: * Scale matrix norm up to SMLNUM
! 349: *
! 350: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
! 351: $ INFO )
! 352: IBSCL = 1
! 353: ELSE IF( BNRM.GT.BIGNUM ) THEN
! 354: *
! 355: * Scale matrix norm down to BIGNUM
! 356: *
! 357: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
! 358: $ INFO )
! 359: IBSCL = 2
! 360: END IF
! 361: *
! 362: IF( M.GE.N ) THEN
! 363: *
! 364: * M > N:
! 365: * Compute the blocked QR factorization of A,
! 366: * using the compact WY representation of Q,
! 367: * workspace at least N, optimally N*NB.
! 368: *
! 369: CALL ZGEQRT( M, N, NB, A, LDA, WORK( 1 ), NB,
! 370: $ WORK( MN*NB+1 ), INFO )
! 371: *
! 372: IF( .NOT.TPSD ) THEN
! 373: *
! 374: * M > N, A is not transposed:
! 375: * Overdetermined system of equations,
! 376: * least-squares problem, min || A * X - B ||.
! 377: *
! 378: * Compute B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS),
! 379: * using the compact WY representation of Q,
! 380: * workspace at least NRHS, optimally NRHS*NB.
! 381: *
! 382: CALL ZGEMQRT( 'Left', 'Conjugate transpose', M, NRHS, N, NB,
! 383: $ A, LDA, WORK( 1 ), NB, B, LDB,
! 384: $ WORK( MN*NB+1 ), INFO )
! 385: *
! 386: * Compute B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
! 387: *
! 388: CALL ZTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS,
! 389: $ A, LDA, B, LDB, INFO )
! 390: *
! 391: IF( INFO.GT.0 ) THEN
! 392: RETURN
! 393: END IF
! 394: *
! 395: SCLLEN = N
! 396: *
! 397: ELSE
! 398: *
! 399: * M > N, A is transposed:
! 400: * Underdetermined system of equations,
! 401: * minimum norm solution of A**T * X = B.
! 402: *
! 403: * Compute B := inv(R**T) * B in two row blocks of B.
! 404: *
! 405: * Block 1: B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
! 406: *
! 407: CALL ZTRTRS( 'Upper', 'Conjugate transpose', 'Non-unit',
! 408: $ N, NRHS, A, LDA, B, LDB, INFO )
! 409: *
! 410: IF( INFO.GT.0 ) THEN
! 411: RETURN
! 412: END IF
! 413: *
! 414: * Block 2: Zero out all rows below the N-th row in B:
! 415: * B(N+1:M,1:NRHS) = ZERO
! 416: *
! 417: DO J = 1, NRHS
! 418: DO I = N + 1, M
! 419: B( I, J ) = ZERO
! 420: END DO
! 421: END DO
! 422: *
! 423: * Compute B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS),
! 424: * using the compact WY representation of Q,
! 425: * workspace at least NRHS, optimally NRHS*NB.
! 426: *
! 427: CALL ZGEMQRT( 'Left', 'No transpose', M, NRHS, N, NB,
! 428: $ A, LDA, WORK( 1 ), NB, B, LDB,
! 429: $ WORK( MN*NB+1 ), INFO )
! 430: *
! 431: SCLLEN = M
! 432: *
! 433: END IF
! 434: *
! 435: ELSE
! 436: *
! 437: * M < N:
! 438: * Compute the blocked LQ factorization of A,
! 439: * using the compact WY representation of Q,
! 440: * workspace at least M, optimally M*NB.
! 441: *
! 442: CALL ZGELQT( M, N, NB, A, LDA, WORK( 1 ), NB,
! 443: $ WORK( MN*NB+1 ), INFO )
! 444: *
! 445: IF( .NOT.TPSD ) THEN
! 446: *
! 447: * M < N, A is not transposed:
! 448: * Underdetermined system of equations,
! 449: * minimum norm solution of A * X = B.
! 450: *
! 451: * Compute B := inv(L) * B in two row blocks of B.
! 452: *
! 453: * Block 1: B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
! 454: *
! 455: CALL ZTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS,
! 456: $ A, LDA, B, LDB, INFO )
! 457: *
! 458: IF( INFO.GT.0 ) THEN
! 459: RETURN
! 460: END IF
! 461: *
! 462: * Block 2: Zero out all rows below the M-th row in B:
! 463: * B(M+1:N,1:NRHS) = ZERO
! 464: *
! 465: DO J = 1, NRHS
! 466: DO I = M + 1, N
! 467: B( I, J ) = ZERO
! 468: END DO
! 469: END DO
! 470: *
! 471: * Compute B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS),
! 472: * using the compact WY representation of Q,
! 473: * workspace at least NRHS, optimally NRHS*NB.
! 474: *
! 475: CALL ZGEMLQT( 'Left', 'Conjugate transpose', N, NRHS, M, NB,
! 476: $ A, LDA, WORK( 1 ), NB, B, LDB,
! 477: $ WORK( MN*NB+1 ), INFO )
! 478: *
! 479: SCLLEN = N
! 480: *
! 481: ELSE
! 482: *
! 483: * M < N, A is transposed:
! 484: * Overdetermined system of equations,
! 485: * least-squares problem, min || A**T * X - B ||.
! 486: *
! 487: * Compute B(1:N,1:NRHS) := Q * B(1:N,1:NRHS),
! 488: * using the compact WY representation of Q,
! 489: * workspace at least NRHS, optimally NRHS*NB.
! 490: *
! 491: CALL ZGEMLQT( 'Left', 'No transpose', N, NRHS, M, NB,
! 492: $ A, LDA, WORK( 1 ), NB, B, LDB,
! 493: $ WORK( MN*NB+1), INFO )
! 494: *
! 495: * Compute B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
! 496: *
! 497: CALL ZTRTRS( 'Lower', 'Conjugate transpose', 'Non-unit',
! 498: $ M, NRHS, A, LDA, B, LDB, INFO )
! 499: *
! 500: IF( INFO.GT.0 ) THEN
! 501: RETURN
! 502: END IF
! 503: *
! 504: SCLLEN = M
! 505: *
! 506: END IF
! 507: *
! 508: END IF
! 509: *
! 510: * Undo scaling
! 511: *
! 512: IF( IASCL.EQ.1 ) THEN
! 513: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
! 514: $ INFO )
! 515: ELSE IF( IASCL.EQ.2 ) THEN
! 516: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
! 517: $ INFO )
! 518: END IF
! 519: IF( IBSCL.EQ.1 ) THEN
! 520: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
! 521: $ INFO )
! 522: ELSE IF( IBSCL.EQ.2 ) THEN
! 523: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
! 524: $ INFO )
! 525: END IF
! 526: *
! 527: WORK( 1 ) = DBLE( LWOPT )
! 528: *
! 529: RETURN
! 530: *
! 531: * End of ZGELST
! 532: *
! 533: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgelst.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack