Annotation of rpl/lapack/lapack/zgetsls.f, revision 1.1
1.1 ! bertrand 1: * Definition:
! 2: * ===========
! 3: *
! 4: * SUBROUTINE ZGETSLS( 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: * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
! 13: * ..
! 14: *
! 15: *
! 16: *> \par Purpose:
! 17: * =============
! 18: *>
! 19: *> \verbatim
! 20: *>
! 21: *> ZGETSLS solves overdetermined or underdetermined complex 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 = 'C' and m >= n: find the minimum norm solution of
! 37: *> an undetermined system A**T * X = B.
! 38: *>
! 39: *> 4. If TRANS = 'C' 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: *> = 'C': the linear system involves A**C.
! 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 COMPLEX*16 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 ZGEQR or ZGELQ.
! 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 COMPLEX*16 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 = 'C'.
! 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 = 'C' and m >= n, rows 1 to M of B contain the
! 106: *> minimum norm solution vectors;
! 107: *> if TRANS = 'C' 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) COMPLEX*16 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 complex16GEsolve
! 158: *
! 159: * =====================================================================
! 160: SUBROUTINE ZGETSLS( 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: COMPLEX*16 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: COMPLEX*16 CZERO
! 183: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
! 184: * ..
! 185: * .. Local Scalars ..
! 186: LOGICAL LQUERY, TRAN
! 187: INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
! 188: $ SCLLEN, MNK, TSZO, TSZM, LWO, LWM, LW1, LW2,
! 189: $ WSIZEO, WSIZEM, INFO2
! 190: DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
! 191: COMPLEX*16 TQ( 5 ), WORKQ
! 192: * ..
! 193: * .. External Functions ..
! 194: LOGICAL LSAME
! 195: INTEGER ILAENV
! 196: DOUBLE PRECISION DLAMCH, ZLANGE
! 197: EXTERNAL LSAME, ILAENV, DLABAD, DLAMCH, ZLANGE
! 198: * ..
! 199: * .. External Subroutines ..
! 200: EXTERNAL ZGEQR, ZGEMQR, ZLASCL, ZLASET,
! 201: $ ZTRTRS, XERBLA, ZGELQ, ZGEMLQ
! 202: * ..
! 203: * .. Intrinsic Functions ..
! 204: INTRINSIC DBLE, MAX, MIN, INT
! 205: * ..
! 206: * .. Executable Statements ..
! 207: *
! 208: * Test the input arguments.
! 209: *
! 210: INFO = 0
! 211: MINMN = MIN( M, N )
! 212: MAXMN = MAX( M, N )
! 213: MNK = MAX( MINMN, NRHS )
! 214: TRAN = LSAME( TRANS, 'C' )
! 215: *
! 216: LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
! 217: IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
! 218: $ LSAME( TRANS, 'C' ) ) ) THEN
! 219: INFO = -1
! 220: ELSE IF( M.LT.0 ) THEN
! 221: INFO = -2
! 222: ELSE IF( N.LT.0 ) THEN
! 223: INFO = -3
! 224: ELSE IF( NRHS.LT.0 ) THEN
! 225: INFO = -4
! 226: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 227: INFO = -6
! 228: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
! 229: INFO = -8
! 230: END IF
! 231: *
! 232: IF( INFO.EQ.0 ) THEN
! 233: *
! 234: * Determine the block size and minimum LWORK
! 235: *
! 236: IF( M.GE.N ) THEN
! 237: CALL ZGEQR( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
! 238: TSZO = INT( TQ( 1 ) )
! 239: LWO = INT( WORKQ )
! 240: CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
! 241: $ TSZO, B, LDB, WORKQ, -1, INFO2 )
! 242: LWO = MAX( LWO, INT( WORKQ ) )
! 243: CALL ZGEQR( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
! 244: TSZM = INT( TQ( 1 ) )
! 245: LWM = INT( WORKQ )
! 246: CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
! 247: $ TSZM, B, LDB, WORKQ, -1, INFO2 )
! 248: LWM = MAX( LWM, INT( WORKQ ) )
! 249: WSIZEO = TSZO + LWO
! 250: WSIZEM = TSZM + LWM
! 251: ELSE
! 252: CALL ZGELQ( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
! 253: TSZO = INT( TQ( 1 ) )
! 254: LWO = INT( WORKQ )
! 255: CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
! 256: $ TSZO, B, LDB, WORKQ, -1, INFO2 )
! 257: LWO = MAX( LWO, INT( WORKQ ) )
! 258: CALL ZGELQ( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
! 259: TSZM = INT( TQ( 1 ) )
! 260: LWM = INT( WORKQ )
! 261: CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
! 262: $ TSZO, B, LDB, WORKQ, -1, INFO2 )
! 263: LWM = MAX( LWM, INT( WORKQ ) )
! 264: WSIZEO = TSZO + LWO
! 265: WSIZEM = TSZM + LWM
! 266: END IF
! 267: *
! 268: IF( ( LWORK.LT.WSIZEM ).AND.( .NOT.LQUERY ) ) THEN
! 269: INFO = -10
! 270: END IF
! 271: *
! 272: END IF
! 273: *
! 274: IF( INFO.NE.0 ) THEN
! 275: CALL XERBLA( 'ZGETSLS', -INFO )
! 276: WORK( 1 ) = DBLE( WSIZEO )
! 277: RETURN
! 278: END IF
! 279: IF( LQUERY ) THEN
! 280: IF( LWORK.EQ.-1 ) WORK( 1 ) = REAL( WSIZEO )
! 281: IF( LWORK.EQ.-2 ) WORK( 1 ) = REAL( WSIZEM )
! 282: RETURN
! 283: END IF
! 284: IF( LWORK.LT.WSIZEO ) THEN
! 285: LW1 = TSZM
! 286: LW2 = LWM
! 287: ELSE
! 288: LW1 = TSZO
! 289: LW2 = LWO
! 290: END IF
! 291: *
! 292: * Quick return if possible
! 293: *
! 294: IF( MIN( M, N, NRHS ).EQ.0 ) THEN
! 295: CALL ZLASET( 'FULL', MAX( M, N ), NRHS, CZERO, CZERO,
! 296: $ B, LDB )
! 297: RETURN
! 298: END IF
! 299: *
! 300: * Get machine parameters
! 301: *
! 302: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
! 303: BIGNUM = ONE / SMLNUM
! 304: CALL DLABAD( SMLNUM, BIGNUM )
! 305: *
! 306: * Scale A, B if max element outside range [SMLNUM,BIGNUM]
! 307: *
! 308: ANRM = ZLANGE( 'M', M, N, A, LDA, WORK )
! 309: IASCL = 0
! 310: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 311: *
! 312: * Scale matrix norm up to SMLNUM
! 313: *
! 314: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
! 315: IASCL = 1
! 316: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 317: *
! 318: * Scale matrix norm down to BIGNUM
! 319: *
! 320: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
! 321: IASCL = 2
! 322: ELSE IF( ANRM.EQ.ZERO ) THEN
! 323: *
! 324: * Matrix all zero. Return zero solution.
! 325: *
! 326: CALL ZLASET( 'F', MAXMN, NRHS, CZERO, CZERO, B, LDB )
! 327: GO TO 50
! 328: END IF
! 329: *
! 330: BROW = M
! 331: IF ( TRAN ) THEN
! 332: BROW = N
! 333: END IF
! 334: BNRM = ZLANGE( 'M', BROW, NRHS, B, LDB, WORK )
! 335: IBSCL = 0
! 336: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
! 337: *
! 338: * Scale matrix norm up to SMLNUM
! 339: *
! 340: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
! 341: $ INFO )
! 342: IBSCL = 1
! 343: ELSE IF( BNRM.GT.BIGNUM ) THEN
! 344: *
! 345: * Scale matrix norm down to BIGNUM
! 346: *
! 347: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
! 348: $ INFO )
! 349: IBSCL = 2
! 350: END IF
! 351: *
! 352: IF ( M.GE.N ) THEN
! 353: *
! 354: * compute QR factorization of A
! 355: *
! 356: CALL ZGEQR( M, N, A, LDA, WORK( LW2+1 ), LW1,
! 357: $ WORK( 1 ), LW2, INFO )
! 358: IF ( .NOT.TRAN ) THEN
! 359: *
! 360: * Least-Squares Problem min || A * X - B ||
! 361: *
! 362: * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
! 363: *
! 364: CALL ZGEMQR( 'L' , 'C', M, NRHS, N, A, LDA,
! 365: $ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
! 366: $ INFO )
! 367: *
! 368: * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
! 369: *
! 370: CALL ZTRTRS( 'U', 'N', 'N', N, NRHS,
! 371: $ A, LDA, B, LDB, INFO )
! 372: IF( INFO.GT.0 ) THEN
! 373: RETURN
! 374: END IF
! 375: SCLLEN = N
! 376: ELSE
! 377: *
! 378: * Overdetermined system of equations A**T * X = B
! 379: *
! 380: * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
! 381: *
! 382: CALL ZTRTRS( 'U', 'C', 'N', N, NRHS,
! 383: $ A, LDA, B, LDB, INFO )
! 384: *
! 385: IF( INFO.GT.0 ) THEN
! 386: RETURN
! 387: END IF
! 388: *
! 389: * B(N+1:M,1:NRHS) = CZERO
! 390: *
! 391: DO 20 J = 1, NRHS
! 392: DO 10 I = N + 1, M
! 393: B( I, J ) = CZERO
! 394: 10 CONTINUE
! 395: 20 CONTINUE
! 396: *
! 397: * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
! 398: *
! 399: CALL ZGEMQR( 'L', 'N', M, NRHS, N, A, LDA,
! 400: $ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
! 401: $ INFO )
! 402: *
! 403: SCLLEN = M
! 404: *
! 405: END IF
! 406: *
! 407: ELSE
! 408: *
! 409: * Compute LQ factorization of A
! 410: *
! 411: CALL ZGELQ( M, N, A, LDA, WORK( LW2+1 ), LW1,
! 412: $ WORK( 1 ), LW2, INFO )
! 413: *
! 414: * workspace at least M, optimally M*NB.
! 415: *
! 416: IF( .NOT.TRAN ) THEN
! 417: *
! 418: * underdetermined system of equations A * X = B
! 419: *
! 420: * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
! 421: *
! 422: CALL ZTRTRS( 'L', 'N', 'N', M, NRHS,
! 423: $ A, LDA, B, LDB, INFO )
! 424: *
! 425: IF( INFO.GT.0 ) THEN
! 426: RETURN
! 427: END IF
! 428: *
! 429: * B(M+1:N,1:NRHS) = 0
! 430: *
! 431: DO 40 J = 1, NRHS
! 432: DO 30 I = M + 1, N
! 433: B( I, J ) = CZERO
! 434: 30 CONTINUE
! 435: 40 CONTINUE
! 436: *
! 437: * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
! 438: *
! 439: CALL ZGEMLQ( 'L', 'C', N, NRHS, M, A, LDA,
! 440: $ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
! 441: $ INFO )
! 442: *
! 443: * workspace at least NRHS, optimally NRHS*NB
! 444: *
! 445: SCLLEN = N
! 446: *
! 447: ELSE
! 448: *
! 449: * overdetermined system min || A**T * X - B ||
! 450: *
! 451: * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
! 452: *
! 453: CALL ZGEMLQ( 'L', 'N', N, NRHS, M, A, LDA,
! 454: $ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
! 455: $ INFO )
! 456: *
! 457: * workspace at least NRHS, optimally NRHS*NB
! 458: *
! 459: * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
! 460: *
! 461: CALL ZTRTRS( 'L', 'C', 'N', M, NRHS,
! 462: $ A, LDA, B, LDB, INFO )
! 463: *
! 464: IF( INFO.GT.0 ) THEN
! 465: RETURN
! 466: END IF
! 467: *
! 468: SCLLEN = M
! 469: *
! 470: END IF
! 471: *
! 472: END IF
! 473: *
! 474: * Undo scaling
! 475: *
! 476: IF( IASCL.EQ.1 ) THEN
! 477: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
! 478: $ INFO )
! 479: ELSE IF( IASCL.EQ.2 ) THEN
! 480: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
! 481: $ INFO )
! 482: END IF
! 483: IF( IBSCL.EQ.1 ) THEN
! 484: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
! 485: $ INFO )
! 486: ELSE IF( IBSCL.EQ.2 ) THEN
! 487: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
! 488: $ INFO )
! 489: END IF
! 490: *
! 491: 50 CONTINUE
! 492: WORK( 1 ) = DBLE( TSZO + LWO )
! 493: RETURN
! 494: *
! 495: * End of ZGETSLS
! 496: *
! 497: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgetsls.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack