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