Annotation of rpl/lapack/lapack/dorm22.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DORM22 multiplies a general matrix by a banded orthogonal matrix.
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DORM22 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorm22.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorm22.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorm22.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
! 22: * $ WORK, LWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER SIDE, TRANS
! 26: * INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * DOUBLE PRECISION Q( LDQ, * ), C( LDC, * ), WORK( * )
! 30: * ..
! 31: *
! 32: *> \par Purpose
! 33: * ============
! 34: *>
! 35: *> \verbatim
! 36: *>
! 37: *>
! 38: *> DORM22 overwrites the general real M-by-N matrix C with
! 39: *>
! 40: *> SIDE = 'L' SIDE = 'R'
! 41: *> TRANS = 'N': Q * C C * Q
! 42: *> TRANS = 'T': Q**T * C C * Q**T
! 43: *>
! 44: *> where Q is a real orthogonal matrix of order NQ, with NQ = M if
! 45: *> SIDE = 'L' and NQ = N if SIDE = 'R'.
! 46: *> The orthogonal matrix Q processes a 2-by-2 block structure
! 47: *>
! 48: *> [ Q11 Q12 ]
! 49: *> Q = [ ]
! 50: *> [ Q21 Q22 ],
! 51: *>
! 52: *> where Q12 is an N1-by-N1 lower triangular matrix and Q21 is an
! 53: *> N2-by-N2 upper triangular matrix.
! 54: *> \endverbatim
! 55: *
! 56: * Arguments
! 57: * =========
! 58: *
! 59: *> \param[in] SIDE
! 60: *> \verbatim
! 61: *> SIDE is CHARACTER*1
! 62: *> = 'L': apply Q or Q**T from the Left;
! 63: *> = 'R': apply Q or Q**T from the Right.
! 64: *> \endverbatim
! 65: *>
! 66: *> \param[in] TRANS
! 67: *> \verbatim
! 68: *> TRANS is CHARACTER*1
! 69: *> = 'N': apply Q (No transpose);
! 70: *> = 'C': apply Q**T (Conjugate transpose).
! 71: *> \endverbatim
! 72: *>
! 73: *> \param[in] M
! 74: *> \verbatim
! 75: *> M is INTEGER
! 76: *> The number of rows of the matrix C. M >= 0.
! 77: *> \endverbatim
! 78: *>
! 79: *> \param[in] N
! 80: *> \verbatim
! 81: *> N is INTEGER
! 82: *> The number of columns of the matrix C. N >= 0.
! 83: *> \endverbatim
! 84: *>
! 85: *> \param[in] N1
! 86: *> \param[in] N2
! 87: *> \verbatim
! 88: *> N1 is INTEGER
! 89: *> N2 is INTEGER
! 90: *> The dimension of Q12 and Q21, respectively. N1, N2 >= 0.
! 91: *> The following requirement must be satisfied:
! 92: *> N1 + N2 = M if SIDE = 'L' and N1 + N2 = N if SIDE = 'R'.
! 93: *> \endverbatim
! 94: *>
! 95: *> \param[in] Q
! 96: *> \verbatim
! 97: *> Q is DOUBLE PRECISION array, dimension
! 98: *> (LDQ,M) if SIDE = 'L'
! 99: *> (LDQ,N) if SIDE = 'R'
! 100: *> \endverbatim
! 101: *>
! 102: *> \param[in] LDQ
! 103: *> \verbatim
! 104: *> LDQ is INTEGER
! 105: *> The leading dimension of the array Q.
! 106: *> LDQ >= max(1,M) if SIDE = 'L'; LDQ >= max(1,N) if SIDE = 'R'.
! 107: *> \endverbatim
! 108: *>
! 109: *> \param[in,out] C
! 110: *> \verbatim
! 111: *> C is DOUBLE PRECISION array, dimension (LDC,N)
! 112: *> On entry, the M-by-N matrix C.
! 113: *> On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
! 114: *> \endverbatim
! 115: *>
! 116: *> \param[in] LDC
! 117: *> \verbatim
! 118: *> LDC is INTEGER
! 119: *> The leading dimension of the array C. LDC >= max(1,M).
! 120: *> \endverbatim
! 121: *>
! 122: *> \param[out] WORK
! 123: *> \verbatim
! 124: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 125: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[in] LWORK
! 129: *> \verbatim
! 130: *> LWORK is INTEGER
! 131: *> The dimension of the array WORK.
! 132: *> If SIDE = 'L', LWORK >= max(1,N);
! 133: *> if SIDE = 'R', LWORK >= max(1,M).
! 134: *> For optimum performance LWORK >= M*N.
! 135: *>
! 136: *> If LWORK = -1, then a workspace query is assumed; the routine
! 137: *> only calculates the optimal size of the WORK array, returns
! 138: *> this value as the first entry of the WORK array, and no error
! 139: *> message related to LWORK is issued by XERBLA.
! 140: *> \endverbatim
! 141: *>
! 142: *> \param[out] INFO
! 143: *> \verbatim
! 144: *> INFO is INTEGER
! 145: *> = 0: successful exit
! 146: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 147: *> \endverbatim
! 148: *
! 149: *
! 150: * Authors:
! 151: * ========
! 152: *
! 153: *> \author Univ. of Tennessee
! 154: *> \author Univ. of California Berkeley
! 155: *> \author Univ. of Colorado Denver
! 156: *> \author NAG Ltd.
! 157: *
! 158: *> \date January 2015
! 159: *
! 160: *> \ingroup complexOTHERcomputational
! 161: *
! 162: * =====================================================================
! 163: SUBROUTINE DORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
! 164: $ WORK, LWORK, INFO )
! 165: *
! 166: * -- LAPACK computational routine (version 3.6.0) --
! 167: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 168: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 169: * January 2015
! 170: *
! 171: IMPLICIT NONE
! 172: *
! 173: * .. Scalar Arguments ..
! 174: CHARACTER SIDE, TRANS
! 175: INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
! 176: * ..
! 177: * .. Array Arguments ..
! 178: DOUBLE PRECISION Q( LDQ, * ), C( LDC, * ), WORK( * )
! 179: * ..
! 180: *
! 181: * =====================================================================
! 182: *
! 183: * .. Parameters ..
! 184: DOUBLE PRECISION ONE
! 185: PARAMETER ( ONE = 1.0D+0 )
! 186: *
! 187: * .. Local Scalars ..
! 188: LOGICAL LEFT, LQUERY, NOTRAN
! 189: INTEGER I, LDWORK, LEN, LWKOPT, NB, NQ, NW
! 190: * ..
! 191: * .. External Functions ..
! 192: LOGICAL LSAME
! 193: EXTERNAL LSAME
! 194: * ..
! 195: * .. External Subroutines ..
! 196: EXTERNAL DGEMM, DLACPY, DTRMM, XERBLA
! 197: * ..
! 198: * .. Intrinsic Functions ..
! 199: INTRINSIC DBLE, MAX, MIN
! 200: * ..
! 201: * .. Executable Statements ..
! 202: *
! 203: * Test the input arguments
! 204: *
! 205: INFO = 0
! 206: LEFT = LSAME( SIDE, 'L' )
! 207: NOTRAN = LSAME( TRANS, 'N' )
! 208: LQUERY = ( LWORK.EQ.-1 )
! 209: *
! 210: * NQ is the order of Q;
! 211: * NW is the minimum dimension of WORK.
! 212: *
! 213: IF( LEFT ) THEN
! 214: NQ = M
! 215: ELSE
! 216: NQ = N
! 217: END IF
! 218: NW = NQ
! 219: IF( N1.EQ.0 .OR. N2.EQ.0 ) NW = 1
! 220: IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
! 221: INFO = -1
! 222: ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'T' ) )
! 223: $ THEN
! 224: INFO = -2
! 225: ELSE IF( M.LT.0 ) THEN
! 226: INFO = -3
! 227: ELSE IF( N.LT.0 ) THEN
! 228: INFO = -4
! 229: ELSE IF( N1.LT.0 .OR. N1+N2.NE.NQ ) THEN
! 230: INFO = -5
! 231: ELSE IF( N2.LT.0 ) THEN
! 232: INFO = -6
! 233: ELSE IF( LDQ.LT.MAX( 1, NQ ) ) THEN
! 234: INFO = -8
! 235: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
! 236: INFO = -10
! 237: ELSE IF( LWORK.LT.NW .AND. .NOT.LQUERY ) THEN
! 238: INFO = -12
! 239: END IF
! 240: *
! 241: IF( INFO.EQ.0 ) THEN
! 242: LWKOPT = M*N
! 243: WORK( 1 ) = DBLE( LWKOPT )
! 244: END IF
! 245: *
! 246: IF( INFO.NE.0 ) THEN
! 247: CALL XERBLA( 'DORM22', -INFO )
! 248: RETURN
! 249: ELSE IF( LQUERY ) THEN
! 250: RETURN
! 251: END IF
! 252: *
! 253: * Quick return if possible
! 254: *
! 255: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
! 256: WORK( 1 ) = 1
! 257: RETURN
! 258: END IF
! 259: *
! 260: * Degenerate cases (N1 = 0 or N2 = 0) are handled using DTRMM.
! 261: *
! 262: IF( N1.EQ.0 ) THEN
! 263: CALL DTRMM( SIDE, 'Upper', TRANS, 'Non-Unit', M, N, ONE,
! 264: $ Q, LDQ, C, LDC )
! 265: WORK( 1 ) = ONE
! 266: RETURN
! 267: ELSE IF( N2.EQ.0 ) THEN
! 268: CALL DTRMM( SIDE, 'Lower', TRANS, 'Non-Unit', M, N, ONE,
! 269: $ Q, LDQ, C, LDC )
! 270: WORK( 1 ) = ONE
! 271: RETURN
! 272: END IF
! 273: *
! 274: * Compute the largest chunk size available from the workspace.
! 275: *
! 276: NB = MAX( 1, MIN( LWORK, LWKOPT ) / NQ )
! 277: *
! 278: IF( LEFT ) THEN
! 279: IF( NOTRAN ) THEN
! 280: DO I = 1, N, NB
! 281: LEN = MIN( NB, N-I+1 )
! 282: LDWORK = M
! 283: *
! 284: * Multiply bottom part of C by Q12.
! 285: *
! 286: CALL DLACPY( 'All', N1, LEN, C( N2+1, I ), LDC, WORK,
! 287: $ LDWORK )
! 288: CALL DTRMM( 'Left', 'Lower', 'No Transpose', 'Non-Unit',
! 289: $ N1, LEN, ONE, Q( 1, N2+1 ), LDQ, WORK,
! 290: $ LDWORK )
! 291: *
! 292: * Multiply top part of C by Q11.
! 293: *
! 294: CALL DGEMM( 'No Transpose', 'No Transpose', N1, LEN, N2,
! 295: $ ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
! 296: $ LDWORK )
! 297: *
! 298: * Multiply top part of C by Q21.
! 299: *
! 300: CALL DLACPY( 'All', N2, LEN, C( 1, I ), LDC,
! 301: $ WORK( N1+1 ), LDWORK )
! 302: CALL DTRMM( 'Left', 'Upper', 'No Transpose', 'Non-Unit',
! 303: $ N2, LEN, ONE, Q( N1+1, 1 ), LDQ,
! 304: $ WORK( N1+1 ), LDWORK )
! 305: *
! 306: * Multiply bottom part of C by Q22.
! 307: *
! 308: CALL DGEMM( 'No Transpose', 'No Transpose', N2, LEN, N1,
! 309: $ ONE, Q( N1+1, N2+1 ), LDQ, C( N2+1, I ), LDC,
! 310: $ ONE, WORK( N1+1 ), LDWORK )
! 311: *
! 312: * Copy everything back.
! 313: *
! 314: CALL DLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
! 315: $ LDC )
! 316: END DO
! 317: ELSE
! 318: DO I = 1, N, NB
! 319: LEN = MIN( NB, N-I+1 )
! 320: LDWORK = M
! 321: *
! 322: * Multiply bottom part of C by Q21**T.
! 323: *
! 324: CALL DLACPY( 'All', N2, LEN, C( N1+1, I ), LDC, WORK,
! 325: $ LDWORK )
! 326: CALL DTRMM( 'Left', 'Upper', 'Transpose', 'Non-Unit',
! 327: $ N2, LEN, ONE, Q( N1+1, 1 ), LDQ, WORK,
! 328: $ LDWORK )
! 329: *
! 330: * Multiply top part of C by Q11**T.
! 331: *
! 332: CALL DGEMM( 'Transpose', 'No Transpose', N2, LEN, N1,
! 333: $ ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
! 334: $ LDWORK )
! 335: *
! 336: * Multiply top part of C by Q12**T.
! 337: *
! 338: CALL DLACPY( 'All', N1, LEN, C( 1, I ), LDC,
! 339: $ WORK( N2+1 ), LDWORK )
! 340: CALL DTRMM( 'Left', 'Lower', 'Transpose', 'Non-Unit',
! 341: $ N1, LEN, ONE, Q( 1, N2+1 ), LDQ,
! 342: $ WORK( N2+1 ), LDWORK )
! 343: *
! 344: * Multiply bottom part of C by Q22**T.
! 345: *
! 346: CALL DGEMM( 'Transpose', 'No Transpose', N1, LEN, N2,
! 347: $ ONE, Q( N1+1, N2+1 ), LDQ, C( N1+1, I ), LDC,
! 348: $ ONE, WORK( N2+1 ), LDWORK )
! 349: *
! 350: * Copy everything back.
! 351: *
! 352: CALL DLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
! 353: $ LDC )
! 354: END DO
! 355: END IF
! 356: ELSE
! 357: IF( NOTRAN ) THEN
! 358: DO I = 1, M, NB
! 359: LEN = MIN( NB, M-I+1 )
! 360: LDWORK = LEN
! 361: *
! 362: * Multiply right part of C by Q21.
! 363: *
! 364: CALL DLACPY( 'All', LEN, N2, C( I, N1+1 ), LDC, WORK,
! 365: $ LDWORK )
! 366: CALL DTRMM( 'Right', 'Upper', 'No Transpose', 'Non-Unit',
! 367: $ LEN, N2, ONE, Q( N1+1, 1 ), LDQ, WORK,
! 368: $ LDWORK )
! 369: *
! 370: * Multiply left part of C by Q11.
! 371: *
! 372: CALL DGEMM( 'No Transpose', 'No Transpose', LEN, N2, N1,
! 373: $ ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
! 374: $ LDWORK )
! 375: *
! 376: * Multiply left part of C by Q12.
! 377: *
! 378: CALL DLACPY( 'All', LEN, N1, C( I, 1 ), LDC,
! 379: $ WORK( 1 + N2*LDWORK ), LDWORK )
! 380: CALL DTRMM( 'Right', 'Lower', 'No Transpose', 'Non-Unit',
! 381: $ LEN, N1, ONE, Q( 1, N2+1 ), LDQ,
! 382: $ WORK( 1 + N2*LDWORK ), LDWORK )
! 383: *
! 384: * Multiply right part of C by Q22.
! 385: *
! 386: CALL DGEMM( 'No Transpose', 'No Transpose', LEN, N1, N2,
! 387: $ ONE, C( I, N1+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
! 388: $ ONE, WORK( 1 + N2*LDWORK ), LDWORK )
! 389: *
! 390: * Copy everything back.
! 391: *
! 392: CALL DLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
! 393: $ LDC )
! 394: END DO
! 395: ELSE
! 396: DO I = 1, M, NB
! 397: LEN = MIN( NB, M-I+1 )
! 398: LDWORK = LEN
! 399: *
! 400: * Multiply right part of C by Q12**T.
! 401: *
! 402: CALL DLACPY( 'All', LEN, N1, C( I, N2+1 ), LDC, WORK,
! 403: $ LDWORK )
! 404: CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Non-Unit',
! 405: $ LEN, N1, ONE, Q( 1, N2+1 ), LDQ, WORK,
! 406: $ LDWORK )
! 407: *
! 408: * Multiply left part of C by Q11**T.
! 409: *
! 410: CALL DGEMM( 'No Transpose', 'Transpose', LEN, N1, N2,
! 411: $ ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
! 412: $ LDWORK )
! 413: *
! 414: * Multiply left part of C by Q21**T.
! 415: *
! 416: CALL DLACPY( 'All', LEN, N2, C( I, 1 ), LDC,
! 417: $ WORK( 1 + N1*LDWORK ), LDWORK )
! 418: CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-Unit',
! 419: $ LEN, N2, ONE, Q( N1+1, 1 ), LDQ,
! 420: $ WORK( 1 + N1*LDWORK ), LDWORK )
! 421: *
! 422: * Multiply right part of C by Q22**T.
! 423: *
! 424: CALL DGEMM( 'No Transpose', 'Transpose', LEN, N2, N1,
! 425: $ ONE, C( I, N2+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
! 426: $ ONE, WORK( 1 + N1*LDWORK ), LDWORK )
! 427: *
! 428: * Copy everything back.
! 429: *
! 430: CALL DLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
! 431: $ LDC )
! 432: END DO
! 433: END IF
! 434: END IF
! 435: *
! 436: WORK( 1 ) = DBLE( LWKOPT )
! 437: RETURN
! 438: *
! 439: * End of DORM22
! 440: *
! 441: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dorm22.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack